June 21, 2013
We studied the basic quadratic sieve algorithm for factoring integers in two previous exercises. Today we examine the multiple polynomial variant of the quadratic sieve, due to Peter Montgomery, which adds considerable power to the basic algorithm; we will be following Figure 2.2 and the surrounding text from Scott Contini’s thesis.
In the basic quadratic sieve we used a polynomial of the form g(x) = (x + b)2 – n where b = ⌈√n⌉. The problem with using a single polynomial is that the values of g(x) increase as x increases, making them less likely to be smooth. Eventually, the algorithm “runs out of gas” when the values of g(x) grow too large.
Montgomery’s suggestion was to use multiple polynomials of the form ga, b(x) = (a x + b)2 − n with a, b integers with 0 < b ≤ a. The graph of ga, b(x) is a parabola, and its values will be smallest when a ≈ √(2n) / m. Thus, we choose b so that b2 − n is divisible by a, say b2 − n = a c for some integer c, and a = q2 for some integer q. Then we calculate ga, b(x) / a = a x2 + 2 b x + c, and, after sieving over the range −m .. m, when ga, b(x) is smooth over the factor base we record the relation ((a x + b) q−1)2 = a x2 + 2 b x + c.
Here is our version of Contini’s algorithm:
Compute startup data: Choose f, the upper bound for factor base primes, and m, which is half the size of the sieve. Determine the factor base primes p < f such that the jacobi symbol (n/p) is 1 indicating that there is a solution to t2 ≡ n (mod p); also include 2 in the factor base. For each factor base prime p, compute and store t, a modular square root of n (mod p). Also compute and store the base-2 logarithm of each prime l = ⌊log2 p⌉ (the floor and ceiling brackets indicate rounding).
Initialize a new polynomial: Find a prime q ≈ √( √(2n) / m) such that the jacobi symbol (n/q) is 1, indicating that n is a quadratic residue mod q, and let a = q2 (mod n). Compute b, a modular square root of n mod a; you will have to compute the square root of n mod q then “lift” the root mod q2 using Hensel’s Lemma. For each odd prime p in the factor base, compute soln1p = a−1 (tmemp − b) and soln2p = a−1 (−tmemp − b).
Perform sieving: Initialize a sieve array of length 2 m + 1 to zeros, with indices from −m to m. For each odd prime p in the factor base, add lp to the locations soln1p + i p for all integers i that satisfy −m ≤ soln1p + i p ≤ m, and likewise for soln2p. For the prime p = 2, sieve only with soln1p.
Trial division: Scan sieve array for locations x that have accumulated a value of a least log(m √n) minus a small error term. Trial divide ga, b(x) by the primes in the factor base. If ga, b(x) factors into primes less than f, then save smooth relation as indicated above. After scanning entire sieve array, if there are more smooth relations than primes in the factor base, then go to linear algebra step. Otherwise, go to initialization step to select a new polynomial.
Linear algebra: Solve the matrix as in Dixon’s method and the continued fraction method. For each null vector, construct relation of form x2 ≡ y2 (mod n) and attempt to factor n by computing gcd(x − y, n). If all null vectors fail to give a non-trivial factorization, go to initialization step to select a new polynomial.
Your task is to write a program to factor integers using the multiple polynomial quadratic sieve as described above. When you are finished, you are welcome to read or run a suggested solution, or to post your own solution or discuss the exercise in the comments below.