Exponentiation rapide (Square & Multiply)

Présentation du problème.

On se donne un monoïde \((G,*)\), un élément \(x\in G\) et \(k\in\N\). On cherche à calculer la puis­san­ce \(x^k\), c'est-à-dire la quantité \begin{equation} \underbrace{x*x*\cdots*x}_{k\ \text{termes} \ x}. \end{equation} avec par convention \(x^0=e\) où \(e\) est l'élément neutre de la loi \(*\). Si \((G,*)\) est de surcroît un groupe, on généralise aux exposants \(k\in\Z\) avec \(x^k=(x^{-1})^{-k}\) si \(k<0\) et où \(x^{-1}\) désigne le symétrique de \(x\) pour la loi \(*\) considérée, i.e. \(x*x^{-1}=x^{-1}*x=e\). Ce problème ne se limite donc pas aux calculs de puissances de nombres réels.

Pour ne citer que deux exemples, l'exponentiation est utilisée intensivement dans de nombreux pro­to­co­les cryptographiques, comme rsa sur le groupe \((\Z/n\Z,\times)\) et en algèbre linéaire où les puissances d'une matrice in­ter­vien­nent dans des contextes variés. Nous supposerons donc dans toute la suite que \(x\) est un élément d'un monoïde \((G,*)\) noté multiplicativement, que \(k\) est un entier naturel et que \(e\) est l'élément neutre pour la loi \(*\).

L'algorithme naïf qui consiste à multiplier la valeur \(x\) par elle-même \(k - 1\) fois dans le monoïde a un coût proportionnel à l'exposant \(k\) si l'on suppose que la multiplication a un coût fixe. Cette dernière hypothèse est extrêmement restrictive selon le monoïde considéré. En effet, il est clair que si l'on multiplie des entiers relatifs, la taille de la représentation du résultat ne reste pas bornée, nous y reviendrons au chapitre Calcul multiprécision.

Écrivez l'algorithme d'exponentiation naïf sur la machine ram en supposant qu'un registre est capable de stocker des valeurs arbitrairement grandes.

L'algorithme.

L'idée sur laquelle repose l'algorithme de l'exponentiation binaire est très naturelle. En calculant le carré \(x^2\) de \(x\), puis le carré de \(x^2\), on obtient \(x^4\) en deux multiplications au lieu des trois nécessaires avec l'algorithme naïf. On peut continuer le processus et obtenir \(x^8\) en trois multiplications au lieu de sept, etc. À première vue, ce procédé ne paraît intéressant que pour les exposants \(k\) qui sont des puissances de 2, mais les règles de calcul sur les exponentielles permettent de dépasser cette apparente limitation. En effet, on rappelle que si \(a\) et \(b\) désignent deux entiers, alors \[x^{a+b}=x^a * x^b.\] Ainsi, comme tout entier \(k\) peut s'écrire de manière unique comme une somme de puis­san­ces de \(2\) (c'est ce qu'exprime son écriture en base 2), il existe donc une séquence \((k_i)\) de \(n\) bits tels que \[k:=\sum_{i=0}^{n-1} k_i\, 2^i,\] on a donc \begin{equation} \label{eq:result} x^k=x^{\small\displaystyle\left(\sum_{i=0}^{n-1}k_i\, 2^i\right)}=\prod_{i=0}^{n-1} x^{k_i\, 2^i}. \end{equation} La table suivante illustre le mécanisme pour k = 87 et l'écriture de k en base 2 est donnée du chiffre le plus significatif à gauche au chiffre le moins significatif à droite.

\(k=87\) : 1 0 1 0 1 1 1
Exposant \(i\) : 6 5 4 3 2 1 0
Puissance de 2 : 64 32 16 8 4 2 1
(S)quare : \(\checkmark\) \(\checkmark\) \(\checkmark\) \(\checkmark\) \(\checkmark\) \(\checkmark\) \(\checkmark\)
(M)ultiply : \(\checkmark\) \(\checkmark\) \(\checkmark\) \(\checkmark\) \(\checkmark\)

