# Multiple polynomial quadratic sieve

The **multiple polynomial quadratic sieve (MPQS)** is a factorization method.

## Contents

## Introduction

Let [math]N[/math] be the number to be factored. This number must not be a perfect power. If somehow we find two integers [math]X[/math] and [math]Y[/math] such that [math]X^2\equiv Y^2\,\pmod N[/math] and [math]X\not\equiv \pm Y\pmod N[/math], then [math]\gcd(X+Y,\, N)[/math] will reveal a proper factor of [math]N[/math].

In order to find these values [math]X[/math] and [math]Y[/math] the method finds *relations* which have the form [math]t^2 \equiv u\,\pmod N[/math] where [math]u[/math] is the product of small prime numbers. The set of these primes is the *factor base*. These relations will be found using sieving as will be explained below.

The relations found are combined using multiplications. The left hand side will always be a square because is a product of squares, so the goal is to have a square at the right hand side. A number is a square when all its prime factors appear an even number of times.

For example: let the number to factor be N = 1817 and we have found the following relations with factor base = {2, 7, 13}:

- [math]45^2\,\equiv \,2^4*7^0*13^1[/math]
- [math]123^2\,\equiv \,2^{10}*7^0*13^1[/math]

Both relations have non-square RHS, since the exponent of 13 is not even. But multiplying them we get:

- [math]84^2\,\equiv \,2^{14}*13^2[/math]
- [math]84^2\,\equiv \,(2^7*13)^2[/math]

Since [math]2^7*13\,=\,1664[/math] we get the factor [math]\gcd(84+1664,\,1817) = 23[/math].

Which relations have to multiplied to find a square in the RHS is a linear algebra problem.

## Sieving

We will sieve in the interval [-M, +M] using several polynomials [math]g_{a,b}(x)\, =\, (ax + b)^2 - N[/math] where [math]b^2 - n\, =\, ac[/math] so [math]g_{a,b}(x)\, =\,a(ax^2+bx+c)[/math] with different values of [math]a[/math] and [math]b[/math]. The advantage here is that we know the factorization of [math]a[/math], so the other factor is small compared to [math]N[/math], thus increasing the probability that all its prime factors are included in the factor base. Using the notation of the previous section, [math]t^2\equiv u\,\pmod N[/math] where [math]t = ax+b[/math] and [math]u = g_{a,b}(x)[/math].

### Generating the factor base

In order to determine the factor base, we have to notice that it includes all the possible factors of [math]g_{a,b}(x)[/math]. Since [math]g_{a,b}(x) = t^2 - N[/math], we get that [math]0 \equiv t^2 - N[/math] or [math]t^2 \equiv N[/math] modulo any prime of the factor base. This means that about half of the primes cannot be in the factor base because [math]N[/math] cannot be a quadratic residue modulo that prime. Notice that the content of this paragraph does not depend on the particular values of [math]a[/math] or [math]b[/math], so the factor base can be computed at the beginning of MPQS. The upper bound for the factor base is determined from experimentation, for example:

- Factoring a 60 digit number: F = 60000
- Factoring a 70 digit number: F = 350000
- Factoring a 80 digit number: F = 900000

The next step consists in computing the modular square roots of [math]N[/math] modulo the different primes in the factor base and computing the logarithms in base 2 of the different primes rounded to the nearest integer, storing these values in the arrays [math]\operatorname{tsqrt}[/math] and [math]\operatorname{tlog}[/math].

### Generating a polynomial

We have to compute the coefficients [math]a[/math] and [math]b[/math]. The first coefficient will be the square of a prime near [math]\sqrt{\sqrt{2N}/M}[/math] such that [math]N[/math] is a quadratic residue modulo this prime, and the second will be a modular square root of [math]a[/math]. All these operations will require multiple precision arithmetic.

For each prime [math]p[/math] in the factor base:

- Compute [math]\operatorname{ainv}[/math] as the inverse of [math]a[/math] modulo [math]p[/math].
- Compute the first solution as [math]\operatorname{soln1} = \operatorname{ainv} * (\operatorname{tsqrt}[p] - b) \bmod p[/math] and add [math]\operatorname{tlog}[p][/math] to all locations [math]\operatorname{soln1} + k*p[/math] of the sieve array.
- Compute the second solution as [math]\operatorname{soln2} = \operatorname{ainv} * (-\operatorname{tsqrt}[p] - b) \bmod p[/math] and add [math]\operatorname{tlog}[p][/math] to all locations [math]\operatorname{soln2} + k*p[/math] of the sieve array.

### Obtaining relations

Walk through the sieve array. For each element greater than about [math]\log(2x\sqrt {N})[/math] minus a small error term (because the logarithms are rounded and we do not sieve with powers of prime numbers), perform a trial division of [math]g_{a,b}(x)/a[/math] by the elements of the factor base. If all prime factors are less than or equal to F, save the relation. If there are enough relations (generally about 5% more relations than elements of the factor base), perform the linear algebra stage. Otherwise more relations must be collected. When the sieve stage is exhausted, change polynomial.

## Large prime variation

A variation that makes the sieving about twice as fast but requiring more memory, named PMPQS, consists in selecting a large prime bound (about 64 F). If after the trial division by the primes in the factor base the quotient is less than the large prime bound, we collect this *partial relation* in a different location than *full relations* (when the final quotient is 1). If the large prime found is the same than the large prime found in a previous partial relation, we can combine them:

Let [math]t^2\equiv u * p\,\pmod{N}[/math] and [math]r^2\equiv v * p\,\pmod{N}[/math] be two partial relations with the same large prime [math]p[/math]. Then a full relation will be:

[math]\large\left(\frac {t*r}{p}\right)^2\,\equiv u*v\,\pmod{N}[/math].

This full relation can be stored with the other full relations found.

At the beginning of PMPQS the number of full relations found in this way will be small, but it increases as the algorithm progresses, and then there will be more relations found from partial relations than without large primes.