How the self-initializing quadratic sieve works

This page assumes comfort with modular arithmetic and a tolerance for prime numbers, but nothing beyond undergraduate mathematics. Every quantity named here — factor_base.txt, relation files, the matrix, the dependencies — is a real text artifact that SIQS.NET writes to disk during a run, so you can follow along inside an actual factorization.

§01 The problem

Take two large primes \(p\) and \(q\), multiply them, and hand someone the product \(N = p \cdot q\). Recovering \(p\) and \(q\) from \(N\) is the integer factorization problem, and for well-chosen primes it is spectacularly harder than the multiplication was. That asymmetry is not a curiosity: the security of RSA rests on it. Factoring the RSA-129 challenge number in 1994 — using precisely the algorithm described on this page — took eight months and around 600 volunteers' machines.

Trial division tests candidate divisors one at a time and needs on the order of \(\sqrt{N}\) steps; for a 90-digit \(N\) that is roughly \(10^{45}\) operations — not happening. Every serious factoring algorithm therefore abandons the search for divisors and instead manufactures a special kind of coincidence.

The quadratic sieve is the simplest genuinely fast algorithm in that family.

For numbers up to roughly 100 digits it held the world record until the general number field sieve overtook it. SIQS — the self-initializing quadratic sieve — is its refined, modern form, and the one SIQS.NET implements.

§02 Fermat's idea, industrialized

In the 1640s Fermat observed that finding a difference of squares factors a number directly. If

then the two parenthesized terms are factors. Fermat's method searches for such \(x\) by walking upward from \(\lceil\sqrt{N}\rceil\) and testing whether \(x^2 - N\) is a perfect square. It works beautifully when \(p\) and \(q\) are close together and abysmally otherwise.

The industrial upgrade, due to Kraitchik in the 1920s, is to relax equality to congruence. We no longer need \(x^2 - y^2\) to equal \(N\) — it is enough that

Then \(N\) divides \(X^2 - Y^2 = (X-Y)(X+Y)\) but divides neither bracket by itself, so the prime factors of \(N\) must split between the two brackets. Computing \(\gcd(X - Y, N)\) and \(\gcd(X + Y, N)\) — each fast — hands us a non-trivial factor. For a random such congruence with \(N\) a product of two primes, the split is non-trivial with probability at least \(1/2\).

But how do you find \(X\) and \(Y\)? Here is the key manoeuvre, made systematic by Dixon: don't look for one perfect coincidence — multiply many small ones together. Suppose we collect many values \(x_i\) for which we completely know the prime factorization of \(x_i^2 \bmod N\):

Each such factored congruence is called a relation. Multiply a subset \(S\) of relations together and the left side is automatically a square — it is \(\bigl(\prod_{i \in S} x_i\bigr)^2\). The right side is a square precisely when every prime's total exponent \(\sum_{i\in S} e_{ij}\) is even.

And that is linear algebra. Reduce every exponent modulo 2, so each relation becomes a 0/1 vector over \(\mathbb{F}_2 = \mathrm{GF}(2)\), one coordinate per prime. A subset with all-even exponent sums is exactly a set of vectors that XOR to zero — a linear dependency. Collect more relations than there are primes in play, and the vectors must be dependent: dependencies are guaranteed by pigeonhole, and each one yields a candidate congruence of squares.

The quadratic sieve is, at heart, a very fast relation factory bolted onto this framework. Three questions remain, and they structure the rest of this page:

• Which primes do we allow on the right-hand side? (§3, the factor base)
• How do we find values whose squares factor over those primes, fast? (§4–5, the sieve and SIQS)
• How do we find the dependency and finish? (§7–9, filtering, linear algebra, square roots)

§03 Smooth numbers and the factor base

A number is called \(B\)-smooth if all of its prime factors are at most \(B\). Smooth numbers are the currency of the quadratic sieve: a relation is precisely a value whose square, reduced mod \(N\), is smooth over our chosen prime list. That list is the factor base.

Only half the primes can play

