8.3. Arithmétique¶
8.3.1. Décomposition d’un entier dans une base¶
L’écriture d’un entier \(n\in\mathbb{N}^*\) dans une base \(b\) où \(b\) est un entier supérieur ou égal à \(2\) est une écriture de la forme
où les \(a_k\) sont des entiers compris entre \(0\) et \(b-1\) et où \(a_p\neq0\). Par exemple, la décomposition en base \(8\) de \(4621\) est
On peut obtenir la liste \([a_0,a_1,\dots,a_p]\) à l’aide d’une suite de divisions euclidiennes par \(b\). Par exemple,
donc
In [1]: def decomp(n, b):
...: l = []
...: while n != 0:
...: l.append(n % b)
...: n //= b
...: return l
...:
In [2]: decomp(24189, 10 )
Out[2]: [9, 8, 1, 4, 2]
In [3]: n, b = 3678, 7
In [4]: l = decomp(n, b)
In [5]: l
Out[5]: [3, 0, 5, 3, 1]
# On vérifie que la liste convient effectivement
In [6]: sum(a * b ** i for i, a in enumerate(l))
Out[6]: 3678
8.3.2. Calcul de PGCD¶
On peut implémenter l’algorithme d’Euclide en Python.
In [7]: def pgcd(a, b):
...: while b != 0:
...: a, b = b, a % b
...: return abs(a) # Le pgcd doit être positif
...:
In [8]: pgcd(30, 12)
Out[8]: 6
In [9]: pgcd(30, -12)
Out[9]: 6
De même, on peut implémenter l’algorithme d’Euclide étendu qui, en plus du pgcd, donne des coefficients d’une relation de Bézout, c’est-à-dire des entiers \(u\) et \(v\) tels que
In [10]: def bezout(a, b):
....: s, t, u, v = 1, 0, 0, 1
....: while b != 0:
....: q = a // b
....: a, s, t, b, u, v = b, u, v, a - q * b, s - q * u, t - q * v
....: return (a, s, t) if a > 0 else (-a, -s, -t) # Le pgcd doit être positif
....:
In [11]: bezout(30, 12)
Out[11]: (6, 1, -2)
In [12]: bezout(30, -12)
Out[12]: (6, -1, -3)
La preuve de l’algorithme est donnée en annexe.
8.3.3. Exponentiation rapide¶
Il s’agit ici de calculer efficacement une puissance entière d’un objet mathématique (nombre ou matrice par exemple). Un algorithme naïf serait le suivant.
In [13]: def exponentiation(x, n):
....: a = 1
....: for _ in range(n):
....: a *= x
....: return a
....:
In [14]: exponentiation(3, 5)
Out[14]: 243
Il est clair que cet algorithme nécessite \(n\) multiplications pour calculer une puissance \(n^\text{ème}\).
On peut proposer un léger raffinement pour éviter une multiplication.
In [15]: def exponentiation2(x, n):
....: if n == 0:
....: return 1
....: a = x
....: for _ in range(n-1):
....: a *= x
....: return a
....:
In [16]: exponentiation2(3, 5)
Out[16]: 243
Mais on peut être beaucoup plus efficace. On remarque que toute puissance \(x^n\) peut en fait s’écrire comme un produit de puissances de la forme \(x^{2^k}\) : en effet, \(n\) peut s’écrire comme une somme d’entiers de la forme \(2^k\).
Il suffit alors de calculer successivement les \(x^{2^k}\) par des élevations au carré puis de multiplier ces puissances entre elles. Par exemple, dans l’exemple précédent, on utilise :
3 multiplications pour calculer successivement \(x^2\), \(x^4\) et \(x^8\) ;
2 multiplications pour effectuer le produit de \(x\), \(x^4\) et \(x^8\).
On effectue en tout 5 multiplications au lieu de 12.
Pour obtenir la décomposition de \(x^n\) en produit de facteurs de la forme \(x^{2^k}\), il suffit de décomposer \(n\) en base \(2\). On utilise donc un algorithme similaire à l’algorithme de décomposition binaire. En effet, en notant \(q\) et \(r\) le quotient et le reste de la division euclidienne de \(n\) par \(2\)
In [17]: def exponentiation_rapide(x, n):
....: a = 1
....: p = x
....: while n > 0:
....: if n % 2 == 1:
....: a *= p
....: p *= p
....: n //= 2
....: return a
....:
In [18]: exponentiation_rapide(3, 5)
Out[18]: 243
On peut également proposer un raffinement pour gagner une multiplication. En effet, à la dernière itération, l’instruction p *= p est inutile puisque cette dernière valeur de p ne sera pas utilisée. De plus, à la dernière itération, n vaut 1 qui est impair donc l’instruction a *= p sera obligatoirement effectuée.
In [19]: def exponentiation_rapide2(x, n):
....: if n == 0:
....: return 1
....: a = 1
....: p = x
....: while n > 1:
....: if n % 2 == 1:
....: a *= p
....: p *= p
....: n //= 2
....: return a * p
....:
In [20]: exponentiation_rapide2(3, 5)
Out[20]: 243
8.3.4. Evaluation de polynômes¶
La méthode naïve pour évaluer un polynôme \(P\) en un scalaire \(x\) consiste à calculer les différentes puissances de \(x\) puis à les multiplier par les coefficients de \(P\) correspondant puis à effectuer la somme de ces produits.
Par exemple, pour évaluer \(5X^3+4X^2-3X+7\) en un scalaire \(x\), on calcule successivement :
les puissances de \(x\), à savoir \(x^2\) et \(x^3\) (2 multiplications) ;
les produits \(-3x\), \(4x^2\) et \(5x^3\) (3 multiplications) ;
la somme de \(7\), \(-3x\), \(4x^2\) et \(5x^3\) (3 additions).
Mais les calculs peuvent être menés plus astucieusement en remarquant que :
On calcule alors successivement :
\(s_1=4+5x\) (1 multiplication et 1 addition) ;
\(s_2=-3+xs_1\) (1 multiplication et 1 addition) ;
\(s_3=7+xs_2\) (1 multiplication et 1 addition).
On a gagné deux multiplications par rapport à la méthode précédente et on comprend bien que le gain sera d’autant plus grand que le polynôme est de degré élevé.
L’algorithme décrit dans l’exemple précédent s’appelle la méthode de Hörner. Pour implémenter cet algorithme, on représente un polynôme par la liste de ses coefficients rangés par ordre décroissant de degré. Par exemple, la liste [1, 2, 3] représente le polynôme \(X^2+2X+3\).
In [21]: def horner(poly, x):
....: s = 0
....: for c in poly:
....: s = s * x + c
....: return s
....:
# On évalue le polynôme X²+2X+3 en 4
In [22]: horner([1, 2, 3], 4)
Out[22]: 27
Si l’on préfère représenter un polynôme par la liste de ses coefficients par ordre de degré croissant, on peut toujours utiliser la fonction reversed qui fait ce que son nom indique [1].
In [23]: def horner(poly, x):
....: s = 0
....: for c in reversed(poly):
....: s = s * x + c
....: return s
....:
# On évalue le polynôme 1+2X+3X² en 4
In [24]: horner([1, 2, 3], 4)
Out[24]: 57