Calcul de grosses puissances de gros nombres dan un anneau modulaire
Nous avons vu que les puissances dans un anneau modulaire sont importants quand on veut trouver de gros nombres premiers. Aussi, les puissances dans un anneau modulaire ont des applications cryptographiques importantes. Bien sûr, une façon triviale de calculer une puissance aN consiste à effectuer (N-1) multiplications de a avec lui-même. Cela n'est cependant pas faisable quand nous avons à faire à des exposants qui sont de très grands nombres. Le nombre d'itérations nécessaires feraient que le calcul prend plus de temps que l'age de l'univers.
Cependant, il y a une façon bien plus rapide de calculer une telle puissance. Considérons: a2 = a * a ; a4 = a2 * a2 ; a8 = a4 * a4 ; a16 = a8 * a8 ; a32 = a16 * a16 etc. Comme nous le constatons, pour calculer a32 il suffit de 5 multiplications, pas de 31. Si nous avons besoin de calculer, disons, a38 on peut le faire ainsi: a32 * a4 * a2 ;2 multiplications supplémentaires. Il s'avère que cette technique est universelle. Pour calculer la puissance N d'un nombre a, il suffit de calculer les puissances de a qui ont des exposants qui sont eux-mêmes les puissances successives de 2 (donc, 1, 2, 4, 8, 16, 32, 64, 128 ...) et de multiplier ensemble certains de ces facteurs. Les facteurs qu'il faut prendre sont donnés par le développement binaire du nombre N: chaque bit "1" indique un facteur correspondant à prendre.
Quel que soit la base dans laquelle est exprimée N, on peut trouver sa représentation binaire en prenant les divisions successives Euclidiennes de N par 2. Les valeurs du reste successives sont les bits, commençant par le lsb. Par exemple, si nous considérons (dans le système décimal) le nombre 23, alors:
- 23 divisé par 2 donne 11, reste 1. 1 est donc le bit lsb de 23.
- 11 divisé par 2 donne 5, reste 1. 1 est donc le deuxième bit lsb de 23.
- 5 divisé par 2 donne 2, reste 1. 1 est donc le troisième bit lsb de 23
- 2 divisé par 2 donne 1, reste 0. 0 est le quatrième bit lsb de 23.
- 1 divisé par 2 donne 0, reste 1. 1 est le bit msb de 23.
La représentation binaire de 23 est donc: 10111
Pour faciliter ces calculs, nous introduisons une fonction auxiliaire : "division Euclidienne par 2"
Nous voulons rappeler: Exemples de code
void divide2(number* x, number* q, int* r)
/* Euclidean division by two */
{
int i;
int rest;
rest = 0;
int divi;
for (i = nofdigits - 1 ; i >= 0 ; i--)
{
divi = rest * basis + x->digit[i];
q->digit[i] = divi / 2;
rest = divi % 2;
}
*r = rest;
} // end divide2
Avec l'aide de cette fonction, nous pouvons introduire finalement la fonction qui calcule la puissance d'un nombre élevée à un autre nombre dans un anneau modulaire:
void power(number* x, number* y, number* xy)
/* number x to the power y */
{
int i;
int j;
number p;
number s;
number zero;
number half;
number full;
int uneven;
int bitsequence[basisbits * nofdigits];
int bitcounter;
for (i = 0 ; i < nofdigits ; i++)
{
xy->digit[i] = 0;
zero.digit[i] = 0;
full.digit[i] = y->digit[i];
}
xy->digit[0] = 1;
bitcounter = 0;
/* we find the binary representation of the power */
do
{
divide2(&full,&half,&uneven);
bitsequence[bitcounter] = uneven;
for (i = 0 ; i < nofdigits ; i++)
{
full.digit[i] = half.digit[i];
}
bitcounter++;
}
while (isbigger(&half,&zero) > 0);
/* the square-and-multiply algorithm */
for (i = bitcounter - 1 ; i >= 0 ; i--)
{
timesnum(xy,xy,&s); // we square
timesnum(&s,x,&p); // we multiply
if (bitsequence[i] == 1)
{
for (j = 0 ; j < nofdigits ; j++)
{
xy->digit[j] = p.digit[j];
}
}
else
{
for (j = 0 ; j < nofdigits ; j++)
{
xy->digit[j] = s.digit[j];
}
}
}
} // end power
Dans cette fonction, nous déterminons d'abord la représentation binaire de l'exposant. En suite, nous appliquons le système de prise de carré et multiplication de la base, en utilisant le développement binaire pour décider si nous "multiplions et prenons le carré" ou si nous prenons seulement le carré.
Le test de primalité de Fermat
Maintenant que nous pouvons calculer des puissances, l'implémentation du test de Fermat n'est pas difficile:
int isprimefermat(number* x, int k)
/* applies the Fermat primality test to the number x with k trials randomly selected in the range 2..x-1 */
/* the function returns 1 when the number x is most probably a prime, and 0 when it is for sure not a prime */
{
number oldmodulo;
number zero;
number one;
number testaraw;
number q;
number testa;
number xminus1;
number testresult;
int kcnt;
int isstillprime;
int i;
srand((unsigned int)time(NULL));
/* we put temporarily the modulus equal to x */
for (i = 0 ; i < nofdigits ; i++)
{
oldmodulo.digit[i] = modulo.digit[i];
modulo.digit[i] = x->digit[i];
zero.digit[i] = 0;
one.digit[i] = 0;
}
one.digit[0] = 1;
difnum(x,&one,&xminus1);
isstillprime = 0;
for (kcnt = 0; kcnt < k ; kcnt++)
{
/* we generate a new random test number that is not zero */
do
{
for (i = 0 ; i < nofdigits ; i++)
{
testaraw.digit[i] = rand() % basis ;
}
dividenum(&testaraw,&modulo,&q,&testa);
}
while (isbigger(&testa,&zero) == 0);
/* we calculate a^(x-1)*/
power(&testa,&xminus1,&testresult);
/* we do the Fermat test */
if (isbigger(&testresult,&one) != 0)
{
isstillprime = 1;
}
} // end of loop over number of random tests
/* we put the original modulo parameter back */
for (i = 0 ; i < nofdigits ; i++)
{
modulo.digit[i] = oldmodulo.digit[i];
}
return 1 - isstillprime;
} // end isprimefermat
Au départ, nous remplaçons l'actuel module actif, par le candidat nombre premier, car tous les calculs du test de Fermat se font dans l'anneau modulaire avec module égal au candidat. A la fin de la routine, nous mettons le module d'origine à sa place.
Un nombre aléatoire est généré. La qualité de ce nombre aléatoire n'est pas parfaite. D'abord, nous utilisons la fonction rand() de la librairie standard C, qui est connue pour ne pas avoir des performances fantastiques comme générateur de nombres aléatoires. Ensuite, nous utilisons simplement l'opération module en C pour restreindre le nombre aléatoire à la base. Finalement, nous "replions" ce nombre dans l'anneau modulaire. Ces méthodes simples peuvent légèrement mettre en cause l'uniformité de la distribution du nombre aléatoire tiré. Le nombre ne doit pas être 0, bien sûr. Si c'est le cas, nous en tirons un autre.
Le test lui-même consiste en fait simplement d'un calcul de puissance, et un test d'égalité à 1.
Le test de primalité Miller-Rabin
Le test de Miller-Rabin est un peu plus compliqué mais n'a besoin de rien d'autre que les opérations arithmétiques que nous avons déjà implémentées:
int isprimerabin(number* x, int k)
/* applies the Miller-Rabin primality test to the number x with k trials randomly selected in the range 2..x-1 */
/* the function returns 1 when the number x is most probably a prime, and 0 when it is for sure not a prime */
{
number oldmodulo;
number zero;
number one;
number testaraw;
number q;
number testa;
number testa2;
number xminus1;
number rabin;
number testresult;
int kcnt;
int isstillprime;
int i;
int u;
int uu;
int r;
srand((unsigned int)time(NULL));
/* we put temporarily the modulus equal to x */
for (i = 0 ; i < nofdigits ; i++)
{
oldmodulo.digit[i] = modulo.digit[i];
modulo.digit[i] = x->digit[i];
zero.digit[i] = 0;
one.digit[i] = 0;
}
one.digit[0] = 1;
/* we obtain the composition of x - 1 = 2^u * r */
difnum(x,&one,&rabin);
difnum(x,&one,&xminus1);
sumnum(&rabin,&zero,&q);
u = 0;
do
{
sumnum(&q,&zero,&rabin);
divide2(&rabin,&q,&r);
u++;
}
while (r == 0);
u--;
/* at this point, the rabin exponent r is in rabin, and u has the power of 2 */
isstillprime = 0;
for (kcnt = 0; kcnt < k ; kcnt++)
{
/* we generate a new random test number that is not zero */
do
{
for (i = 0 ; i < nofdigits ; i++)
{
testaraw.digit[i] = rand() % basis ;
}
dividenum(&testaraw,&modulo,&q,&testa);
}
while (isbigger(&testa,&zero) == 0 || isbigger(&testa,&one) == 0 || isbigger(&testa,&xminus1) == 0);
/* we start the Miller-Rabin test */
power(&testa,&rabin,&testresult);
/* the first part of the test verifies whether this test results in a 1 or a x-1 */
if (isbigger(&testresult,&one) != 0 && isbigger(&testresult,&xminus1) != 0)
{
uu = 0;
do
{
uu++;
// we take the square of the a^m and put it back in testa
timesnum(&testresult,&testresult,&testa2);
sumnum(&testa2,&zero,&testresult);
if (isbigger(&testresult,&one) == 0)
{
isstillprime = 1;
}
} // end of loop over squaring
while (uu < u && isbigger(&testresult,&xminus1) != 0);
// final test on last entity
if (isbigger(&testresult,&xminus1) != 0)
{
isstillprime = 1;
}
} // end need for a loop test
} // next random number
/* we put the original modulo parameter back */
for (i = 0 ; i < nofdigits ; i++)
{
modulo.digit[i] = oldmodulo.digit[i];
}
return 1 - isstillprime;
} // end isprimefermat
La première partie du code trouve la décomposition de p - 1 comme une puissance de 2 fois un nombre impair. Ensuite, nous appliquons l'algorithme de Miller-Rabin, en utilisant une seule puissance, et un nombre de prises au carré. D'abord, il faut aussi générer un nombre aléatoire de test, comme dans le cas du test de Fermat. Cette fois, il faut aussi éviter que ce nombre soit 1 ou p - 1.