The sieve generates candidate values from a polynomial whose values \(V(x)\) satisfy \(T(x)^2 \equiv A\cdot V(x) \pmod{N}\). If a prime \(p\) divides such a \(V(x)\), then \(N\) is a square modulo \(p\) — a quadratic residue. The Legendre symbol packages this up:

Each odd prime has a 50/50 chance, so the factor base is built by generating all primes up to a bound \(B\) with the sieve of Eratosthenes and keeping only those with \(\left(\frac{N}{p}\right) = 1\). For each survivor, SIQS.NET precomputes the two modular square roots \(r\) with \(r^2 \equiv N \pmod p\) (via Tonelli–Shanks), because those roots tell the siever exactly where \(p\) will strike. All of this lands in factor_base.txt: index, prime, both roots, and a scaled logarithm.

The multiplier trick

A delightfully counter-intuitive optimization. Sometimes \(N\) has bad luck: the small primes 3, 5, 7, 11 — the ones that divide sieve values most often — happen to be non-residues, leaving an anaemic small-prime section. The fix is to factor a different number: pick a small square-free multiplier \(k\) and run the algorithm on \(kN\). A factor of \(kN\) that isn't a factor of \(k\) gives a factor of \(N\) immediately, and a good \(k\) can make many more small primes usable.

SIQS.NET scores each candidate \(k \in \{1, 2, 3, 5, 7, \dots, 47\}\) with the classical Knuth–Schroeppel function:

where \(g(p, kN) = \frac{2\log p}{p-1}\) when \(kN\) is a residue mod \(p\), \(\frac{\log p}{p}\) when \(p \mid k\), and \(0\) otherwise, with a separate case \(g_2\) for the prime 2 depending on \(kN \bmod 8\). The winner becomes the multiplier, the pipeline works with \(\widetilde{N} = kN\) throughout, and only the very last GCD step returns to the original \(N\).

How big should the factor base be?

A larger bound \(B\) makes smoothness easier but demands more relations in total. The optimal bound grows subexponentially with \(N\); SIQS.NET uses a digit-calibrated formula anchored to \(L = \sqrt{\ln N \,\ln\ln N}\) (see §11), tuned so a 20-digit input gets a bound near a thousand and a 90-digit input a bound in the millions. For the 90-digit class the factor base holds roughly 90,000 primes; the tuned 110-digit profile pushes the bound to 40 million.

§04 The sieve: why it's fast

We need many values \(V(x)\) that are smooth over the factor base. The naive approach — trial-divide every candidate by every factor-base prime — spends almost all its time proving that non-smooth numbers are non-smooth. The sieve inverts the logic: instead of asking “what divides this value?”, ask “which values does this prime divide?” — and answer it for free with arithmetic progressions.

The marking is done with logarithms. Allocate one byte per candidate. For each factor-base prime \(p\), walk its two progressions and add \(\log p\) (scaled to a byte) at each hit. A fully smooth value receives log-contributions summing to about \(\log|V(x)|\), so a single scan flags every position whose sum exceeds

where \(\varepsilon\) is a configured error margin absorbing rounding, skipped tiny primes, and the possibility of one moderate leftover factor (which §6 turns into a feature).

§05 Self-initialization: many polynomials, cheaply

The single polynomial \(Q(x) = (x + \lceil\sqrt{N}\rceil)^2 - N\) has a weakness: its values grow. Near \(x = 0\) they are small and pleasantly likely to be smooth; a million steps out they are a million times larger and the smoothness probability has collapsed. Sieving one polynomial forever means grinding an ever-drier quarry.

