Mathématiques et Python - Exponentiation rapide de nombres réels et de polynômes en Python

Un billet blog de Denis Hulo

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 :


Testons maintenant notre fonction avec ces quelques lignes :


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 :


Le code de la méthode permettant de multiplier deux polynômes devient alors :


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 :


Testons maintenant l'opérateur de puissance « ** » pour les polynômes avec la méthode classique :



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 :


Testons notre fonction avec ces quelques lignes :



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]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):[/1]...