Lipsum.dev

Nombres premiers, Miller-Rabin et OpenSSL

27 juin 2021 — Cryptographie, Python

La commande ci-dessous permet de générer un nombre premier sur 256 bits :

$ openssl prime -generate -bits 256
97883500005213530030432141140344527501513253652173335444305361382001588159937

Cette fonctionnalité d’OpenSSL est utilisée notamment pour générer des clefs privées RSA (openssl genrsa, on pourra consulter cette entrée au sujet de RSA).

Nous allons présenter ici les quelques concepts mis en oeuvre dans l’algorithme de génération de nombres premiers utilisé par OpenSSL.

Test de primalité

Les nombres premiers sont les entiers n>1n > 1 qui n’ont que deux diviseurs : 11 et nn lui-même.

Dans un premier temps, nous nous intéressons au problème de décision suivant :

Étant donné un nombre entier positif nn, nn est-il premier ?

La fonction basique ci-dessous, en Python, permet de déterminer si un nombre donné est premier :

def is_prime(n):
if n <= 1:
return False
for i in range(2, n):
if n % i == 0:
return False
return True
print([i for i in range(100) if is_prime(i)])

On affiche ainsi les nombres premiers inférieurs à 100 : 2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61, 67, 71, 73, 79, 83, 89, 97.

Cependant, cette fonction montre ses limites dès qu’on veut tester des nombres dont les facteurs non triviaux sont « grands » :

is_prime(4208312609**2) # Hangs forever

Test de Fermat

Le Petit théorème de Fermat, énoncé vers 1640, peut se formuler ainsi :

Si pp est un nombre premier et aa un entier quelconque, alors apaa^p - a est un multiple de pp.

Exercice :

  • En remarquant que apa^p correspond au nombre de mots de pp lettres que l’on peut composer avec un alphabet de aa caractères et en comptant séparément les mots constitués d’une même lettre (e.g. BBBBB) et les autres, retrouver le théorème ci-dessus.
  • En supposant de plus 1a<p1 \leq a \lt p, montrer que l’égalité du théorème se simplifie en ap1=1modpa^{p-1} = 1 \mod p.

Ceci nous donne un critère de non-primalité d’un entier, que l’on peut illustrer par le programme suivant :

def fermat_test(n, max_witness=100):
witnesses = range(2, min(n, max_witness + 1))
for a in witnesses:
if pow(a, n - 1, n) != 1:
print(f"{n} is not a prime ({a} is a witness)")
return
print(f"{n} might be a prime")

Par exemple :

>>> fermat_test(4208312609) # 4208312609 is prime
4208312609 might be a prime
>>> fermat_test(4208312609**2)
17709895015068386881 is not a prime (2 is a witness)
>>> fermat_test(5472940991761) # 5472940991761 = 199 * 241 * 863 * 132233
5472940991761 might be a prime
Limitations

Le nombre 54729409917615472940991761 (identifié dans cet article) est un nombre de Carmichael : ce sont les nombres nn qui ne sont pas premiers mais vérifient an1=1modna^{n - 1} = 1 \mod n pour tout entier aa premier avec nn.

Ils sont relativement rares, mais suffisent à mettre en défaut l’algorithme ci-dessus, car il n’est pas facile de trouver des témoins.

Par exemple, 5472940991761=199×241×863×1322335472940991761 = 199 \times 241 \times 863 \times 132233 et chacun de ces quatre facteurs est un témoin pour le test de Fermat (car an11modna^{n-1} \neq 1 \mod n si aa divise nn…).

L’identification de ces facteurs n’est cependant pas un problème simple, puisque cela suffirait à montrer la non-primalité de 54729409917615472940991761

💡 Le test de Fermat est utilisé dans les tests unitaires de CPython, comme test sommaire de la primalité du paramètre utilisé par la fonction int.__hash__.

Test de Miller-Rabin

Si p>2p > 2 est un nombre premier, on peut écrire p1=2r×dp - 1 = 2^r \times d avec r1r \geq 1 et dd impair.

D’après le petit théorème de Fermat, si 1a<p1 \leq a \lt p, alors ap1=a2rd=1modpa^{p - 1} = a^{2^r d} = 1 \mod p.

