Algorithme rho de Pollard (logarithme discret)

Considérons un groupe cyclique ${\displaystyle G}$ d'ordre ${\displaystyle n}$ engendré par ${\displaystyle \alpha }$. Le but est de calculer ${\displaystyle \gamma }$ tel que ${\displaystyle \alpha ^{\gamma }=\beta }$, où ${\displaystyle \beta }$ appartient à ${\displaystyle G}$. L'algorithme cherche des entiers ${\displaystyle a}$, ${\displaystyle b}$, ${\displaystyle A}$, et ${\displaystyle B}$ tels que ${\displaystyle \alpha ^{a}\beta ^{b}=\alpha ^{A}\beta ^{B}}$. Une fois de tel entier trouvé, l'inconnue ${\displaystyle \gamma }$ est l'une des solutions de l'équation ${\displaystyle (B-b)\gamma =(a-A){\pmod {n}}}$, autrement dit ${\displaystyle \gamma =(B-b)^{-1}(a-A){\pmod {n}}}$, où ${\displaystyle (B-b)^{-1}}$est l'inverse modulaire de ${\displaystyle (B-b)}$.

Pour trouver les entiers ${\displaystyle a}$, ${\displaystyle b}$, ${\displaystyle A}$, et ${\displaystyle B}$ l'algorithme utilise l'algorithme du lièvre et de la tortue pour trouver un cycle dans les séquences ${\displaystyle x_{i}=\alpha ^{a_{i}}\beta ^{b_{i}}}$. Autrement dit les entiers ${\displaystyle a}$, ${\displaystyle b}$ sont des valeurs prises par ${\displaystyle a_{i}}$, ${\displaystyle b_{i}}$ et les entiers ${\displaystyle A}$ et ${\displaystyle B}$ sont les valeurs prises par ${\displaystyle a_{2i}}$, ${\displaystyle b_{2i}}$.

On définit la suite ${\displaystyle (x_{i})_{i\in \mathbb {N} }}$ par ${\displaystyle x_{i+1}=f(x_{i})}$ où ${\displaystyle f}$ est une fonction supposée pseudo-aléatoire : ainsi on a des chances d'entrer dans une boucle de longueur approximative ${\displaystyle {\sqrt {\frac {\pi n}{8}}}}$ après ${\displaystyle {\sqrt {\frac {\pi n}{8}}}}$ étapes[1]. Un moyen de définir une telle fonction est d'utiliser les règles suivantes : diviser ${\displaystyle G}$ en trois sous-ensembles disjoints de tailles approximativement égales : ${\displaystyle S_{0}}$, ${\displaystyle S_{1}}$, et ${\displaystyle S_{2}}$. Si ${\displaystyle x_{i}}$ appartient à ${\displaystyle S_{0}}$ alors multiplier par deux ${\displaystyle a}$ et ${\displaystyle b}$ ; si ${\displaystyle x_{i}\in S_{1}}$ alors incrémenter ${\displaystyle a}$, si ${\displaystyle x_{i}\in S_{2}}$ alors incrémenter ${\displaystyle b}$.

Soit G un groupe cyclique d'ordre n. Soient ${\displaystyle \alpha ,\beta \in G}$. Soit une partition ${\displaystyle G=S_{0}\cup S_{1}\cup S_{2}}$.

Soit une fonction ${\displaystyle f:G\to G}$ définie par

${\displaystyle f(x)={\begin{cases}\beta x&x\in S_{0}\\x^{2}&x\in S_{1}\\\alpha x&x\in S_{2}\end{cases}}}$

et soient enfin les deux fonctions ${\displaystyle g:G\times \mathbb {Z} \to \mathbb {Z} }$ et ${\displaystyle h:G\times \mathbb {Z} \to \mathbb {Z} }$ définies par

${\displaystyle {\begin{aligned}g(x,a)&={\begin{cases}a&x\in S_{0}\\2a\ ({\bmod {\ }}n)&x\in S_{1}\\a+1\ ({\bmod {\ }}n)&x\in S_{2}\end{cases}}\\h(x,b)&={\begin{cases}b+1\ ({\bmod {\ }}n)&x\in S_{0}\\2b\ ({\bmod {\ }}n)&x\in S_{1}\\b&x\in S_{2}\end{cases}}\end{aligned}}}$

