The quadratic sieve, implemented: factoring RSA semiprimes in Python

MAJOR credit to NachiketUN on github for serving as the inspiration for the program itself. I couldn't have programmed it without his detailed, streamlined implementation of a similar Quadratic Sieve Algorithm.

What are RSA and the Quadratic Sieve Algorithm?

RSA has been one of the most prominent asymmetric-key (meaning one key is public and one is private) encryption algorithms since its inception in 1977 by MIT colleagues Ron Rivest, Adi Shamir, and Leonard Adleman. As we'll see, RSA is predicated on the difficulty of factoring very large semi-prime numbers. Semi-prime numbers are defined as numbers with only 4 factors: 1, itself, and two primes \(p\) and \(q\). Looking at the semi-prime \(N = 12193263122374638001\) (that's 20 digits!), it's tough to 'see' that its two prime factors are \(p = 1234567891\) and \(q = 9876543211\). A computer could exhaustively check divisibility of \(N\) by every integer from 2 to \(\sqrt{12193263122374638001} = 3.49\) billion to find \(p\) and \(q\), but there ought to be a better way...

... And there is. It's called the Quadratic Sieve (QS) Algorithm, and it's currently the second most efficient way to factor large semi-primes (behind the General Number Field Sieve, which is far more complex). Of course, once quantum computers advance far enough, Shor's Algorithm will blow it out of the water, but for now, QS remains a useful algorithm for factoring, hence my decision to code it in Python to break the RSA encryption scheme.

NOTE: RSA uses semi-primes upwards of 300 decimal digits, my QS Algorithm can only factor semi-primes of no more than 25 decimal digits in reasonable time. If I had a supercomputer, I could've managed the upper limit of QS of about 100 decimal digits, but no one is breaking \(real\) RSA encryption as it's used on the internet.

Sit back, relax, and enjoy some math and computer science!

RSA Code

p = 1234567891
q = 9876543211
N = p*q

def Totient(p, q):
    return (p-1)*(q-1)

def mod_inverse(a, mod):
    return pow(a, -1, mod)

def KeyGeneration(p, q):
    phi = Totient(p, q)
    e = 65537  # select a public key 1 < e < phi s.t. gcd(e, phi) = 1
    d = mod_inverse(e, phi)  # private key = modular inverse of e
    return e, d

def encrypt(m, e, N):
    return pow(m, e, N)

def decrypt(c, d, N):
    return pow(c, d, N)

def main():
    print("Welcome to the RSA Encryption Program!")
    print()
    message = int(input("Enter a number: "))
    public_key, private_key = KeyGeneration(p, q)
    c = encrypt(message, public_key, N)
    print(f"Encrypted Message: {c}")
    print(f"Decrypted Message: {decrypt(c, private_key, N)}")

main()

Explanation

Again, the core of RSA is the difficulty of factoring a large semi-prime \(N\), in this case, \(N = 12193263122374638001\). Since we know the prime factors of \(N\), \(p = 1234567891\) and \(q = 9876543211\), we can compute \(\phi(p, q)\) from Euler's Totient function. \(\phi(p, q)\) represents all of the integers less than \(N\) that are coprime with \(N\). Then we select a public key \(e\) for \(1 < e < \phi\) such that \(e\) is coprime with \(\phi(p, q)\). 65537 is a popular choice for \(e\) because it's prime. Encrypting an input message \(m\) in integer form entails raising \(m\) to the power of \(e\) and then taking the remainder of the resulting value modulo \(N\). We'll call this encrypted number \(c\). If we know \(p\) and \(q\), we can find the private key (used for decryption) by finding the modular inverse of \(e\) mod \(N\) such that \(ed \equiv 1 \pmod{\phi(p,q)}\), given in the code by \(\texttt{pow(a, -1, mod)}\).

Quadratic Sieve

Abstract

To factor a semi-prime \(N\), you can find \(a\) and \(b\), \(a \ne \pm b\), such that

\[a^2 - b^2 = (a-b)(a+b) \equiv 0 \pmod{N}.\]

This implies that \((a-b)(a+b) = kN\) for \(k \in \mathbb{N}\), meaning that \(\gcd(a \pm b, N)\) will yield a factor of \(N\). The chances that this factor is nontrivial are around \(\tfrac{1}{2}\).

So, we have \(a^2 - b^2 \equiv 0 \pmod{N} \rightarrow a^2 \equiv b^2 \pmod{N}\).

