## Multiple Polynomial Quadratic Sieve

### 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 *g*_{a, b}(*x*) = (*a x* + *b*)^{2} − *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* ≈ √(2*n*) / *m*. Thus, we choose *b* so that *b*^{2} − *n* is divisible by *a*, say *b*^{2} − *n* = *a c* for some integer *c*, and *a* = *q*^{2} for some integer *q*. Then we calculate *g*_{a, b}(*x*) / *a* = *a x*^{2} + 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*^{−1})^{2} = *a x*^{2} + 2 *b x* + *c*.

Here is our version of Contini’s algorithm:

Compute startup data: Choosef, the upper bound for factor base primes, andm, which is half the size of the sieve. Determine the factor base primesp<fsuch that the jacobi symbol (n/p) is 1 indicating that there is a solution tot^{2}≡n(modp); also include 2 in the factor base. For each factor base primep, compute and storet, a modular square root ofn(modp). Also compute and store the base-2 logarithm of each primel= ⌊log_{2}p⌉ (the floor and ceiling brackets indicate rounding).

Initialize a new polynomial: Find a primeq≈ √( √(2n) /m) such that the jacobi symbol (n/q) is 1, indicating thatnis a quadratic residue modq, and leta=q^{2}(modn). Computeb, a modular square root ofnmoda; you will have to compute the square root ofnmodqthen “lift” the root modq^{2}using Hensel’s Lemma. For each odd primepin the factor base, computesoln1=_{p}a^{−1}(tmem−_{p}b) andsoln2=_{p}a^{−1}(−tmem−_{p}b).

Perform sieving: Initialize a sieve array of length 2m+ 1 to zeros, with indices from −mtom. For each odd primepin the factor base, addlto the locations_{p}soln1+_{p}i pfor all integersithat satisfy −m ≤soln1+_{p}i p≤m, and likewise forsoln2. For the prime_{p}p= 2, sieve only withsoln1._{p}

Trial division: Scan sieve array for locationsxthat have accumulated a value of a least log(m √n) minus a small error term. Trial divideg_{a, b}(x) by the primes in the factor base. Ifg_{a, b}(x) factors into primes less thanf, 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 formx^{2}≡y^{2}(modn) and attempt to factornby 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.