La fonction g (respectivement h) permet de calculer les valeurs a_i (respectivement b_i). L'algorithme est le suivant :

 Entrées : ${\displaystyle \alpha }$: un générateur de G, ${\displaystyle \beta }$: un élément de G
 Sortie : un entier x tel que ${\displaystyle \alpha }$x = ${\displaystyle \beta }$, ou "échec"
 a0 ← 0, b0 ← 0, x0 ← 1
 i ← 1
 boucle
     xi ← f(xi-1), 
     ai ← g(xi-1, ai-1), 
     bi ← h(xi-1, bi-1)
     x2i ← f(f(x2i-2)), 
     a2i ← g(f(x2i-2), g(x2i-2, a2i-2)), 
     b2i ← h(f(x2i-2), h(x2i-2, b2i-2))
     si xi = x2i alors
         r ← bi - b2i
         si r = 0 retourner "échec"
         retourner r−1(a2i - ai) mod n
     sinon # xi ≠ x2i
         i ← i+1

Considérons par exemple le groupe engendré par 2 modulo ${\displaystyle N=1019}$ (l'ordre du groupe est ${\displaystyle n=1018}$). L'algorithme est implémenté par le programme suivant écrit en C++ :

 #include <stdio.h>
 
 const int n = 1018, N = n + 1;  /* N = 1019 -- prime     */
 const int alpha = 2;            /* generator             */
 const int beta = 5;             /* 2^{10} = 1024 = 5 (N) */
 
 void new_xab( int& x, int& a, int& b ) {
   switch( x%3 ) {
   case 0: x = x*x     % N;  a =  a*2  % n;  b =  b*2  % n;  break;
   case 1: x = x*alpha % N;  a = (a+1) % n;                  break;
   case 2: x = x*beta  % N;                  b = (b+1) % n;  break;
   }
 }
 
 int main(void) {
   int x=1, a=0, b=0;
   int X=x, A=a, B=b;
   for(int i = 1; i < n; ++i ) {
     new_xab( x, a, b );
     new_xab( X, A, B );
     new_xab( X, A, B );
     printf( "%3d  %4d %3d %3d  %4d %3d %3d\n", i, x, a, b, X, A, B );
     if( x == X ) break;
   }
   return 0;
 }

Les résultats sont les suivants (édités):

 i     x   a   b     X   A   B
------------------------------
 1     2   1   0    10   1   1
 2    10   1   1   100   2   2
 3    20   2   1  1000   3   3
 4   100   2   2   425   8   6
 5   200   3   2   436  16  14
 6  1000   3   3   284  17  15
 7   981   4   3   986  17  17
 8   425   8   6   194  17  19
..............................
48   224 680 376    86 299 412
49   101 680 377   860 300 413
50   505 680 378   101 300 415
51  1010 681 378  1010 301 416

Ce qui donne ${\displaystyle 2^{681}5^{378}=1010=2^{301}5^{416}{\pmod {1019}}}$ et ainsi ${\displaystyle (416-378)\gamma =681-301{\pmod {1018}}}$, pour lequel ${\displaystyle \gamma _{1}=10}$ est une solution comme on l'avait prévu. Comme ${\displaystyle n=1018}$ n'est pas premier, il y a une autre solution ${\displaystyle \gamma _{2}=519}$, pour laquelle ${\displaystyle 2^{519}=1014=-5{\pmod {1019}}}$ convient.

Le temps d'exécution moyen est en ${\displaystyle {\mathcal {O}}({\sqrt {n}})}$. Si cet algorithme est utilisé avec l'algorithme de Pohlig-Hellman, le temps d'exécution des deux algorithmes combinés est en ${\displaystyle {\mathcal {O}}({\sqrt {p}})}$, où ${\displaystyle p}$ est le plus grand facteur premier de ${\displaystyle n}$.