Programming Praxis

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)² – 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 g_{a, b}(x) = (a x + b)² − n with a, b integers with 0 < b ≤ a. The graph of g_{a, b}(x) is a parabola, and its values will be smallest when a ≈ √(2n) / m. Thus, we choose b so that b² − n is divisible by a, say b² − n = a c for some integer c, and a = q² for some integer q. Then we calculate g_{a, b}(x) / a = a x² + 2 b x + c, and, after sieving over the range −m .. m, when g_{a, b}(x) is smooth over the factor base we record the relation ((a x + b) q⁻¹)² = a x² + 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 t² ≡ 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 = ⌊log₂ 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 = q² (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 q² using Hensel’s Lemma. For each odd prime p in the factor base, compute soln1_p = a⁻¹ (tmem_p − b) and soln2_p = a⁻¹ (−tmem_p − 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 l_p to the locations soln1_p + i p for all integers i that satisfy −m ≤ soln1_p + i p ≤ m, and likewise for soln2_p. For the prime p = 2, sieve only with soln1_p.

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 g_{a, b}(x) by the primes in the factor base. If g_{a, 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 x² ≡ y² (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.

Pages: 1 2 3