16.4.2 The quadratic sieve algorithm

We now describe a practical improvement to Algorithm SEF. This algorithm,

known as the quadratic sieve, is faster in practice than Algorithm SEF;

however, its analysis is somewhat heuristic.

First, let us return to the simpli¬ed version of Algorithm SEF, where

we collect relations of the form (16.6). Furthermore, instead of choosing

the values ±i at random, we will choose them in a special way, as we now

describe. Let

√

n :=

˜ n,

and de¬ne the polynomial

F := (X + n)2 ’ n ∈ Z[X].

˜

In addition to the usual “smoothness parameter” y, we need a “sieving

parameter” z, whose choice will be discussed below. We shall assume that

both y and z are of the form exp[(log n)1/2+o(1) ], and our ultimate choices

of y and z will indeed satisfy this assumption.

For all s = 1, 2, . . . , z , we shall determine which values of s are “good,”

in the sense that the corresponding value F (s) is y-smooth. For each good

s, since we have F (s) ≡ (s + n)2 (mod n), we obtain one relation of the form

˜

(16.6), with ±i := [s + n]n . If we ¬nd at least k + 1 good values of s, then

˜

we can apply Gaussian elimination as usual to ¬nd a square root β of 1 in

Zn . Hopefully, we will have β = ±1, allowing us to split n.

Observe that for 1 ¤ s ¤ z, we have

1 ¤ F (s) ¤ z 2 + 2zn1/2 ¤ n1/2+o(1) .

Now, although the values F (s) are not at all random, we might expect

heuristically that the number of good s up to z is roughly equal to σ z, ˆ

1/2 ] is

where σ is the probability that a random integer in the interval [1, n

ˆ

y-smooth, and by Theorem 16.7, we have

σ = exp[(’1/4 + o(1))(log n/ log y) log log n].

ˆ

354 Subexponential-time discrete logarithms and factoring

If our heuristics are valid, this already gives us an improvement over Al-

gorithm SEF, since now we are looking for y-smooth numbers near n1/2 ,

which are much more common than y-smooth numbers near n. But there

is another improvement possible; namely, instead of testing each individual

number F (s) for smoothness using trial division, we can test them all at

once using the following “sieving procedure”:

Create a vector v[1 . . . z ], and initialize v[s] to F (s), for 1 ¤

s ¤ z. For each prime p up to y, do the following:

1. Compute the roots of the polynomial F modulo p.

This can be done quite e¬ciently, as follows. For p = 2,

F has exactly one root modulo p, which is determined

by the parity of n. For p > 2, we may use the fa-

˜

miliar quadratic formula together with an algorithm for

computing square roots modulo p, as discussed in Exer-

cise 13.3. A quick calculation shows that the discrimi-

nant of F is n, and thus, F has a root modulo p if and

only if n is a quadratic residue modulo p, in which case

it will have two roots (under our usual assumptions, we

cannot have p | n).

2. Assume that the distinct roots of F modulo p lying in

the interval [1, p] are ri , for i = 1, . . . , vp .

Note that vp = 1 for p = 2 and vp ∈ {0, 2} for p > 2.

Also note that F (s) ≡ 0 (mod p) if and only if s ≡

ri (mod p) for some i = 1, . . . , vp .

For i = 1, . . . , vp , do the following:

s ← ri

while s ¤ z do

repeat v[s] ← v[s]/p until p v[s]

s←s+p

At the end of this sieving procedure, the good values of s may be identi¬ed

as precisely those such that v[s] = 1. The running time of this sieving

procedure is at most len(n)O(1) times

z 1

= O(z log log y) = z 1+o(1) .

=z

p p

p¤y p¤y

Here, we have made use of Theorem 5.10, although this is not really nec-

essary ” for our purposes, the bound p¤y (1/p) = O(log y) would su¬ce.

16.4 Practical improvements 355

Note that this sieving procedure is a factor of k 1+o(1) faster than the method

for ¬nding smooth numbers based on trial division. With just a little extra

book-keeping, we can not only identify the good values of s, but we can also

compute the factorization of F (s) into primes.

Now, let us put together all the pieces. We have to choose z just large

enough so as to ¬nd at least k + 1 good values of s up to z. So we should

choose z so that z ≈ k/ˆ ” in practice, we could choose an initial estimate

σ

for z, and if this choice of z does not yield enough relations, we could keep

doubling z until we do get enough relations. Assuming that z ≈ k/ˆ , the

σ

cost of sieving is (k/ˆ )1+o(1) , or

σ

exp[(1 + o(1))(1/4)(log n/ log y) log log n + log y].

The cost of Gaussian elimination is still O(k 3 ), or

exp[(3 + o(1)) log y].

Thus, if T is the running time of the entire algorithm, we have

T ¤ exp[(1 + o(1)) max{(1/4)(log n/ log y) log log n + log y, 3 log y}].

Let µ := log y, A := (1/4) log n log log n, S1 := A/µ + µ and S2 := 3µ, and

let us ¬nd the value of µ that minimizes max{S1 , S2 }. Using a little calculus,

one ¬nds that S1 is minimized at µ = A1/2 . For this value of µ, we have

S1 = 2A1/2 and S2 = 3A1/2 > S1 , and so this choice of µ is a bit larger

than optimal. For µ < A1/2 , S1 is decreasing (as a function of µ), while S2

is always increasing. It follows that the optimal value of µ is obtained by

setting

A/µ + µ = 3µ

and solving for µ. This yields µ = (A/2)1/2 . So setting

√

y = exp[(1/(2 2))(log n log log n)1/2 ],

we have

√

T ¤ exp[(3/(2 2) + o(1))(log n log log n)1/2 ].

Thus, we have reduced the constant in the exponent from 2, for Algorithm

√

SEF (using the more accurate smoothness density estimates), to 3/(2 2) ≈

1.061.

We mention one ¬nal improvement. The matrix to which we apply Gaus-

sian elimination in stage 2 is “sparse”; indeed, since any integer less than

n has O(log n) prime factors, the total number of non-zero entries in the

356 Subexponential-time discrete logarithms and factoring

matrix is k 1+o(1) . In this case, there are special algorithms for working with

such sparse matrices, which allow us to perform stage 2 of the factoring

algorithm in time k 2+o(1) , or

exp[(2 + o(1)) log y].

This gives us

T ¤ exp[(1 + o(1)) max{(1/4)(log n/ log y) log log n + log y, 2 log y}],

and setting

y = exp[(1/2)(log n log log n)1/2 ]

yields

T ¤ exp[(1 + o(1))(log n log log n)1/2 ].

Thus, this improvement reduces the constant in the exponent from

√

3/(2 2) ≈ 1.061 to 1. Moreover, the special algorithms designed to work

with sparse matrices typically use much less space than ordinary Gaussian

elimination ” even if the input to Gaussian elimination is sparse, the inter-

mediate matrices will not be. We shall discuss in detail later, in §19.4, one

such algorithm for solving sparse systems of linear equations.

The quadratic sieve may fail to factor n, for one of two reasons: ¬rst, it

may fail to ¬nd k + 1 relations; second, it may ¬nd these relations, but in

stage 2, it only ¬nds a trivial square root of 1. There is no rigorous theory

to say why the algorithm should not fail for one of these two reasons, but

experience shows that the algorithm does indeed work as expected.

16.5 Notes