We can define a function \(Q(x) = (x + \lfloor\sqrt{N}\rfloor)^2 - N\) and compute \(Q(x_1), Q(x_2), \ldots, Q(x_k)\) for integers \(x_k\) over some sieving interval \([1, I]\). Ideally, we find a \(Q(x_i)\) such that

\[Q(x_i) = (x_i + \lfloor\sqrt{N}\rfloor)^2 - N = (p_1^{e_1}p_2^{e_2}\cdots p_n^{e_n})^2 \ \rightarrow\ (x_i + \lfloor\sqrt{N}\rfloor)^2 \equiv (p_1^{e_1}\cdots p_n^{e_n})^2 \pmod{N} \ \rightarrow\ a^2 \equiv b^2 \pmod{N}.\]

However, \(x_i\) that satisfy this condition are rare, so instead we typically rely on a sieving process that identifies so-called smooth \(Q(x_j)\) that factor completely over some pre-defined prime factor base, and then a composition process that finds a perfect square of prime factors modulo \(N\) for some combination of \((x_j + \lfloor\sqrt{N}\rfloor)^2\) values (each associated with a smooth \(Q(x_j)\)). In symbolic notation (for integers \(e_n\)),

\[Q(x_1)Q(x_2)\cdots Q(x_j) = (p_1^{e_1}\cdots p_n^{e_n})^2 \ \rightarrow\ \prod (x_j + \lfloor\sqrt{N}\rfloor)^2 \equiv (p_1^{e_1}\cdots p_n^{e_n})^2 \pmod{N} \ \rightarrow\ a^2 \equiv b^2 \pmod{N}.\]

Generating Factor Base and Sieving Interval

Generating a factor base over which to factor values of \(Q(x)\) is essential for QS because it forms the basis for deciding which values of \(Q(x)\) to eliminate from consideration for possible combination to produce a congruence of squares modulo \(N\).

In order to limit the complexity of the calculations to follow, we'll only consider prime factors up to a smoothness bound \(B\) in our factor base. \(B\) is arbitrary, but must be large enough to generate a sufficient amount of smooth \(Q(x)\) that are completely factorable over the factor base. Otherwise, we may not find a combination of \(Q(x)\) that yields a congruence of squares modulo \(N\).

First, we need an efficient way to generate primes up to \(B\), so we'll use the Sieve of Eratosthenes, which generates the first \(n\) primes by excluding all multiples of integers from 2 up to \(\tfrac{n}{2}\). Here's the Python implementation:

def sieve_of_eratosthenes(n):
    if n < 2:
        return []
    nums, isPrime = [], []
    for i in range(0, n+1):
        nums.append(i); isPrime.append(True)
    isPrime[0] = False; isPrime[1] = False
    for j in range(2, int(n/2)):
        if isPrime[j]:
            for i in range(2*j, n+1, j):
                isPrime[i] = False
    primes = [nums[i] for i in range(n+1) if isPrime[i]]
    return primes

We can also ignore many primes lower than the smoothness bound, though. This is because for a prime to divide \(Q(x)\),

\[(x + \lfloor\sqrt{N}\rfloor)^2 - N \equiv 0 \pmod{p} \ \rightarrow\ (x + \lfloor\sqrt{N}\rfloor)^2 \equiv N \pmod{p} \ \rightarrow\ s^2 \equiv N \pmod{p}\]

for some integer \(s\). If there is one solution to this congruence, there are two, since if \(s_1^2 \equiv N \pmod{p}\) then there must be an \(s_2 = p - s_1\) such that \(s_2^2 = (p - s_1)^2 = p^2 - 2ps_1 + s_1^2 \equiv s_1^2 \equiv N \pmod{p}\).

Therefore, only primes \(p\) for which \(N\) is a quadratic residue \(\pmod{p}\) should be included in the factor base. We can determine whether or not \(N\) is quadratic residue \(\pmod{p}\) by evaluating the Legendre symbol \(\left(\tfrac{N}{p}\right) = N^{(p-1)/2} \pmod{p}\). The following Python function checks this:

def quad_residue(N, p):
    l = 1
    q = (p-1)//2
    x = q**l
    if x == 0:
        return 1
    N = N % p
    z = 1
    while x != 0:
        if x % 2 == 0:
            N = (N**2) % p
            x //= 2
        else:
            x -= 1
            z = (z*N) % p
    return z