Le nombre a2r1da^{2^{r-1} d} est donc solution de l’équation x2=1modp{x^2 = 1 \mod p}, ou de façon équivalente (x+1)(x1)=0modp{(x + 1)(x - 1) = 0 \mod p}.

On en déduit, car pp est premier, que a2r1d=±1modpa^{2^{r-1} d} = \pm 1 \mod p.

Ainsi, par récurrence finie on montre que l’une des deux propositions suivantes est vraie :

  • ad=1modpa^d = 1 \mod p
  • Il existe un entier i0,r1i \in \llbracket 0, r - 1 \rrbracket tel que a2id=1modpa^{2^i d} = -1 \mod p

Le test de Miller-Rabin consiste à chercher, pour un entier nn et un témoin aa donné, si cette propriété est mise en défaut :

COMPOSITE = 'COMPOSITE'
NO_WITNESS = 'NO_WITNESS'
def miller_rabin_test(n, witnesses):
assert n >= 3
# Find greatest r s.t. 2^r | n - 1
r = 1
while (n - 1) % 2**(r + 1) == 0:
r += 1
d = (n - 1) // 2**r
for a in witnesses:
# Check first equality
x = pow(a, d, n)
if x in (1, n - 1):
continue
# Check successive squares
for i in range(1, r):
x = pow(x, 2, n)
if x == n - 1:
break
if x == 1:
return COMPOSITE
else:
return COMPOSITE
return NO_WITNESS

L’algorithme présenté ci-dessus est celui implémenté dans OpenSSL.

Il permet de discriminer rapidement les cas ci-dessous :

>>> miller_rabin_test(4208312609**2, [2])
'COMPOSITE'
>>> miller_rabin_test(5472940991761, [2])
'COMPOSITE'

Cependant, il existe toujours pour certains témoins des « faux-positifs » (on parle de pseudo-premiers forts), comme ci-dessous :

>>> miller_rabin_test(47197, [3])
'NO_WITNESS'
Témoins de Miller-Rabin

La fiabilité des tests de Fermat et de Miller-Rabin réside dans la capacité à trouver, pour un nombre composé, un témoin.

Dans le cas du test de Fermat, on a vu que pour certains nombres (par exemple ceux de Carmichael), il n’est pas garanti que l’on puisse trouver efficacement un témoin.

Dans le cas du test de Miller-Rabin, on a le résultat suivant pour un entier n>4n > 4 impair et non premier :

La proportion des entiers entre 22 et n2n-2 qui sont témoins de Miller-Rabin pour nn est supérieure à 34\frac{3}{4}.

On pourra consulter ce document par Keith Conrad pour une preuve du résultat.

Ceci nous donne un critère permettant de contrôler la probabilité de détection d’un nombre composite, lorsque les témoins sont choisis selon une loi uniforme.

Par exemple :

from random import randrange
def is_probable_prime(n, rounds=10):
if n == 1:
return False
if n in (2, 3):
return True
return miller_rabin_test(
n,
(randrange(2, n - 2) for _ in range(rounds))
) != COMPOSITE

On vérifie ensuite :

>>> is_probable_prime(4208312609)
True
>>> is_probable_prime(4208312609**2)
False

Dans le second cas, la probabilité que 420831260924208312609^2 ne soit pas détecté comme composite est inférieure à (14)10(\frac{1}{4})^{10}.

Pour des nombres sur 2048 bits, OpenSSL teste jusqu’à 128 nombres aléatoire, ce qui permet d’obtenir une probabilité d’erreur inférieure à 22562^{-256}.

💡 Complexité de PRIMES

Dans un formalisme théorique, le problème de décision PRIMES consiste à trouver un algorithme qui, pour un nombre nn représenté par bb chiffres binaires, détermine s’il est premier.

Algorithme déterministe

L’algorithme déterministe (fonction is_prime en début d’article), consistant à tester successivement tous les entiers entre 22 et nn, effectue un nombre d’opération en O(nα)O(n^\alpha) (où α\alpha est une constante qui dépendra des détails d’implémentation de la division).

Comme n2bn \leq 2^b, ce nombre d’opérations exprimé en fonction de nn est un O(2α×b)O(2^{\alpha \times b}).

Ainsi, le problème PRIMES est décidé en temps exponentiel par cet algorithme.

Algorithme randomisé