The fix (Montgomery's, refined by Alford–Pomerance into “self-initialization”) is to use many polynomials, each sieved only over a short interval \([-M, M]\) where its values stay small. SIQS.NET uses the general quadratic family

with coefficients engineered so that \(B^2 \equiv \widetilde{N} \pmod{A}\). A one-line calculation then gives the congruence that makes relations possible:

Each smooth value of \(A \cdot V(x)\) yields a relation with the known value \(T(x)\) on the square side — and since we choose \(A\) ourselves as a product of factor-base primes, its factorization contributes known exponents for free.

Choosing A

Two constraints shape \(A\). First, magnitude: over \([-M, M]\) the values are minimized in the worst case when \(A \approx \frac{\sqrt{2\widetilde{N}}}{M}\), which balances \(|V(\pm M)|\) against \(|V(0)|\). Second, structure: \(A\) must be a product of \(s\) distinct odd factor-base primes \(q_1 \cdots q_s\) (with \(s\) between 3 and 10) — because that structure is exactly what makes many \(B\) values available. SIQS.NET builds a window of primes near \(A^{1/s}\), scores each subset by \(\bigl|\log A - \log A_{\text{target}}\bigr|\), and consumes candidates in deterministic best-first order.

One A, a family of Bs — the self-initializing part

Fix \(A = q_1 q_2 \cdots q_s\). We need \(B\) with \(B^2 \equiv \widetilde{N} \pmod A\). By the Chinese Remainder Theorem it suffices to satisfy the congruence modulo each \(q_i\) separately — and each \(q_i\) already has its two stored square roots. Build one basis term per prime:

then every sign pattern \(\varepsilon \in \{\pm 1\}^s\) gives a valid coefficient

because squaring mod \(q_i\) kills every term except \((\pm B_i)^2 \equiv \widetilde N\). One \(A\) therefore yields \(2^s\) polynomials — of which half are kept, since \(\varepsilon\) and \(-\varepsilon\) are mirror images. SIQS.NET enumerates \(2^{s-1}\) per \(A\) using the first half of a binary-reflected Gray code.

The Gray code is the punchline. Successive patterns differ in one bit, so stepping from one polynomial to the next changes \(B\) by \(\pm 2B_v\) for a single index \(v\) — and the sieve's starting positions shift by a precomputed per-prime constant:

Initializing a polynomial from scratch costs a modular inverse per prime; switching within a family costs one addition per prime. That is the entire meaning of self-initializing: the family amortizes its setup so thoroughly that polynomial switching becomes almost free.

§06 Large primes: profiting from near-misses

Most threshold survivors fail exact trial division not by much — they factor over the base except for one leftover prime a bit beyond the bound. Discarding these is wasteful, and the error margin was chosen deliberately so they surface. A survivor whose cofactor is a single prime \(\ell\) with \(B < \ell \le \ell_{\max}\) is recorded as a partial relation: a relation with one extra prime attached.

One partial is useless; two partials sharing the same \(\ell\) can be multiplied so \(\ell\) appears squared, contributes even parity, and vanishes from the matrix. In practice partials contribute a large share of the final matrix at nearly no extra sieving cost.

For a two-large-prime (2LP) relation the leftover cofactor is a product \(\ell_1 \ell_2\), split by SQUFOF or Pollard-rho. Combining these becomes a graph problem: vertices are large primes (plus a special vertex “1” for single-prime partials), and every partial is an edge. Any cycle is a set of partials in which every large prime is touched an even number of times — multiply along the cycle and the large primes cancel.

§07 Filtering: preparing the matrix

Sieving ends with hundreds of thousands of raw relations. Filtering turns them into the smallest clean matrix that still guarantees dependencies:

• Deduplicate. Parallel workers and overlapping polynomials can find the same relation twice; duplicates create fake dependencies that factor nothing.
• Combine partials. Find large-prime cycles (§6) and merge each into a synthetic combined_partial relation, recording the cancelled primes — the square-root phase will need them.
• Prune singletons. A prime appearing in exactly one relation can never cancel; that relation is dead weight. Remove it — which may create new singletons — and iterate to a fixed point.
• Trim the surplus. Keep only modestly more rows than columns (a few dozen excess rows suffice), preferring to discard the heaviest rows to keep the matrix sparse.

The output is a sparse 0/1 matrix — filtered_matrix.txt lists, for each relation, only the indexes of its odd-exponent columns (column 0 is a virtual sign column tracking negative values). Alongside it, relations_filtered.txt preserves the full arithmetic that the parity-only matrix deliberately forgets.

§08 Linear algebra over GF(2)

Now find a non-empty set of rows that XOR to zero. For a 90-digit factorization the matrix has on the order of \(10^5\) rows and columns; Gaussian elimination would work in principle but suffers fill-in — it densifies the matrix until it no longer fits in memory, with \(O(n^3)\) bit-operations on top.

SIQS.NET instead uses Block Lanczos (Montgomery, 1995). Two ideas power it:

• The matrix is never modified. Lanczos-type methods touch \(M\) only through products \(v \mapsto Mv\) and \(v \mapsto M^{\mathsf T}v\), one XOR per stored nonzero. Sparsity is preserved forever; memory stays flat.
• Work on 64 vectors at once. Over \(\mathrm{GF}(2)\), a machine word holds 64 independent coordinates, and one hardware XOR processes all of them — a 64-fold speedup that also emits up to 64 nullspace vectors per run.

The iteration converges after roughly \(n/64\) matrix passes, for a total cost \(O(n \cdot w)\) XOR-word operations on a matrix with \(w\) nonzeros per row. Each recovered nullspace vector selects a subset of filtered relations whose parities all cancel: dependencies.txt lists them, typically several dozen per run. Each is an independent lottery ticket for the final phase.

§09 The square root: extracting the factors

A dependency asserts: over this set \(S\) of relations, every prime's total exponent is even. Cash it in:

where \(t_i\) is each relation's stored square-side value, \(E_j\) the verified-even total exponent of prime \(p_j\), and each combined partial contributes one copy of its cancelled large prime \(\ell\). By construction \(X^2 \equiv Y^2 \pmod{\widetilde N}\), hence also mod \(N\).

The two GCDs

are computed against the original \(N\). Each dependency independently succeeds with probability at least \(1/2\); with sixty-odd dependencies banked, the chance that all fail is smaller than \(2^{-60}\). In practice the first or second attempt splits \(N\), the factors are verified by multiplication, and factors.txt is written.

§10 A worked run, end to end

SIQS.NET's test suite pins a complete generated artifact set for the 29-digit target

with every cross-phase file checked against the contracts described above. The run:

The final GCDs deliver:

and indeed \(123456789012419 \times 678901234567937 = N\). Note the economics even at this toy scale: 2555 near-misses were collected almost free and condensed into 760 extra matrix rows — the large-prime variation more than doubled the effective yield of the sieve.

§11 Complexity and parameter choices

Every tension in the algorithm resolves into one quantity. Under standard heuristics on the density of smooth numbers, the optimized cost of the quadratic sieve is

with the optimal factor-base bound near \(L(N)^{1/2}\). This is subexponential: dramatically better than trial division's \(\exp\bigl(\tfrac12\ln N\bigr)\), yet still growing fast enough that each additional 10 digits roughly quadruples the work at these sizes. (The number field sieve improves the exponent to \((\ln N)^{1/3}(\ln\ln N)^{2/3}\), which is why it wins beyond ~100 digits — at the price of vastly more machinery.)

SIQS.NET derives every default from the digit count \(D\) of the target, each calibrated by trial sieving on real targets:

Every one of these is overridable from the qs command line, which is precisely how the defaults were found: run --trial-sieve-percent samples across a parameter sweep, compare raw relation throughput, commit the winner.

§12 Where the engineering lives

The mathematics above fits on a napkin; making it fast is where the implementation earns its keep. A few highlights of what you'll find in the source:

• Cache-blocked, SIMD-filled sieving. The interval is processed in cache-sized blocks (256 KiB–1 MiB); small primes use AVX2/scalar direct fills, larger primes are routed through per-block bucket hit lists, and a middle band is rediscovered by resieving.
• Parallel and distributed sieving. Polynomial families are independent, so they parallelize across cores and across machines: the Overlord leases disjoint family ranges to clients that rebuild the factor base themselves and verify the job before contributing.
• Deterministic replay. Every phase is deterministic given its inputs; with --parallelism 1, whole runs are byte-for-byte reproducible — the property the test suite's pinned end-to-end fixture is built on.
• Text artifacts everywhere. Each phase's output is a documented UTF-8 file with a versioned header. Debugging a factorization is reading files, not attaching debuggers.
• Resumability. A run's workspace records enough state that qs --resume <run-dir> can pick up an interrupted job.

The best way to make all of this concrete is to watch it happen: factor something, then open the run directory and read the artifacts alongside this page.