And the final function of this section generates the factor base:

def factor_base_gen(B):
    factor_base = []
    primes = sieve_of_eratosthenes(B)
    for p in primes:
        if quad_residue(N, p) == 1:
            factor_base.append(p)
    return factor_base

Quadratic Sieving Process

The sieving process itself involves finding values of \(Q(x_s)\) that are completely factorable over the factor base, meaning that they have no other prime factors besides those that occupy the factor base. The function below, \(\texttt{smooth\_relations}\), attempts to find these so-called smooth values of \(Q(x)\):

def smooth_relations(factor_base, N, I):
    Q_values = []
    rootN = int(sqrt(N))
    for x in range(rootN, rootN + I):
        Q_values.append(pow(x, 2) - N)

    Q_values_OG = Q_values.copy()
    if factor_base[0] == 2:
        i = 0
        while Q_values[i] % 2 != 0:
            i += 1
        for j in range(i, len(Q_values), 2):
            while Q_values[j] % 2 == 0:
                Q_values[j] //= 2

    for p in factor_base[1:]:
        residues = STonelli(N, p)
        for r in residues:
            for i in range((r - rootN) % p, len(Q_values), p):
                while Q_values[i] % p == 0:
                    Q_values[i] //= p

    xlist, smooth_nums, indices = [], [], []
    T = 1
    for i in range(len(Q_values)):
        if len(smooth_nums) >= len(factor_base) + T:
            break
        if Q_values[i] == 1:
            smooth_nums.append(Q_values_OG[i])
            xlist.append(i + rootN)
            indices.append(i)
    return smooth_nums, xlist, indices

Let's analyze this function piece by piece.

The first block creates an empty list to store values of \(Q(x)\) and evaluates \(Q(x) = (x + \lfloor\sqrt{N}\rfloor)^2 - N\) for all \(x\) in the sieving interval \([0, I]\).

The next block is the first step in determining which values of \(Q(x)\) are smooth. The code divides every even value of \(Q(x)\) in our list by 2 until the result is odd, given that 2 is in our factor base. This process effectively removes every factor of 2 from all the elements of our list that contain them. We take a copy to ensure that our manipulations of the list containing values of \(Q(x)\) don't make us lose these values.

The next segment continues the process of determining which \(Q(x)\) are smooth by dividing out odd primes \(p\) from values of \(Q(x)\) that contain them. The \(\texttt{STonelli}\) function, which is listed below, leverages the fact that all odd primes \(p\) in our factor base satisfy the condition that \(N\) is a quadratic residue modulo \(p\) to find periodic \(Q(x)\) that contain a certain prime \(p\). In symbolic notation, since \(s^2 \equiv N \pmod{p}\) has exactly two solutions for \(s\), \(s_1\) and \(s_2 = p - s_1\), then there must also be \(x_{s_1} = s_1 - \lfloor\sqrt{N}\rfloor\), \(x_{s_2} = s_2 - \lfloor\sqrt{N}\rfloor\) such that \(Q(x_{s_1} + kp) \equiv 0 \pmod{p}\) and \(Q(x_{s_2} + kp) \equiv 0 \pmod{p}\) for \(k \in \mathbb{Z}\). Once \(z\) = \(\texttt{(r-rootN) \% p}\) is one of \(x_{s_1}\) or \(x_{s_2}\) in the code, is found for each prime \(p\), we can iterate through values of \(Q(x)\), going \(p\) places each iteration, and divide by \(p\) until a given value of \(Q(x)\) is no longer divisible by \(p\).

The \(\texttt{STonelli}\) function is given below:

def STonelli(n, p):
    assert quad_residue(n, p) == 1, "not a square (mod p)"
    q = p - 1
    s = 0
    while q % 2 == 0:
        q //= 2; s += 1
    if s == 1:
        r = pow(n, (p+1)//4, p)
        return r, p-r
    for z in range(2, p):
        if p - 1 == quad_residue(z, p):
            break
    c = pow(z, q, p)
    r = pow(n, (q+1)//2, p)
    t = pow(n, q, p)
    m = s; t2 = 0
    while (t - 1) % p != 0:
        t2 = (t*t) % p
        for i in range(1, m):
            if (t2 - 1) % p == 0:
                break
            t2 = (t2*t2) % p
        b = pow(c, 1 << (m - i - 1), p)
        r = (r*b) % p
        c = (b*b) % p
        t = (t*c) % p
        m = i
    return (r, p - r)

The final segment iterates through the values of \(Q(x)\) from the copied list and determines which ones have been completely factored over the factor base, which would imply that all divisors have been divided out, resulting in a value of 1. Note too that the number of smooth \(Q(x)\) collected doesn't need to exceed the size of factor base itself past some tolerance interval \(T\), since finding a combination of \(Q(x_i)\) becomes possible as soon as the number of smooth \(Q(x)\) is the same as the number of factors in the factor base.

Now that we have all smooth values of \(Q(x)\), we either need to find a \(Q(x_i)\) that is itself a perfect square — which is unfortunately quite unlikely — or find a combination of \(Q(x)\) whose product is a perfect square modulo \(N\). The latter is doable with linear algebra.

Matrix Construction

We will construct a matrix with rows corresponding to each smooth \(Q(x)\), and columns corresponding to each factor in the prime factor base. The values occupying each position in the matrix will reflect the exponent on a given prime factor of a given smooth \(Q(x)\) modulo 2. This has the effect of depicting whether or not the exponent on a prime factor of a given \(Q(x)\) is even or odd. We want to find a combination of rows that yield 0's across the entire column space, which would imply that all exponents of all prime factors of this composite value are even, yielding a perfect square modulo \(N\). The following function \(\texttt{exponent\_matrix}\) creates said matrix with the help of a \(\texttt{factor}\) function that records each prime factor and its number of occurrences in a list for each smooth \(Q(x)\).

def factor(s, factor_base):
    factors = []
    for p in factor_base:
        while s % p == 0:
            factors.append(p); s //= p
    return factors

def exponent_matrix(smooth_nums, factor_base):
    A = []
    for s in smooth_nums:
        exp_vector = [0]*(len(factor_base))
        s_factors = factor(s, factor_base)
        for i in range(len(factor_base)):
            if factor_base[i] in s_factors:
                exp_vector[i] = s_factors.count(factor_base[i]) % 2
        if 1 not in exp_vector:
            return True, s
        A.append(exp_vector)
    return [False], transpose(A)

NOTES: (1) returns \(\texttt{True}\) if a column with all 0's is found (in which case we have a solution), (2) Transposition of matrix \(A\) is done to aid calculations in the next step, Gaussian Elimination.

Gaussian Elimination

We have our exponent matrix modulo 2, now we need to solve for a combination of rows corresponding to smooth values of \(Q(x)\) that, when multiplied, yield a solution. Gaussian Elimination over GF(2) (i.e. working modulo 2) does this by first identifying pivot elements in the matrix that can be used to eliminate variables and reduce the system to a simpler form. This step is crucial as it helps us isolate the rows (i.e., smooth numbers) that can be combined to potentially form a perfect square.

def gaussian_elimination(A):
    marks = [False]*len(A[0])
    for i in range(len(A)):
        row = A[i]
        for num in row:
            if num == 1:
                j = row.index(num)
                marks[j] = True
                for k in chain(range(0, i), range(i+1, len(A))):
                    if A[k][j] == 1:
                        for i in range(len(A[k])):
                            A[k][i] = (A[k][i] + row[i]) % 2
                break
    A = transpose(A)
    sol_rows = []
    for i in range(len(marks)):
        if marks[i] == False:
            free_row = [A[i], i]
            sol_rows.append(free_row)
    if not sol_rows:
        return "No solution found. Need more smooth numbers."
    print("Found {} potential solutions".format(len(sol_rows)))
    return sol_rows, marks, A

Explanation:

A boolean array \(\texttt{marks}\) of length equal to the number of columns in matrix \(A\) is initialized to \(\texttt{False}\). This array tracks which columns have been used to eliminate other 1s in the same column, effectively marking the pivot columns.
For each row \(i\) in \(A\), we search for the first occurrence of 1, which serves as our pivot element. Upon finding a pivot at column \(j\), \(\texttt{marks[j]}\) is set to \(\texttt{True}\), signifying that column \(j\) has been used for elimination.
Using the pivot element found, the algorithm then proceeds to eliminate other 1s in the same column across different rows. For any row \(k\) with a 1 in the pivot column, we add (modulo 2) the pivot row to row \(k\), which effectively zeroes out the 1 in row \(k\)'s pivot column.
After processing all rows, matrix \(A\) is transposed. This step is vital as it facilitates the identification of free variables (columns that were not marked), which after transposition, become solution rows indicating combinations of smooth numbers.
If there are no solution rows, more smooth numbers are needed to populate matrix rows.

Following the Gaussian Elimination and the transposition of \(A\), we now look for columns (which are originally rows before transposition) that were not marked by any pivot. These unmarked columns, now rows in the transposed matrix, represent the free variables in our system. Each free variable offers a potential solution vector, which is a combination of smooth numbers that, when multiplied together, might yield a perfect square modulo \(N\).

Solving for Rows with the solve_row Function

After Gaussian Elimination and identifying potential solution vectors, the \(\texttt{solve\_row}\) function isolates specific combinations of smooth numbers that may lead to a successful factorization of \(N\). This function takes the solution rows identified earlier and pinpoints which rows in the original matrix correspond to a viable solution vector.

def solve_row(sol_rows, A, marks, K=0):
    solution_vec, indices = [], []
    free_row = sol_rows[K][0]
    for i in range(len(free_row)):
        if free_row[i] == 1:
            indices.append(i)
    for r in range(len(A)):
        for i in indices:
            if A[r][i] == 1 and marks[r]:
                solution_vec.append(r); break
    solution_vec.append(sol_rows[K][1])
    return solution_vec

The \(\texttt{solve\_row}\) function operates as follows:

It first selects a "solution row" based on an index \(K\). This row corresponds to a free variable in the transposed matrix, indicating a potential combination of original rows (smooth numbers) that sum to zero modulo 2.
The function then extracts indices within the selected solution row where the value is 1. These indices represent the original matrix's rows that, when combined, could yield a perfect square.
With the indices identified, \(\texttt{solve\_row}\) scans the entire matrix to find rows that share these indices. It constructs a solution vector comprising rows that contain a 1 in any of the marked positions. This vector signifies a group of smooth numbers that, collectively, are likely to produce a square when multiplied together.
The index of the free variable itself is also appended to the solution vector, ensuring that the corresponding smooth number is included in the final product.

With our solution vector, we perform one final step to find the factors of \(N\).

Computing Factors From Solution Vector

This is the final step in the sieving process! And don't worry, it's not too hard. The task is left to the \(\texttt{solve}\) function, which aims to find \(a = (p_1^{e_1}\cdots p_n^{e_n})\) and \(b = \prod (x_j + \lfloor\sqrt{N}\rfloor)\) in the congruence \(a^2 \equiv b^2 \pmod{N}\). This implies that \((a-b)(a+b) = kN\), meaning that \(\gcd(a \pm b, N)\) will yield a factor of \(N\).

def solve(solution_vec, smooth_nums, xlist, N):
    solution_nums = [smooth_nums[i] for i in solution_vec]
    x_nums = [xlist[i] for i in solution_vec]

    square = 1
    for n in solution_nums:
        square *= n
    a = int(sqrt(square))

    b = 1
    for n in x_nums:
        b *= n

    factor = gcd(b - a, N)
    return factor

With our solution vector, we know exactly which smooth values of \(Q(x)\) should be multiplied to yield a perfect square modulo \(N\). We extract these values and their corresponding \(x + \lfloor\sqrt{N}\rfloor\) values, take the square root of the product to get \(a\), multiply the \(x\)-values to get \(b\), and compute \(\gcd(b - a, N)\) (yields same factor as \(\gcd(a-b, N)\)), yielding a (hopefully non-trivial) factor of \(N\).

Breaking RSA

Now that we have an algorithm to factor large semi-primes, let's test it on \(N = 12193263122374638001\). Here are the terminal logs:

abelobsenz@Abes-MacBook-Air desktop % python3 QuadraticSieveAlgorithm.py
Enter a semiprime number:  12193263122374638001
Enter an appropriate smoothness bound:  100000
Enter an appropriate interval for smoothness testing (should be F + 1 smooth numbers):  100000
Factoring 12193263122374638001 now
Generating B-smooth factor base...
Looking for smooth numbers...
Building exponent matrix...
Unlucky, you need Gaussian elimination...
Found 1407 potential solutions
Solving congruence of squares...
Didn't work. Trying different solution vector...
Didn't work. Trying different solution vector...
Found factors!
Factors of N: (9876543211, 1234567891)

What to do with these factors, what to do? I think we should find the private key in RSA by computing the inverse of the public key (65537) mod \(\phi(9876543211, 1234567891)\)!

It's 6751626181311963136. He he, time to cause some trouble.