Séance 8 : factorisation d’entiers – crible quadratique

Contenu

Séance 8 : factorisation d’entiers – crible quadratique#

Exercice#

Dans cet exercice, on propose d’implémenter (de façon naïve et non-optimisée) l’algorithme de crible quadratique qui trouve un diviseur propre d’un entier \(N\).

Question 1. Écrire une fonction calcule_base_facteurs(N, B) qui calcule une base de facteurs plus petits que B dans le but de factoriser un entier N.

def calcule_base_facteurs(N, B):
    L = [ 2 ]
    p = 3
    while p <= B:
        if jacobi_symbol(N, p) == 1:
            L.append(p)
        p = next_prime(p)
    return L
    
print(calcule_base_facteurs(N=369713, B=21))

[2, 7, 11, 19]

Question 2. Écrire une fonction de quad_sieve(FB, N, A) qui prend en entrée une base de facteurs FB, un entier N à factoriser et une borne A, et qui effectue l’étape d’effritement (crible).

def quad_sieve(FB, N, A): 
    S = ceil(sqrt(N))
    Q = lambda x: x**2 + 2*S*x + S**2 - N
    T = [ Q(x) for x in range(A) ]
    for p in FB:
        a = mod(N, p)
        if a.is_square():
            roots = a.sqrt(all=True)
            roots = [Integer(r-S) for r in roots ]
            roots = [r for r in roots if r < A ]
            e = 1
            while roots != []:
                for x in roots:
                    while x < A:
                        T[x] = T[x]//p
                        x = x + p**e
                e += 1
                a = mod(N, p**e)
                roots = a.sqrt(all=True)
                roots = [Integer(r-S) for r in roots ]
                roots = [r for r in roots if r < A ]
    return [ [x,Q(x)] for x in range(A) if T[x] == 1 ]
    
N = 369713
FB = calcule_base_facteurs(N=369713, B=21)
A = ceil(sqrt(N))//3
L = quad_sieve(FB, N, A)
print(L)

[[6, 8512], [8, 10976], [24, 30976], [106, 141512], [120, 161728]]

Question 3. Écrire une fonction de build_matrix(FB, sieve) qui prend en entrée une base de facteurs FB et une liste d’entiers effrités sieve, et qui construit la matrice associée.

def build_matrix(FB, res_sieve):
    n = len(FB)
    m = len(res_sieve)
    M = zero_matrix(ZZ, m, n)
    for i in range(m):
        r = res_sieve[i][1]
        for j in range(n):
            p = FB[j]
            while (r != 0) and ((r % p) == 0):
                M[i,j] = M[i,j] + 1
                r = r//p
    return M
    
M = build_matrix(FB, L)
print(M)

[6 1 0 1]
[5 3 0 0]
[8 0 2 0]
[3 2 0 2]
[6 1 0 2]

Question 4. Assembler les fonctions pour obtenir une fonction qui, étant donné un entier N en entrée, retourne un diviseur de N grâce à l’algorithme du crible quadratrique.

def quadratic_sieve(N, verbose=False):

    # Initialisation des constantes
    sqN = ceil(sqrt(N))
    Q = lambda x: x**2 + 2*sqN*x + sqN**2 - N
    A = sqN//3 # A = sqN
    B = floor(2**(0.5*sqrt(log(N, 2) * log(log(N, 2), 2))))
    B_large_enough = False
    
    # Collecte de relations (en avoir + que de facteurs)
    #  et solution du système 
    while not(B_large_enough):
        FB = calcule_base_facteurs(N, B)
        s = len(FB)
        sieve = quad_sieve(FB, N, A)
        M = build_matrix(FB, sieve)
        M2 = matrix(GF(2), M)
        B_large_enough = (M2.kernel().dimension() > 0)
        if not(B_large_enough) and verbose:
            print("on augmente B")
            B *= 2
        
    # Reconstruction d'un diviseur
    factor_found = False
    while not(factor_found):
        u = M2.left_kernel().random_element()
        a = 1
        b = 1 
        vec = vector(ZZ, [0 for j in range(s)])
        for i in range(len(u)):
            if u[i] == 1:
                a = (a * (sieve[i][0] + sqN)) % N
                vec += M[i]
        vec = vec/2
        for j in range(s):
            b = (b * FB[j]**vec[j]) % N
        p, q = gcd(a-b, N), gcd(a+b, N)
        factor_found = (p not in  [1, N]) or (q not in  [1, N])
        if not(factor_found) and verbose:
            print("facteur non trouvé, on réessaie")
    return p, q

Question 5. Testez votre fonction avec \(N =\) par exemple.

# N1 = 73217
N2 = 369713
# N3 = 7029871

# print(quadratic_sieve(N1))
print(quadratic_sieve(N2))
# print(quadratic_sieve(N3))

(457, 809)
Contenu