Si l’on autorise l’utilisation d’une source de nombres aléatoires, on parle d’algorithme randomisé.

Le test de Miller-Rabin présenté plus haut nous a permis de définir une fonction is_probable_prime qui renvoie en temps polynomial — je vous laisse le justifier — un résultat avec une probabilité d’erreur inférieure à 14\frac{1}{4} (pour une passe de Miller-Rabin).

Du point de vue théorique, ce résultat montre que PRIMES appartient à la classe de complexité BPP.

Primes est dans P

Des chercheurs de l’Indian Institute of Technology Kanpur ont montré en 2002, l’existence d’un algorithme déterministe qui décide de la primalité en temps polynomial en bb. Ce résultat théorique fort ne donne pas d’implémentation utilisable en pratique.

Génération de nombres premiers

Nous avons évoqué le test de la primalité d’un entier donné. Comment passer de ce problème à celui de la génération de nombre premiers d’une taille fixe ?

Une première approche consiste à générer aléatoirement des nombres impairs sur bb bits — dans l’intervalle 2b1+1,2b1\llbracket 2^{b-1} + 1, 2^b - 1 \rrbracket — jusqu’à en trouver un qui passe le test de primalité.

La question qui se pose est de savoir combien d’essais-erreurs sont nécessaires pour trouver un nombre premier avec cette procédure ?

Commençons par une approche empirique :

import statistics
def gen_prime_trials(bits):
i = 1
while True:
rnd = randrange(2**(bits - 1) + 1, 2**bits, 2)
if is_probable_prime(rnd):
return i
i += 1
print("b\tTrials\tTrials / b")
for b in range(2, 64, 5):
avg = statistics.mean(gen_prime_trials(b) for _ in range(10**4))
print(f"{b}\t{float(avg):.3}\t{avg / b:.3}")

On constate que le nombre moyen d’essais nécessaires semble être (asymptotiquement) proportionnel au nombre de bits :

b Trials Trials / b
2 1.0 0.5
7 2.46 0.351
12 3.99 0.333
17 5.76 0.339
22 7.49 0.34
27 9.16 0.339
32 11.0 0.344
37 12.5 0.339
42 14.4 0.342
47 16.2 0.345
52 17.8 0.342
57 20.0 0.351
62 21.3 0.343
Théorème des nombres premiers

Le théorème des nombres premiers est un résultat fondamental qui peut se formuler de la façon suivante :

En désignant par π(n)\pi(n) le compte des nombres premiers pnp \leq n, on a l’équivalence asymptotique : π(n)nln(n){\pi(n) \thicksim \cfrac{n}{\ln(n)}}.

Dans l’intervalle 2b1,2b1\llbracket 2^{b-1}, 2^b - 1 \rrbracket, on aura donc « de l’ordre de » 2b1ln(2b)\frac{2^{b-1}}{\ln(2^b)} nombres premiers.

Pour un tirage aléatoire parmi des nombres impairs de cet intervalle, ceci nous donne une probabilité 2b12b2ln(2b)=2bln(2)\frac{2^{b-1}}{2^{b-2} \ln(2^b)} = \frac{2}{b \ln(2)}.

En considérant l’espérance de la loi géométrique associée, on en déduit le nombre moyen d’essais :

N(b)=bln(2)20.346×bN(b) = \cfrac{b \ln(2)}{2} \approx 0.346 \times b

On explique ainsi le résultat empirique précédent.

💡 Remarque subtile

La procédure gen_prime(b) effectue des appels successifs à la fonction is_probable_prime.

Pour un tirage uniforme d’un nombre nn de taille bb, on sait borner la probabilité P(b)P(b) que is_probable_prime(n) soit un faux-positif.

Néanmoins, la probabilité d’obtenir un nombre composite avec notre procédure gen_prime(b) est supérieure à P(b)P(b) — voyez-vous pourquoi ?

L’article Average case error estimates for the strong probable prime test s’intéresse à l’estimation et au contrôle de cette probabilité.

Limitations et implémentations pratiques

La procédure décrite ci-dessus présente des limitations pratiques : pour générer un nombre premier sur 20482048 bits (minimum recommandé aujourd’hui), il faut en moyenne 0.346×20487090.346 \times 2048 \approx 709 essais.

