User

I. Introduction

Dans 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 aⁿ, on peut naïvement effectuer n−1 multiplications en réalisant le produit :

aⁿ = a × a × ... × a

Cependant, l'exemple suivant montre que l'on peut calculer a³² en effectuant seulement 5 multiplications au lieu de 31 :

a³² = a²⁵ = ((((a²)²)²)²)²

Si l'exposant est différent d'une puissance de 2, avec par exemple a⁴¹, on peut alors écrire :

a⁴¹ = a × a⁸ × a³² = a × a²³ × a²⁵

En remarquant que a⁸ est également calculé dans a³², 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 aⁿ, 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 n² produits. La multiplication (X + 2) × (X + 2) nécessite ainsi 4 multiplications entre termes.

Si on prend maintenant :

(X + 2)⁸ = (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)⁸ = ((X + 2) × (X + 2)) × (X + 2) × (X + 2) × (X + 2) × (X + 2) × (X + 2) × (X + 2)

La 1^re multiplication nécessite donc 4 opérations de multiplication entre termes.

En poursuivant le développement :

(X + 2)⁸ = ((X² + 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)⁸ = (X + 2)²³ = (((X + 2)²)²)²

Si on effectue le développement en respectant la règle de priorité des parenthèses, on obtient successivement :

(X + 2)⁸ = (((X + 2)²)²)²

2×2=4 opérations de multiplication entre termes pour la 1^re produit.

Puis :

(X + 2)⁸ = ((X² + 4X + 4))²)²

3×3=9 opérations de multiplication pour le 2^e produit.

Et enfin :

(X + 2)⁸ = (X⁴ + 8X³ + 24X² + 32X + 16)²

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 x^k de x, pour tout k, tel que 1 ≤ k < n.

• Si n est pair alors xⁿ = (x²)^n/2. Il suffit alors de calculer y^n/2 pour y = x².
• Si n est impair et n > 1, alors xⁿ = x(x²)^{(n – 1)/2}. Il suffit de calculer y^{(n – 1)/2} pour y = x² et de multiplier le résultat par x.

Cette remarque nous amène à l'algorithme récursif suivant qui calcule xⁿ pour un entier strictement positif n :



On obtient ainsi la fonction d'exponentiation rapide de nombres réels :

Code Python :

Sélectionner tout - Visualiser dans une fenêtre à part

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 - Visualiser dans une fenêtre à part

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 - Visualiser dans une fenêtre à part

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 - Visualiser dans une fenêtre à part

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 X³ + 2 = 2 + 0.X + 0.X² + X³ 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 - Visualiser dans une fenêtre à part

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 - Visualiser dans une fenêtre à part

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 p^k de p, pour tout k, tel que 1 ≤ k < n.

• Si n est pair alors pⁿ = (p²)^n/2. Il suffit alors de développer q^n/2 pour q = p².
• Si n est impair et n > 1, alors pⁿ = p(p²)^{(n – 1)/2}. Il suffit de développer q^{(n – 1)/2} pour q = p² et de multiplier le résultat par p.

Cette remarque nous amène à l'algorithme récursif suivant qui développe pⁿ 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 - Visualiser dans une fenêtre à part

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 - Visualiser dans une fenêtre à part

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 - Visualiser dans une fenêtre à part

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)ⁿ 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