I. IntroductionDans un précédent billet, on a parlé des avantages du « diviser pour régner » dans le développement d'un produit de petits polynômes, on souhaite maintenant montrer comment réaliser une exponentiation rapide de nombres réels et de polynômes.
Pour cela, on va d'abord décrire cette méthode à l'aide d'exemples simples, pour ensuite l'implémenter en Python.
II. Principe de l'exponentiation rapide
II-A. Nombres réels
Si on souhaite calculer an, on peut naïvement effectuer n−1 multiplications en réalisant le produit :
an = a × a × ... × a
Cependant, l'exemple suivant montre que l'on peut calculer a32 en effectuant seulement 5 multiplications au lieu de 31 :
a32 = a25 = ((((a2)2)2)2)2
Si l'exposant est différent d'une puissance de 2, avec par exemple a41, on peut alors écrire :
a41 = a × a8 × a32 = a × a23 × a25
En remarquant que a8 est également calculé dans a32, et en conservant sa valeur, on a en fait besoin de seulement 7 multiplications au lieu de 40.
Alors que l'algorithme naïf demande de l'ordre de n multiplications pour calculer an, l'algorithme d'exponentiation rapide a besoin seulement de l'ordre de log2(n) multiplications.
II-B. Polynômes
Une multiplication entre deux polynômes nécessite de multiplier chaque terme du multiplicateur par chaque terme du multiplicande.
Si ces deux polynômes ont n termes, cela exige n2 produits. La multiplication (X + 2) × (X + 2) nécessite ainsi 4 multiplications entre termes.
Si on prend maintenant :
(X + 2)8 = (X + 2) × (X + 2) × (X + 2) × (X + 2) × (X + 2) × (X + 2) × (X + 2) × (X + 2)
En développant le produit de gauche à droite, on obtient successivement :
(X + 2)8 = ((X + 2) × (X + 2)) × (X + 2) × (X + 2) × (X + 2) × (X + 2) × (X + 2) × (X + 2)
La 1re multiplication nécessite donc 4 opérations de multiplication entre termes.
En poursuivant le développement :
(X + 2)8 = ((X2 + 4X + 4) × (X + 2)) × (X + 2) × (X + 2) × (X + 2) × (X + 2) × (X + 2)
On a ainsi besoin jusque là de 4 + 6 opérations de multiplication.
etc..
On obtient finalement 4 + 6 + 8 + 10 + 12 + 14 + 16 = 70 opérations de multiplication en tout pour ce produit.
Considérons maintenant la même puissance de (X + 2) réarrangée :
(X + 2)8 = (X + 2)23 = (((X + 2)2)2)2
Si on effectue le développement en respectant la règle de priorité des parenthèses, on obtient successivement :
(X + 2)8 = (((X + 2)2)2)2
2×2=4 opérations de multiplication entre termes pour la 1re produit.
Puis :
(X + 2)8 = ((X2 + 4X + 4))2)2
3×3=9 opérations de multiplication pour le 2e produit.
Et enfin :
(X + 2)8 = (X4 + 8X3 + 24X2 + 32X + 16)2
5×5=25 opérations pour le dernier produit.
Ce qui nous fait au total 4 + 9 + 25 = 38 opérations de multiplication entre termes au lieu de 70.
III. Implémentation en Python
III-A. Exponentiation rapide de nombres réels
L'algorithme est décrit sur la page Wikipedia :
Soit n un entier strictement supérieur à 1, supposons que l'on sache calculer, pour chaque réel x, toutes les puissances xk de x, pour tout k, tel que 1 ≤ k < n.
- Si n est pair alors xn = (x2)n/2. Il suffit alors de calculer yn/2 pour y = x2.
- Si n est impair et n > 1, alors xn = x(x2)(n – 1)/2. Il suffit de calculer y(n – 1)/2 pour y = x2 et de multiplier le résultat par x.
Cette remarque nous amène à l'algorithme récursif suivant qui calcule xn pour un entier strictement positif n :
On obtient ainsi la fonction d'exponentiation rapide de nombres réels :
| Code Python : | Sélectionner tout |
1 2 3 4 5 6 7 8 9 10 11 12 13 14 | def puissance(x,n): # exponentiation rapide de x^n # si n=1 alors on renvoie x if n==1: return x else: # sinon # si n est pair if (n % 2 == 0): # appel récursif : puissance(x^2,n/2) return puissance(x**2,n//2) else: # sinon si n est impair # appel récursif : x*puissance(x^2,(n-1)/2) return x*puissance(x**2,n//2) |
Testons maintenant notre fonction avec ces quelques lignes :
| Code Python : | Sélectionner tout |
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 | # nombre réel x = 12 # exposant n = 8 # exponentiation du polynôme xn = puissance(x,n) print("x^n = {0}^{1}".format(x,n)) print() print("x^n = " + str(xn)) print() print("nombre de multiplications de l'ordre de log2({0})".format(n)) |
Le code affiche :
x^n = 12^8
x^n = 429981696
nombre de multiplications de l'ordre de log2(8)
III-B. Puissance de polynômes
On utilise à nouveau notre classe Polynome, comme dans le billet précédent, en ajoutant un attribut permettant de compter le nombre d'opérations de multiplication entre termes nécessaires pour développer le produit de polynômes :
| Code Python : | Sélectionner tout |
1 2 3 4 5 6 7 8 9 10 11 12 13 14 | class Polynome: def __init__(self, liste_coefs=[0], nbre_operations=0): # méthode constructeur de la classe # on définit la liste des coefficients du polynôme [a0, a1, ..., an] self.coefs = liste_coefs # on initialise le nombre d'opérations de multiplication entre termes self.nombre_operations = nbre_operations # suppression si nécessaire des zéros en queue de liste de coefficients. Exemple : [2, 3, 1, 0, 0] -> [2, 3, 1] self.reduire() ... |
Le code de la méthode permettant de multiplier deux polynômes devient alors :
| Code Python : | Sélectionner tout |
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 | class Polynome: ... def __mul__(self, other): # méthode permettant de redéfinir l'opérateur « * » pour 2 polynômes : (2 + X) * (1 + 2x) = 1 + 3x + 2x^2 # initialisation de la liste des coefficients qu'avec des zéros liste_coefs=[0]*(len(self.coefs)+len(other.coefs)-1) # exemple : [0, 0, 0] for i1 in range(len(self.coefs)): # parcours des indices des coefs du polynôme n°1 for i2 in range(len(other.coefs)): # parcours des indices des coefs du polynôme n°2 # multiplication des coefficients d'indices i1 et i2 liste_coefs[i1+i2] = liste_coefs[i1+i2] + self.coefs[i1]*other.coefs[i2] poly = Polynome(liste_coefs) # création de l'objet Polynome basé sur la liste # ajout du nombre d'opérations de multiplication entre termes effectuées lors de la multiplication entre polynômes poly.compteur_operations = self.compteur_operations + other.compteur_operations + len(self.coefs)*len(other.coefs) return poly # renvoie le polynôme résultat de la multiplication |
Rappel important : les objets Polynome de notre classe sont représentés par une liste de coefficients dont la longueur est égale au degré du polynôme + 1:
Par exemple le polynôme X3 + 2 = 2 + 0.X + 0.X2 + X3 sera représenté par [2, 0, 0, 1].
III-B-1. Méthode classique
On doit d'abord réécrire le code de la méthode permettant de développer une puissance de polynôme afin de pouvoir comptabiliser le nombre de multiplications entre termes :
| Code Python : | Sélectionner tout |
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 | class Polynome: ... def __pow__(self, n): # méthode permettant de redéfinir l'opérateur de puissance : self ** n if n==0: # si n=0 return Polynome([1]) # on renvoie le polynôme égal à 1 # copie de self initialisé avec le même nombre d'opérations de multiplication entre termes poly = Polynome(self.coefs, self.nombre_operations) # copie de self avec un nombre d'opérations de multiplication à 0 pol = Polynome(self.coefs) for i in range(n-1): # on multiplie n fois poly par pol à l'aide de l'opérateur * poly = poly*pol # équivalent à : poly = poly.__mul__(pol) return poly # renvoie le polynôme résultat de l'opération (self ** n) |
Testons maintenant l'opérateur de puissance « ** » pour les polynômes avec la méthode classique :
| Code Python : | Sélectionner tout |
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 | # création du polynôme p = 2 + X = X + 2 p = Polynome([2,1]) # exposant n = 8 # exponentiation du polynôme pn = p**n print("p^n = ({0})^{1}".format(p,n)) print() print("p^n = " + str(pn)) print() print("nombre d'opérations = " + str(pn.nombre_operations)) |
Le code affiche :
p^n = (2 + X)^8
p^n = 256 + 1024⋅X + 1792⋅X^2 + 1792⋅X^3 + 1120⋅X^4 + 448⋅X^5 + 112⋅X^6 + 16⋅X^7 + X^8
nombre d'opérations = 70
III-B-2. Exponentiation rapide de polynômes
L'algorithme d'exponentiation rapide de polynômes est le même que celui décrit sur la page Wikipedia pour les nombres réels :
Soit n un entier strictement supérieur à 1, supposons que l'on sache développer, pour chaque polynôme p, toutes les puissances pk de p, pour tout k, tel que 1 ≤ k < n.
- Si n est pair alors pn = (p2)n/2. Il suffit alors de développer qn/2 pour q = p2.
- Si n est impair et n > 1, alors pn = p(p2)(n – 1)/2. Il suffit de développer q(n – 1)/2 pour q = p2 et de multiplier le résultat par p.
Cette remarque nous amène à l'algorithme récursif suivant qui développe pn pour un entier strictement positif n :
On peut alors écrire en Python la fonction d'exponentiation rapide de polynômes :
| Code Python : | Sélectionner tout |
1 2 3 4 5 6 7 8 9 10 11 12 13 14 | def puissance_polynome(p,n): # exponentiation rapide de p^n # si n=1 alors on renvoie p if n==1: return p else: # sinon # si n est pair if (n % 2 == 0): # appel récursif : puissance_polynome(p^2,n/2) return puissance_polynome(p**2,n//2) else: # sinon si n est impair # appel récursif : p*puissance_polynome(p^2,(n-1)/2) return p*puissance_polynome(p**2,n//2) |
Testons notre fonction avec ces quelques lignes :
| Code Python : | Sélectionner tout |
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 | # création du polynôme p = 2 + X = X + 2 p = Polynome([2,1]) # exposant n = 8 # exponentiation du polynôme pn = puissance_polynome(p,n) print("p^n = ({0})^{1}".format(p,n)) print() print("p^n = " + str(pn)) print() print("nombre d'opérations = " + str(pn.nombre_operations)) |
Le code affiche :
p^n = (2 + X)^8
p^n = 256 + 1024⋅X + 1792⋅X^2 + 1792⋅X^3 + 1120⋅X^4 + 448⋅X^5 + 112⋅X^6 + 16⋅X^7 + X^8
nombre d'opérations = 38
III-C. Module complet
On donne pour finir le contenu du module complet :
| Code Python : | Sélectionner tout |
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 | class Polynome: def __init__(self, liste_coefs=[0], nbre_operations=0): # méthode constructeur de la classe # on définit la liste des coefficients du polynôme [a0, a1, ..., an] self.coefs = liste_coefs # on initialise le nombre d'opérations de multiplication entre termes self.nombre_operations = nbre_operations # suppression si nécessaire des zéros en queue de liste de coefficients. Exemple : [2, 3, 1, 0, 0] -> [2, 3, 1] self.reduire() def __str__(self): # permet d'afficher le polynôme sous la forme 1 + 2x + 3x^2 s="" # initialisation de la chaîne de caractères # on vérifie d’abord si le degré du polynôme est nul if (len(self.coefs)-1==0): return str(self.coefs[0]) else: # sinon if self.coefs[0]!=0: s=str(self.coefs[0]) + " + " for i in range(1, len(self.coefs)): # parcours des indices des coefficients du polynôme : [a1, a2, ..., an] if self.coefs[i]!=0: # si le coefficient de degré i n'est pas nul if self.coefs[i]!=1: # si le coefficient de degré i est différent de 1 s+="{}⋅X^{} + ".format(self.coefs[i],i) else: s+="X^{} + ".format(i) # élimination des caractères en trop s = s[:-3].replace("+ -", "- ").replace("X^1 ","X ").replace(" 1⋅X"," X") if s[-2:]=="^1": s = s[:-2] if s[:3]=="1⋅X": s = s[3:] return s # on retourne l'expression du polynôme def degre(self): # retourne le degré du polynôme return (len(self.coefs)-1) def __add__(self, other): # méthode permettant de redéfinir l'opérateur « + » pour 2 polynômes : (1 + 2x + x^2) + (2 + X) = 2 + 3x + x^2 # p1 = self, p2 = other if len(other.coefs) >len(self.coefs): # si degré de p2 > degré de p1 liste_coefs = other.coefs[:]; n = len(self.coefs) # on copie les coefs du polynôme de degré le plus élevé et la longueur de la liste de coefs la plus petite. else: liste_coefs = self.coefs[:]; n = len(other.coefs) # sinon, ... for i in range(n): # parcours des indices de liste_coefs liste_coefs[i] = self.coefs[i] + other.coefs[i] # addition des coefficients de degré i return Polynome(liste_coefs) # renvoie le polynôme résultat de l'addition def __sub__(self, other): # méthode permettant de redéfinir l'opérateur « + » pour 2 polynômes : (1 + 2x + x^2) + (2 + X) = 2 + 3x + x^2 # p1 = self, p2 = other if len(other.coefs) >len(self.coefs): # si degré de p2 > degré de p1 liste_coefs = other.coefs[:]; n = len(self.coefs) # on copie les coefs du polynôme de degré le plus élevé et la longueur de la liste de coefs la plus petite. else: liste_coefs = self.coefs[:]; n = len(other.coefs) # sinon, ... for i in range(n): # parcours des indices de liste_coefs liste_coefs[i] = self.coefs[i] - other.coefs[i] # addition des coefficients de degré i return Polynome(liste_coefs) # renvoie le polynôme résultat de l'addition def reduire(self): # tant que le dernier élément de la liste est nul while round(self.coefs[-1],12) == 0 and len(self.coefs)>1: self.coefs.pop() # supprimer le dernier élément for i in range(len(self.coefs)): self.coefs[i] = round(self.coefs[i],12) def __mul__(self, other): # méthode permettant de redéfinir l'opérateur « * » pour 2 polynômes : (2 + X) * (1 + 2x) = 1 + 3x + 2x^2 # initialisation de la liste des coefficients qu'avec des zéros liste_coefs=[0]*(len(self.coefs)+len(other.coefs)-1) # exemple : [0, 0, 0] for i1 in range(len(self.coefs)): # parcours des indices des coefs du polynôme n°1 for i2 in range(len(other.coefs)): # parcours des indices des coefs du polynôme n°2 # multiplication des coefficients d'indices i1 et i2 liste_coefs[i1+i2] = liste_coefs[i1+i2] + self.coefs[i1]*other.coefs[i2] poly = Polynome(liste_coefs) # création de l'objet Polynome basé sur la liste # ajout du nombre d'opérations de multiplication entre termes effectuées lors de la multiplication entre polynômes poly.nombre_operations = self.nombre_operations + other.nombre_operations + len(self.coefs)*len(other.coefs) # return poly # renvoie le polynôme résultat de la multiplication def __pow__(self, n): # méthode permettant de redéfinir l'opérateur de puissance : self ** n if n==0: # si n=0 return Polynome([1]) # on renvoie le polynôme égal à 1 # copie de self initialisé avec le même nombre d'opérations de multiplication entre termes poly = Polynome(self.coefs, self.nombre_operations) # copie de self avec un nombre d'opérations de multiplication à 0 pol = Polynome(self.coefs) for i in range(n-1): # on multiplie n fois poly par pol à l'aide de l'opérateur * poly = poly*pol # équivalent à : poly = poly.__mul__(pol) return poly # renvoie le polynôme résultat de l'opération (self ** n) def __eq__(poly1, other): # méthode permettant de redéfinir l'opérateur « == » pour 2 polynômes return (poly1.coefs==other.coefs) # renvoie True si les 2 listes de coefficients sont égales def eval(self,x): # méthode permettant d'évaluer le polynôme en fonction de x # intialisation des variables valeur_polynome = self.coefs[0] # valeur_polynome = a0 power=1 for coef in self.coefs[1:]: # parcours des coefficients du polynôme à partir de a1 : a1, a2, ..., an power = power*x # calcul de la puissance de x pour ai : power = x^i valeur_polynome += coef*power # valeur_polynome = valeur_polynome + ai*x^i return valeur_polynome # renvoie la valeur du polynôme def eval_horner(self,x): # méthode permettant d'évaluer le polynôme en fonction de x # intialisation de la variable valeur_polynome = self.coefs[-1] # valeur_polynome = an for coef in reversed(self.coefs[:-1]): # parcours des coefficients du polynôme à partir de an-1 : an-1, ..., a1, a0 valeur_polynome = valeur_polynome*x + coef # valeur_polynome = valeur_polynome*x + ai return valeur_polynome # renvoie la valeur du polynôme def puissance(x,n): # exponentiation rapide de x^n # si n=1 alors on renvoie x if n==1: return x else: # sinon # si n est pair if (n % 2 == 0): # appel récursif : puissance(x^2,n/2) return puissance(x**2,n//2) else: # sinon si n est impair # appel récursif : x*puissance(x^2,(n-1)/2) return x*puissance(x**2,n//2) def puissance_polynome(p,n): # exponentiation rapide de p^n # si n=1 alors on renvoie p if n==1: return p else: # sinon # si n est pair if (n % 2 == 0): # appel récursif : puissance_polynome(p^2,n/2) return puissance_polynome(p**2,n//2) else: # sinon si n est impair # appel récursif : p*puissance_polynome(p^2,(n-1)/2) # return p*puissance_polynome(p**2,n//2) return Polynome(p.coefs)*puissance_polynome(p**2,n//2) # pour calculer le nombre d'opérations de mult. print() print("I-A. Exponentiation de nombres réels") print() # nombre entier x = 12 # exposant n = 8 # exponentiation du polynôme xn = x**n print("x^n = {0}^{1}".format(x,n)) print() print("x^n = " + str(xn)) print() print("nombre de multiplications = " + str(n-1)) print();print() print("I-B. Exponentiation rapide de nombres réels") print() # exponentiation du polynôme xn = puissance(x,n) print("x^n = {0}^{1}".format(x,n)) print() print("x^n = " + str(xn)) print() print("nombre de multiplications de l'ordre de log2({0})".format(n)) print();print() print("II-A. Exponentiation de polynômes") print() # création du polynôme p = 2 + X = X + 2 p = Polynome([2,1]) # exposant n = 8 # exponentiation du polynôme pn = p**n print("p^n = ({0})^{1}".format(p,n)) print() print("p^n = " + str(pn)) print() print("nombre d'opérations = " + str(pn.nombre_operations)) print();print() print("II-B. Exponentiation rapide de polynômes") print() # exponentiation rapide de polynômes # pn = ((p**2)**2)**2 pn = puissance_polynome(p,n) print("p^n = ({0})^{1}".format(p,n)) print() print("p^n = " + str(pn)) print() print("nombre d'opérations = " + str(pn.nombre_operations)) |
IV. Complément : Formule du multinôme de Newton
On peut également se baser sur la formule du multinôme pour développer une puissance de polynômes.
Par exemple, on peut réaliser plus rapidement le développement de (X + 2)n en utilisant la formule du binôme de Newton :
Libre à chacun d'écrire la fonction Python basée sur cette formule.
V. Conclusion
Nous avons donc pu montrer les avantages de l'exponentiation rapide dans le calcul de puissances de nombres réels, puis dans le développement de puissances de polynômes, pour enfin le vérifier en Python à l'aide de la classe Polynome.
Sources :
https://fr.wikipedia.org/wiki/Exponentiation_rapide
https://www.bibmath.net/dico/index.p...ionrapide.html