Chacun de ces essais comporte en particulier les opérations suivantes, qui sont gourmandes en resources et en temps :

  • la génération de 20482048 bits aléatoires
  • le test de Miller-Rabin

Afin de générer des nombres premiers plus rapidement, la plupart des implémentations utilisent des heuristiques qui diminuent le nombre d’appels à ces deux opérations.

OpenSSL génère par exemple un nombre aléatoire incrémenté successivement et testé par une méthode de crible (fonction probable_prime), la librarie crypto de Go intègre également une procédure similaire.

Voici une version simplifiée, en Python, de cette heuristique :

SMALL_PRIMES = [
2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61, 67,
71, 73, 79, 83, 89, 97, 101, 103, 107, 109, 113, 127, 131, 137, 139, 149,
151, 157, 163, 167, 173, 179, 181, 191, 193, 197, 199, 211, 223, 227,
229, 233, 239, 241, 251, 257, 263, 269, 271, 277, 281, 283, 293, 307,
311, 313, 317, 331, 337, 347, 349, 353, 359, 367, 373, 379, 383, 389,
397, 401, 409, 419, 421, 431, 433, 439, 443, 449, 457, 461, 463, 467,
479, 487, 491, 499
]
def has_small_factor(n):
for p in SMALL_PRIMES:
if p**2 > n:
return False
if n % p == 0:
return True
return False
def gen_prime(bits):
while True:
rnd = randrange(2**(bits - 1) + 1, 2**bits, 2)
delta = 0
while rnd + delta < 2**bits:
if has_small_factor(rnd + delta):
delta += 2
elif is_probable_prime(rnd + delta):
return rnd + delta
else:
break

Ce programme, bien que simpliste et écrit en Python contrairement à OpenSSL (en C), permet déjà de générer en un temps raisonnable des nombres premiers sur 2048 bits :

>>> gen_prime(2048)
2001654208410281697166113376152877952059774812539457567141276111875935594130
7177853087753246463390886375884900098872655420132639282144229658911423922253
3397994790176285299341649345738609245279453160430138010182085938848124936370
9920643660264696071791799127181964339284714740892034337350851716475260350807
6643843155606155367699307063482046551001080336218108256721430098736604339134
3264153045851641530643543343048302123376397042344217184508522527154966446290
5224684588239663684021039550926018871922710560835679520722287506767855874235
5313318698844909063214755154640315114450028082650226711129793844459698307019
042102663

💡 Comme indiqué dans le Handbook of Applied Cryptography, un théorème de Mertens permet de donner une approximation de la probabilité qu’un nombre composé soit éliminé par le crible has_small_factor.

Biais d’échantillonnage

L’heuristique présentée précédemment introduit un biais potentiel : le tirage aléatoire effectué par rnd = randrange(...) est uniforme, mais la procédure incrémentale consistant à chercher parmi les nombres de la forme rnd + delta introduit un biais..

On peut facilement le constater en générant une série de petits nombres premiers :

>>> c = Counter(gen_prime(10) for _ in range(10**4))
>>> c.most_common(3)
[(907, 403), (541, 316), (967, 291)]
>>> c.most_common()[-3:]
[(601, 36), (883, 35), (619, 34)

Ci-dessus, parmi les 7575 nombres premiers entre 292^9 et 2102^{10}, on remarque que des nombres comme 601601, 883883 ou 619619 reviennent relativement peu souvent (on s’attendrait à ce qu’ils apparaissent environ 130130 fois en moyenne pour une loi uniforme).

Ces nombres sont des nombres premiers jumeaux (successeurs respectifs de 599599, 881881 et 617617).

Au contraire, 907907 est généré plus souvent — le nombre premier qui le précède, 887887, est relativement éloigné.

Ce biais, qui a été discuté sur la mailing list OpenSSL, contribue à faire baisser l’entropie effective du générateur de nombres premiers.

Pour aller plus loin…

Cet article et cette entrée de blog suggèrent que les collisions observées en pratique sont principalement dues à des défauts d’implémentation de générateurs de nombres aléatoires.

Nous avons volontairement passé sous silence ce sujet crucial, il nécessiterait un article à part entière…

Le code source des fonctions présentées est disponible ici.



Maths et applications, avec les mains et avec du code 💻
Suivre sur Twitter