Écriture binaire de l'entier \(k=87\).
et ainsi \begin{equation} x^{87}=x^{\color{yellow}64}*x^{\color{yellow}16}* x^{\color{yellow}4}*x^{\color{yellow}2}*x^{\color{yellow}1}. \end{equation}

On peut reconstituer \(x^k\) avec deux opérateurs uniquement, l'élévation au carré (noté S pour Square) et la multiplication par x (noté M pour Multiply). En effet, en initialisant une variable \(R\) à l'élément neutre \(e\) pour la multiplication, il suffit de lire successivement les bits de l'écriture binaire de l'entier \(k\) (notés k[i] dans l'algorithme) dans l'ordre des puissances décroissantes de \(2\), et pour chaque bit lu :

  1. On élève \(R\) au carré;
  2. Si le bit est égal à 1, on multiplie \(R\) par \(x\).

L'algorithme ci-dessous formalise ce processus en toute généralité. Le symbole \(*\) désigne la loi de composition du monoïde considéré et \(e\) l'élément neutre pour cette loi : 
SquareMultiply(x,k):valeur;
données
   x: valeur  // dans le monoide considéré
   k: entier de n bits  // exposant
variables
   R: valeur 
   i: entier 
debut
   R ← e
   i ← 0
   TQ (i < n) FAIRE  
      R ← R * R
      SI (k[n - (i + 1)] > 0) ALORS
         R ← R * x
      FSI
      i ← i + 1
   FTQ
   RENVOYER R
fin
Pour calculer \(x^{87}\), on applique donc la séquence d'instructions SM, S, SM, S, SM, SM, SM. Le tableau ci-dessous fournit le contenu de la variable \(R\) à après l'exécution des instructions de la boucle pour la valeur \(k=87\) :

 \(i\)   \(k[n-(i+1)]\)  (S) (M) \(R\) (après opération)
\(e\)
0 1 \(\checkmark\) \(\checkmark\) \(x\)
1 0 \(\checkmark\) \(x^2\)
2 1 \(\checkmark\) \(\checkmark\) \(x^4\times x\)
3 0 \(\checkmark\) \(x^8\times x^2\)
4 1 \(\checkmark\) \(\checkmark\) \(x^{16}\times x^4\times x\)
5 1 \(\checkmark\) \(\checkmark\) \(x^{32}\times x^8\times x^2\times x\)
6 1 \(\checkmark\) \(\checkmark\) \(x^{64}\times x^{16}\times x^4\times x^2\times x\)
Évolution des variables dans l'algorithme d'ex­po­nen­tia­tion bi­nai­re.

Il ne faut pas décomposer \(k\) en binaire, les entiers sont déjà codés en binaire en machine et il existe toujours une manière efficace d'accéder aux chiffres binaires d'un entier sans avoir à faire de décomposition binaire inutile et coûteuse. Le chapitre sur le modèle de donnée ensemble reviendra sur cette question.

La preuve d'arrêt est immédiate, une seule instruction modifie la variable \(i\) dans la boucle, elle est incrémentée, ses valeurs constituent donc une suite strictement croissante de premier terme \(0\), et comme \(N\) est archimédien, la condition d'arrêt \(\color{yellow}i\geq n\) sera donc sa­tis­fai­te après \(n+1\) itérations.

Considérons le prédicat \(P(i)\) suivant : \begin{equation} \label{eq:Rp} \text{après \(i\) itérations de la boucle}\ \ R=\prod_{j=0}^{i-1}x^{k_{n-(j+1)}2^{i-(j+1)}}. \end{equation} Avant d'entrer la boucle, \(P(0)\) est vrai puisque par convention le produit vide est égal à l'élément neutre \(e\) pour l'opé­ra­teur \(*\) considéré. Supposons que \(P(i)\) soit vrai à l'entrée de la boucle. Après l'élévation au carré à la ligne #12, on a \[ R=\prod_{j=0}^{i-1}x^{k_{n-(j+1)}2^{i-(j+1)+1}}. \] En notant que \[x^{k_{n-(i+1)}}= \begin{cases}x&\text{si}\ k_{n-(i+1)}=1\\ 1&\text{si}\ k_{n-(i+1)}=0 \end{cases} \] Le bloc conditionnel #13 à #15 évite simplement le calcul \(x^{k_{n-(i+1)}}\) et ne réalise une multiplication par \(x\) que si c'est nécessaire. Après son exécution, on a donc \[ R=x^{k_{n-(i+1)}}\prod_{j=0}^{i-1}x^{k_{n-(j+1)}2^{i-(j+1)+1}}. \] En réécrivant l'expression (noter que \(2^{i-(i+1)+1}=2^0=1\)), on obtient \[ R=\prod_{j=0}^{i}x^{k_{n-(j+1)}2^{i-(j+1)+1}} \] et \(P(i+1)\) est vrai. Comme \(i=n\) en sortant de la boucle, on a \(P(n)\), soit \begin{align*} R &=\prod_{j=0}^{n-1}x^{k_{n-(j+1)}2^{n-(j+1)}}\\ &=\prod_{i=0}^{n-1} x^{k_i\times 2^i}\quad(\text{en réindexant}\ i:=n-(j+1)) \end{align*} ce qui est bien le résultat attendu \((\ref{eq:result})\).

Quel est le nombre minimal de multiplications permettant de calculer \(x^{31}\) ? Comparer ce nombre au nombre de multiplications effectuées par l'algorithme SquareMultiply.
Écrivez l'algorithme d'exponentiation binaire sur la machine ram en supposant qu'un registre est capable de stocker des valeurs arbitrairement grandes.

Complexité.

Plutôt que la complexité au sens strict, nous allons évaluer le coût \(C(n)\) de l'algorithme exprimé en nombre de produits en fonction du nombre de chiffres de l'exposant \(k\). Ce choix est légimité par le fait que l'al­go­ri­thme s'applique à n'importe quel monoïde, un calcul de complexité strict imposerait d'intégrer le coût de chaque produit, or l'algorithme auxiliaire qui effectue ce produit dépend de la nature du monoïde considéré, un produit d'entiers n'est pas équi­va­lent à un produit de matrices par exemple.

Nous verrons à ce sujet dans le cadre du calcul multiprécision que les opérations arithmétiques sur les entiers n'ont pas un coût constant et que ce coût dépend de la taille des nom­bres que l'on manipule. L'algorithme de la multiplication étudié dans les classes du primaire a un coût quadratique en la taille des opérandes par exemple.

La preuve du lemme suivant est à faire en exercice (voir plus bas), ne la dépliez qu'après avoir tenté votre chance…

Soit \(N\) et \(b\) deux entiers naturels avec \(b\geq 2\). Le nombre de chiffres de l'écriture de \(N\) en base \(b\) est égal à \(\lfloor\log_{b}N\rfloor+1.\)
On note \(n\) le nombre de chiffres de l'écriture de l'entier \(N\) en base \(b\) et \(N_i,\ i\in[0,n-1]\) ses chiffres de manière à ce que \begin{equation*} N=\sum_{i=0}^{n-1}N_ib^i. \end{equation*} Puisque \(N\) a \(n\) chiffres en base \(b\), son chiffre le plus significatif \(N_{n-1}\) n'est pas nul, \(N\) est donc supérieur à \(b^{n-1}\). D'autre part, le plus grand nombre à \(n\) chiffres en base \(b\) est le nombre dont les \(n\) chiffres sont tous égaux à \(b-1\) qui est égal à \(b^n-1\). On en déduit l'encadrement : \begin{equation} b^{n-1}\leq N < b^n. \end{equation} La fonction logarithme en base \(b\) est une fonction croissante, on en déduit que \begin{equation} n-1\leq \log_b N < n. \end{equation} Si l'on se rappelle que la partie entière \(\lfloor x\rfloor\) d'un nombre réel \(x\) est l'unique entier \(e\) tel que \(e\leq x < e+1\), alors \(n-1=\lfloor\log_b N\rfloor\) et le résultat en découle.
Soit \(q\) un entier naturel non-nul. On appelle poids \(q\)-aire d'un entier naturel \(N\) le nombre de ses chiffres non-nuls dans son écriture en base \(q\) et que l'on note \(\text{wt}_q(N)\).

Notons \(n\) le nombre de bits de l'exposant \(k\), c'est-à-dire \(n:=\lfloor\log_{2}k\rfloor+1\) d'après le lemme ci-dessus. Nous savons que le nombre de multiplications utilisées pour les élévations au carré (S) est exactement \(n\) et que le nombre de mul­ti­pli­ca­tions pour l'opération (M) est égal au poids binaire \(\text{wt}_2(k)\) de l'entier \(k\). Comme le nombre de bits \(n\) est fixé, le bit de poids fort de \(k\) est nécessairement égal à \(1\) et donc \(1\leq\text{wt}_2(k)\leq n\). Au total, le nombre de multiplications pour traiter l'exposant \(k\) de \(n\) bits est égal à \(n+\text{wt}_2(k)\) et on a l'en­ca­dre­ment suivant : \begin{equation} n+1\leq C(n)\leq 2n. \end{equation} si le nombre de multiplications \(C(n)\) est exprimé en fonction du nombre de bits de \(k\). On en déduit directement la proposition

Le nombre de produits effectués par l'algorithme SquareMultiply dans le meilleur des cas et dans le pire est donné respectivement par les fonctions suivantes : \begin{align} \check{C}(n)&=n+1.\\ \hat{C}(n)&=2n. \end{align}

Il reste à évaluer le cas moyen. On sait que quelle que soit l'instance \(k\) de taille \(n\) que l'on considère pour l'exposant, le nombre de multiplications pour calculer les carrés (S) est toujours égal à \(n\). D'autre part le bit de poids fort de \(k\) étant nécessairement égal à 1, il y a toujours une multiplication (M) à la première étape. Ce sont les \(n-1\) bits de poids faible qui vont fixer le nombre de mul­ti­pli­ca­tions qu'il reste à effectuer. Il y a donc \(2^{n-1}\) instances différentes à étudier, correspondant aux dif­férentes valeurs possibles des \(n-1\) bits de poids faible.

Le nombre d'opérations (M) étant égal au poids binaire de \(k\), on peut calculer la moyenne en par­ti­tion­nant l'ensemble \(K\) des séquences binaires de longueur \(n-1\) en \(n\) classes selon les différents poids possibles :

\[K_p:=\left\{u\in\{0,1\}^{n-1}\mid \text{wt}_2(u) = p\right\},\quad p\in\ab{0}{n-1}.\]

On a donc à ce stade

\[\overline{C}(n)=n+1+\frac{1}{2^{n-1}} \sum_{p=0}^{n-1}p|K_p|. \]

On veut dans un premier temps calculer les cardinaux \(|K_p|\). On rappelle que si \(p\) et \(n\) désignent deux entiers naturels avec \(0\leq p\leq n\), le nombre de parties à \(p\) éléments dans un ensemble à \(n\) élements, noté \(\binom{n}{p}\) est égal à

\[\binom{n}{p}=\frac{n!}{(n-p)!p!}.\]

En montre aisément que

Il existe exactement \(\binom{n}{p}\) séquences binaires de longueur \(n\) et de poids binaire \(p\).

On en déduit que \(|K_p|=\binom{n}{p}\), le coût moyen devient

\[\overline{C}(n)=n+1+\frac{1}{2^{n-1}} \color{yellow}{\sum_{p=0}^{n-1}p\binom{n-1}{p}}. \]

Il nous reste à calculer la somme ci-dessus.

Soient \(p\) et \(n\) deux entiers na effectués par l'algorithme Square­Multiply est donné par la fonction suivante : \begin{equation}\label{eq:moyenne} \overline{C}(n)=\frac {3n+1}{2}\quad\text{où}\quad n=\lfloor\log_2 k\rfloor+1 \end{equation}
Démontrez le lemme 2. Rappel : on admettra que pour tout nombre réel \(x\geq 0\), il existe un unique entier naturel \(N\) tel que \(N\leq x < N+1\), cet entier est appelé la partie entière de \(x\), notée \(\lfloor x\rfloor\).
Démontrez le lemme 5.
Soient \(p\) et \(n\) deux entiers naturels tels que \(0\leq p\leq n\). Écrivez un algorithme qui calcule le nombre de combinaisons de \(p\) parmi \(n\) en exploitant les trois identités suivantes : \begin{equation*} \binom{n}{p}=\binom{n}{n-p},\qquad\binom{n+1}{p+1}=\frac{n+1}{p+1}\binom{n}{p},\qquad\binom{n}{0}=1. \end{equation*} Comparez sa complexité avec celle du calcul direct qui s'appuie sur la fonction factorielle dont on intègrera le coût.
Soit \(n\) un entier naturel tel que \(n\geq 1\) et \(B\) une matrice \(2^n\times n\) dont les lignes sont cons­ti­tuées par les \(2^n\) séquences binaires distinctes de longueur \(n\). Démontrez que chaque colonne de cette matrice contient exactement \(2^{n-1}\) valeurs \(0\) et autant de valeurs \(1\). En déduire le résultat \((\ref{eq:moyenne})\).
En considérant que l'algorithme Square & Multiply s'applique à l'exponentiation d'un nombre entier et que la multiplication est réalisée à l'aide de l'algorithme naïf étudié au primaire, calculez sa complexité dans le pire des cas en notant \(n\) la taille de l'exposant et \(m\) le nombre de chiffres de l'entier \(X\).

Travaux pratiques

Le langage utilisé dans tous les sujets qui suivent est le \(C\).
Définissez un type ullong pour des entiers unsigned long long puis un type tmat pour coder une matrice \(2\times 2\) de nombres entiers (de type ullong) (attention, une fonction en C ne peut pas renvoyer un tableau. Si vous voulez utiliser un tableau, rangez le dans une structure ou faites un malloc). Écrivez une fonction tmat produit(tmat A, tmat B) qui renvoie le produit matriciel \(A\times B\).
Écrivez une fonction tmat puissance(tmat M, uint n) qui calcule \(M^n\) avec l'al­go­ri­thme naïf qui effectue \(n-1\) produits.
Écrivez une fonction tmat SM(tmat M, uint n) qui calcule \(M^n\) avec l'al­go­ri­thme SquareMultiply et renvoie le résultat.
La suite de Fibonacci est une suite récurrente de premiers termes \(F_0:=0\) et \(F_1:=1\) définie par la relation \begin{equation} \forall n\in{\mathbb N}\quad F_{n+2}:=F_{n+1}+F_n. \end{equation}
  1. Écrivez une fonction ullong Fib(uint n) qui pour un entier \(n\) en entrée, calcule le terme \(F_n\) de la suite de Fibonacci.
  2. On définit la matrice \(2\times 2\) : \begin{equation} \Phi:=\begin{pmatrix} 1&1\\1&0 \end{pmatrix}, \end{equation} et on admet que \begin{equation} \forall n\in{\mathbb N}^*\quad \Phi^n=\begin{pmatrix} F_{n+1}&F_n\\F_n&F_{n-1}. \end{pmatrix} \end{equation}
  3. Écrivez une fonction ullong FibSM(uint n) qui exploite cette relation et utilise l'algorithme SquareMultiply pour calculer \(F_n\).
  4. Quelle est la plus grande valeur de \(n\) pour laquelle on peut calculer \(F_n\) sans dépassement de capacité d'un registre ?
On veut comparer les temps d'exécution des deux fonctions Fib et FibSM. Pour mesurer la durée d'exécution d'un programme exécutable a.out, utilisez

/usr/bin/time --format="%U" a.out >/dev/null

Mesurez le temps d'exécution de vos deux programmes pour des valeurs croissantes de \(n\) sans dépasser la capacité d'un ullong. Tracez les graphes des temps d'exécution en fonction de \(n\) de ces deux algorithmes avec la commande Gnuplot.