IdentifiantMot de passe
Loading...
Mot de passe oublié ?Je m'inscris ! (gratuit)
Voir le flux RSS

User

[Actualité] Exponentiation rapide de nombres réels et de polynômes en Python

Noter ce billet
par , 24/01/2024 à 09h47 (7588 Affichages)
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

Nom : exponentiation_rapide1.png
Affichages : 6855
Taille : 4,9 Ko


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 :

Nom : algo_exponentiation1.png
Affichages : 2066
Taille : 7,1 Ko

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 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 - 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 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 :

Nom : algo_exponentiation2.png
Affichages : 2077
Taille : 7,8 Ko

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)n en utilisant la formule du binôme de Newton :

Nom : formule_binome.png
Affichages : 1979
Taille : 3,9 Ko

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

Envoyer le billet « Exponentiation rapide de nombres réels et de polynômes en Python » dans le blog Viadeo Envoyer le billet « Exponentiation rapide de nombres réels et de polynômes en Python » dans le blog Twitter Envoyer le billet « Exponentiation rapide de nombres réels et de polynômes en Python » dans le blog Google Envoyer le billet « Exponentiation rapide de nombres réels et de polynômes en Python » dans le blog Facebook Envoyer le billet « Exponentiation rapide de nombres réels et de polynômes en Python » dans le blog Digg Envoyer le billet « Exponentiation rapide de nombres réels et de polynômes en Python » dans le blog Delicious Envoyer le billet « Exponentiation rapide de nombres réels et de polynômes en Python » dans le blog MySpace Envoyer le billet « Exponentiation rapide de nombres réels et de polynômes en Python » dans le blog Yahoo

Mis à jour 07/02/2024 à 06h38 par User

Catégories
Programmation , Python , Algorithmique

Commentaires