Outils pour utilisateurs

Outils du site


factorisation_shor

Une curiosité: factorisation par l'algorithme de shor avec un ordinateur non-quantique!

Problématique

La factorisation d'un grand nombre entier dans un temps “raisonnable” est une chose difficile, mais qui comporte un énorme enjeu: celui qui trouvera un algorithme permettant de faire ça rapidement rendra obsolète tout ce qui est cryptage de type RSA (clé publique / clé privée)!

Dans mes recherches d'amateur sur le sujet, j'ai étudié par curiosité l'algorithme de Shor, considéré comme très dangereux pour l'avenir du RSA, à part qu'il faut un ordinateur quantique pour le faire tourner…

Alors, j'ai essayé avec un “ordinateur normal” (non-quantique) pour voir… Et ça marche, mais pour l'instant sans la performance qu'il faudrait.

Liens de référence (il y en a beaucoup d'autres):

Vous pouvez même trouver ici les travaux originaux de Shor lui-même (1995/96):

Voilà le début de mon étude.

Code de base

def rshor(a, P):
    """algorithme de shor: calcul de la période de (a**k)%P"""
    k = 0
    f = 1
    while True:
        k += 1
        f = pow(a, k, P)
        if f == 1:
            break
    return k  # période trouvée
 
def shor(P):
    """algorithme de shor: factorisation de P"""
    while True:
        a = random.randint(2, P-1)
        print u"=====> Essai avec a=", a
 
        x = pgcd(a,P)
        if x!=1:
            n1 = x
            n2 = P//n1
            print u"Solution: n1=", n1, "n2=", n2
            return [n1, n2]
        else:
            r = rshor(a, P)  # calcul de la periode de f(k) = (a**k)%P par la fonction separee rshor()
            print u"période r=", r
 
            if r&1!=1:
                print u"calcul de a**(r//2)"
                x = pow(a, r//2)
                if x%P!=1:
                    print u"calcul de n1"
                    n1 = pgcd(x+1, P)  #n1 = pgcd(a**(r//2)+1, P)
                    print u"calcul de n2"
                    n2 = pgcd(x-1, P)  #n2 = pgcd(a**(r//2)-1, P)
                    if (n1!=1) and (n2!=1):
                        print "Solution: n1=", n1, "n2=", n2
                        return [n1, n2]
            # si on arrive ici, on boucle pour essayer un autre "a"

Il est nécessaire d'avoir aussi le calcul du PGCD (merci Euclide!):

def pgcd(a,b):
    """pgcd(a,b): calcul du 'PGCD' entre les 2 nombres entiers a et b"""
    a, b = abs(a), abs(b)
    while b:
        r = a%b
        a, b = b, r
    return a

Il est aussi nécessaire d'avoir le calcul de l'exponentiation modulaire (a**b)%n. Avec Python v2.5, la fonction intégrée pow(a, b, n) fait ça très bien, et si elle fonctionne chez vous, il n'y a pas de raison de la remplacer. Si vous avez une version plus ancienne, vous pouvez utiliser la fonction suivante:

def lpowmod(a, b, n):
    """exponentiation modulaire: (a**b)%n"""
    r = 1
    while b>0:
        if b&1==0:
            b = b>>1
        else:
            r = (r*a)%n
            b = (b-1)>>1
        a = (a*a)%n
    return r

Résultats et Performance

Commençons par un exemple simple, celui utilisé comme test dans l'ordinateur quantique: factorisation de P = 15. On devrait trouver en principe 3 x 5.

L'algorithme tire au hasard une valeur de a entre 2 et P-1 (bornes comprises), par exemple a = 7

Si pgcd(a,P) différent de 1, on a trouvé une solution, et on arrête. Mais ce n'est pas le cas ici: pgcd(7,15)=1

On va donc calculer f = (a**k)%P, c'est à dire f = (7**k)%15, pour des valeurs de k croissantes:

k=0 f=1
k=1 f=7
k=2 f=4
k=3 f=13
k=4 f=1
k=5 f=7
k=6 f=4
k=7 f=13
k=8 f=1
k=9 f=7
# etc...

On remarque que les valeurs se retrouvent à partir de k=4: la fonction f = (a**k)%P est périodique!

Le jeu va consister à trouver la période r de cette fonction, et c'est la fonction python “rshor(a,P) qui va s'en charger.

Puisque la fonction f = (a**k)%P est égale à 1 pour k = 0, on va chercher le 1er k>0 qui va de nouveau aboutir à f=1: on aura trouvé la période r = k. Ici, on trouve la période r = 4.

On a ensuite quelques test pour vérifier que cette période, associée à la valeur de a, conduira bien à une solution non triviale.

Et si oui, les facteurs sont obtenus par:

n1 = pgcd(a**(r//2)+1, P)
n2 = pgcd(a**(r//2)-1, P)

c'est à dire ici:

print pgcd(7**(4//2)+1, 15)
5
print pgcd(7**(4//2)-1, 15)
3

Ce sont bien les 2 facteurs attendus: 5×3=15

Testons maintenant le code avec une autre valeur plus importante: P = 31313

Nombre à factoriser:
P= 31313 (produit de  181  et de  173)
=====> Essai avec a= 5811
période r= 172
calcul de a**(r//2)
calcul de n1
calcul de n2
=====> Essai avec a= 3859
période r= 7740
calcul de a**(r//2)
calcul de n1
calcul de n2
=====> Essai avec a= 30757
période r= 3870
calcul de a**(r//2)
calcul de n1
calcul de n2
Solution: n1= 181 n2= 173
181 173 31313 31313 0.0185447902482

La fonction a échoué sur une première valeur de a = 5811, puis sur a = 3859 et a réussi avec a = 30757 (au bout de 2/100 ème de seconde).

Cette dernière valeur de “a” a permis de trouver une période r=3870, qui a fournit les 2 solutions n1= 181 et n2= 173 (dont le produit est bien 31313).

Si on recommence avec d'autres valeurs, on voit que le temps d'exécution change pas mal en fonction de la pertinence du choix du “a”. Quelquefois, la réponse est quasi immédiate, d'autre fois, cela peut dépasser la dizaine de seconde pour des valeurs de cette taille.

Essayez par exemple le code d'essai suivant:

P1 = random.randint(1000,9999)
P2 = random.randint(1000,9999)
P = P1*P2
 
tps = time.clock()
x = shor(P)
tps = time.clock()-tps
print x[0], x[1], x[0]*x[1], P, tps

Si on essaye des valeurs plus grandes, le temps d'exécution augmente beaucoup trop pour rendre cet algorithme utilisable dans l'état.

En fait, la première chose qui coince, c'est le calcul de la période r = rshor(a,P). Et c'est la partie dévolue normalement à l'ordinateur quantique: ce n'est pas par hasard…

On verra aussi qu'une deuxième chose consomme beaucoup de temps pour les grandes valeurs de a et de r: le calcul de a**(r//2)

Calcul amélioré de la période r

On va seulement modifier la fonction rshor(a,P) pour essayer de trouver la période r plus rapidement pour des nombres à factoriser plus grands.

Voilà la nouvelle fonction pour le calcul de la période r:

def rshor(a, P):
    """recherche de la période r de la fonction f=pow(a,k,P)"""
 
    # calcul des premières valeurs de f(a,k,P) avec k de 1 à imax-1:
    M = [1] # car on sait que pow(a,0,P) = 1
    imax = 10000  # taille du motif de départ M
    for k in xrange(1, imax):
        f = pow(a, k, P)
        if f == 1:
            return k # on a trouvé la période r=k
        M.append(f)
 
    # définir un pas
    pas = (P-imax)//100000
    if pas<1:
        pas = 1
 
    k = imax
    while k<P:
        f = pow(a, k, P) # on calcule la valeur de f correspondante
        i = -1
        while True:
            try:
                i = M.index(f,i+1) # on cherche si f(k) se trouve dans le motif de départ
            except:
                break  # non, on passe au k suivant
 
            if pow(a,k-i,P)==1:
                r = k-i  # on a trouvé la période ou un multiple de la période
 
                # tentative de trouver un r plus petit si c'est un multiple de la période
                for j in xrange(10, 1, -1):
                    if r%j==0:
                        k = r//j
                        if pow(a,k,P)==1:
                            return k
 
                return r
        k += pas
    return 0  # c'est une condition d'échec: on n'a pas trouvé de r qui convienne

Le principe est assez simple:

on commence par calculer les 10000 premières valeurs de la fonction f = pow(a,k,P) pour k=0 à (imax-1)=9999.

Bien sûr, si on trouve un f=1 avant d'arriver au bout, on renvoie la valeur de r=k car on l'a trouvé!

Si ce n'est pas le cas, on va calculer la fonction pour un certain nombre de k commençant à imax, se terminant à P, par pas de “pas”. Ce pas est fixé a priori à (P-imax)//10000, mais il ne doit pas être inférieur à 1.

A chaque valeur de k ainsi calculé, on vérifie si le f(k) correspondant ne se trouverait pas dans le motif de départ M[]. S'il s'y trouve à l'indice i, on peut en déduire que la période est probablement à l'indice k-i. On vérifie que f(k-i)==1, et si ce n'est pas le cas, on cherche s'il ne serait pas plus loin dans le motif. Si ce n'est pas le cas, on passe au k suivant.

Petite tentative d'amélioration, si par hasard on avait trouvé non pas la période mais un de ses multiples, on cherche à trouver la période en essayant r//10, r//9, r//8, etc… jusqu'à r//2.

En ce qui concerne les résultats, ça s'améliore un peu, mais ce n'est pas encore ça.

Par exemple:

Nombre à factoriser:
P= 4948277 (produit de  2087  et de  2371 )
=====> Essai avec a= 4490746
période r= 2471910
calcul de a**(r//2)
calcul de n1
calcul de n2
=====> Essai avec a= 1457519
période r= 823970
calcul de a**(r//2)
calcul de n1
calcul de n2
Solution: n1= 2087 n2= 2371
2087 2371 4948277 102.754698033

On a bien trouvé le résultat, mais en 102 secondes.

Autre exemple:

Nombre à factoriser:
P= 5611561 (produit de  2621  et de  2141 )
=====> Essai avec a= 3948837
période r= 280340
calcul de a**(r//2)
calcul de n1
calcul de n2
Solution: n1= 2621 n2= 2141
2621 2141 5611561 4.17664580239

Là, le bon résultat [2621, 2141] a été trouvé en 4 secondes avec un seul essai a=3948837 qui a donné la période r=280340.

Les limites “raisonnables” (quelques minutes?) correspondent à un produit de nombres à 12 ou 16 bits, ce qui est déjà pas mal.

Et, curieusement, avec des nombres de cette taille, le calcul de la période se passe plutôt bien, et c'est le calcul de a**(r//2) qui dure trop de temps!

On est donc encore très loin de factoriser un nombre de 2048 bits (composé d'environ 600 chiffres)…

Manifestement cette méthode nécessite un accélérateur à 2 endroits précis:

  • le calcul de la période r
  • le calcul du a**(r//2)

En ce qui concerne le calcul de la période, si on n'a pas d'ordinateur quantique sous la main (:-D), peut-être qu'une FFT appliquée à f(k) serait une bonne solution?

En ce qui concerne le calcul de a**(r//2), il faut se rappeler que ça se fait avec des nombres très grands, et que le résultat est énorme. Pour prendre le dernier exemple, essayer donc de calculer 3948837**(280340//2), et imaginez un peu ce que ça donnerait pour des nombres de plusieurs centaines de chiffres…


Amusez-vous bien!

factorisation_shor.txt · Dernière modification: 2009/01/06 08:20 de tyrtamos

Outils de la page