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 > 1$ qui n’ont que deux diviseurs : $1$ et $n$ lui-même.

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

Étant donné un nombre entier positif $n$ , $n$ 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 $p$ est un nombre premier et $a$ un entier quelconque, alors $a^p - a$ est un multiple de $p$ .

Exercice :

En remarquant que $a^p$ correspond au nombre de mots de $p$ lettres que l’on peut composer avec un alphabet de $a$ 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 $1 \leq a \lt p$ , montrer que l’égalité du théorème se simplifie en $a^{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 $5472940991761$ (identifié dans cet article) est un nombre de Carmichael : ce sont les nombres $n$ qui ne sont pas premiers mais vérifient $a^{n - 1} = 1 \mod n$ pour tout entier $a$ premier avec $n$ .

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 \times 241 \times 863 \times 132233$ et chacun de ces quatre facteurs est un témoin pour le test de Fermat (car $a^{n-1} \neq 1 \mod n$ si $a$ divise $n$ …).

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

💡 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 > 2$ est un nombre premier, on peut écrire $p - 1 = 2^r \times d$ avec $r \geq 1$ et $d$ impair.

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

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

On en déduit, car $p$ est premier, que $a^{2^{r-1} d} = \pm 1 \mod p$ .

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

$a^d = 1 \mod p$
Il existe un entier $i \in \llbracket 0, r - 1 \rrbracket$ tel que $a^{2^i d} = -1 \mod p$

Le test de Miller-Rabin consiste à chercher, pour un entier $n$ et un témoin $a$ 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 > 4$ impair et non premier :

La proportion des entiers entre $2$ et $n-2$ qui sont témoins de Miller-Rabin pour $n$ est supérieure à $\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 $4208312609^2$ ne soit pas détecté comme composite est inférieure à $(\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 à $2^{-256}$ .

💡 Complexité de PRIMES

Dans un formalisme théorique, le problème de décision PRIMES consiste à trouver un algorithme qui, pour un nombre $n$ représenté par $b$ 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 $2$ et $n$ , effectue un nombre d’opération en $O(n^\alpha)$ (où $\alpha$ est une constante qui dépendra des détails d’implémentation de la division).

Comme $n \leq 2^b$ , ce nombre d’opérations exprimé en fonction de $n$ est un $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 à $\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 $b$ . 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 $b$ bits — dans l’intervalle $\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 $\pi(n)$ le compte des nombres premiers $p \leq n$ , on a l’équivalence asymptotique : ${\pi(n) \thicksim \cfrac{n}{\ln(n)}}$ .

Dans l’intervalle $\llbracket 2^{b-1}, 2^b - 1 \rrbracket$ , on aura donc « de l’ordre de » $\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é $\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) = \cfrac{b \ln(2)}{2} \approx 0.346 \times b$

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

💡 Une petite subtilité

Pour un tirage uniforme d’un nombre $n$ de taille $b$ , on a montré qu’on pouvait contrôler la probabilité $P(b)$ que is_probable_prime(n) donne un faux-positif.

La procédure gen_prime_trials(b) effectue des appels successifs à la fonction is_probable_prime, la probabilité d’obtenir un nombre composite avec cette procédure est donc supérieure à $P(b)$ .

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 $2048$ bits (minimum recommandé aujourd’hui), il faut en moyenne $0.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 $2048$ 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 $75$ nombres premiers entre $2^9$ et $2^{10}$ , on remarque que des nombres comme $601$ , $883$ ou $619$ reviennent relativement peu souvent (on s’attendrait à ce qu’ils apparaissent environ $130$ fois en moyenne pour une loi uniforme).

Ces nombres sont des nombres premiers jumeaux (successeurs respectifs de $599$ , $881$ et $617$ ).

Au contraire, $907$ est généré plus souvent — le nombre premier qui le précède, $887$ , 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