How the self-initializing quadratic sieve works

A careful walkthrough of the algorithm inside SIQS.NET — from Fermat's four-century-old idea to the sparse linear algebra that finishes the job. ⏱ ~30-minute read

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.

1. 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 algorithm in that family that is genuinely fast, and for numbers up to roughly 100 digits it remained the world record holder 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.

2. 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

— each a fast operation — hands us a non-trivial factor. For a random such congruence with \(N\) a product of two primes, the probability that the split is non-trivial is at least \(1/2\). Collect a handful of independent congruences and success is all but guaranteed.

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. If we multiply a subset \(S\) of relations together, 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. So the search for a congruence of squares reduces to a search for a subset of relations whose exponent vectors sum to all-even values.

And that is linear algebra. Reduce every exponent modulo 2, so each relation becomes a 0/1 vector over the field \(\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 linearly 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)

3. 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 quadratic sieve generates its candidate values from a polynomial (we'll meet it shortly) whose values \(V(x)\) satisfy \(T(x)^2 \equiv A\cdot V(x) \pmod{N}\) for known \(T\) and \(A\). If a prime \(p\) divides such a \(V(x)\), then \(N\) is congruent to a square modulo \(p\) — in the standard vocabulary, \(N\) must be a quadratic residue mod \(p\). The Legendre symbol packages this up:

For a random \(N\), 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\) — about half survive. For each survivor, SIQS.NET also precomputes and stores the two modular square roots \(r\) with \(r^2 \equiv N \pmod p\) (via the Tonelli–Shanks algorithm), because those roots tell the siever exactly where \(p\) will strike. All of this lands in factor_base.txt, one line per prime: index, prime, both roots, and a scaled logarithm (more on that soon).

The multiplier trick

Here is a delightfully counter-intuitive optimization. Sometimes \(N\) itself has bad luck: the small primes 3, 5, 7, 11 — the ones that divide sieve values most often — happen to be non-residues, leaving the factor base with an anaemic small-prime section. The fix is to factor a different number: pick a small square-free multiplier \(k\) and run the whole algorithm on \(kN\) instead. 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 — expected log-contribution of each small prime to sieve values, minus a penalty for how much bigger \(k\) makes the numbers:

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 analysis \(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 (more relations per polynomial) but demands more relations in total — you need more rows than columns, and columns are primes. The optimal bound grows subexponentially with \(N\); SIQS.NET uses a smooth digit-calibrated formula anchored to the quantity \(L = \sqrt{\ln N \,\ln\ln N}\) (the same quantity that governs the whole algorithm's runtime — see §11), tuned so that a 20-digit input gets a bound near a thousand and a 90-digit input gets 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.

4. 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 (the overwhelming majority) 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.

Consider the basic quadratic-sieve polynomial \(Q(x) = (x + \lceil\sqrt{N}\rceil)^2 - N\). If \(p\) divides \(Q(x)\), then \((x + \lceil\sqrt{N}\rceil)^2 \equiv N \pmod p\) — so \(x + \lceil\sqrt{N}\rceil\) must be congruent to one of the two precomputed square roots \(\pm r\) of \(N\) mod \(p\). In other words, the solutions form two arithmetic progressions with common difference \(p\). Once you find the first hit of each, every later hit is a constant stride away. No divisions, no remainders: just “jump ahead by \(p\), mark; jump ahead by \(p\), mark.”

The marking is done with logarithms. Allocate an array over the sieve interval, one byte per candidate \(x\), initialized to zero. For each factor-base prime \(p\), walk its two progressions and add \(\log p\) (scaled to fit a byte) at each hit. A value that is fully smooth will have received log-contributions summing to approximately \(\log|V(x)|\) — so after all primes have been walked, a single linear scan flags every position whose accumulated sum exceeds

where \(\varepsilon\) is a configured error margin that absorbs rounding, skipped tiny primes, and the possibility of one moderate leftover factor (which §6 will turn into a feature). Only the flagged survivors — a tiny fraction — are handed to exact trial division, which reconstructs the full factorization, confirms smoothness, and emits the relation. The expensive certainty of division is spent only where the cheap evidence of addition says it's worth it.

5. Self-initialization: many polynomials, cheaply

The single polynomial \(Q(x)\) has a weakness: its values grow. Near \(x = 0\), \(Q(x) \approx 2x\sqrt{N}\) is comparatively small and pleasantly likely to be smooth; a million steps out, the values 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}\) (making \(C\) an exact integer). A one-line calculation then gives the congruence that makes relations possible:

So 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 polynomial's values are minimized in the worst case when

which balances \(|V(\pm M)|\) against \(|V(0)|\) so the whole interval sieves well. 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 depending on the size of \(N\)) — because that structure is exactly what makes many \(B\) values available. SIQS.NET builds a window of factor-base primes near the ideal size \(A^{1/s}\), scores each \(s\)-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\), being a factor-base prime, already has its two stored square roots of \(\widetilde N\). Build one basis term per prime:

where \(r_i\) is the stored root. 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 choice of \(A\) therefore yields \(2^s\) polynomials — of which half are kept: \(\varepsilon\) and \(-\varepsilon\) produce mirror-image polynomials whose relations are arithmetically dependent (over the symmetric interval they generate guaranteed-trivial square roots), so SIQS.NET enumerates exactly one representative of each pair, \(2^{s-1}\) polynomials per \(A\), using the first half of a binary-reflected Gray code.

The Gray code is the punchline, not a detail. Successive sign 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 for every prime \(p\) shift by a precomputed per-prime constant:

Initializing a polynomial from scratch costs a modular inverse per factor-base prime — expensive when repeated tens of thousands of times. 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, and the siever can afford to keep every interval short and every value small.

6. 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 sieve threshold's error margin was chosen deliberately so they surface. A survivor whose cofactor is a single prime \(\ell\) with \(B < \ell \le \ell_{\max}\) (typically \(\ell_{\max}\) is a few hundred times \(B\)) is recorded as a partial relation: a relation with one extra prime attached.

One partial is useless — its \(\ell\) appears to an odd power. But two partials sharing the same \(\ell\) can be multiplied: \(\ell\) then appears squared, contributes even parity, vanishes from the matrix, and the pair behaves exactly like one full relation. By the birthday paradox, collisions among large primes arrive far more often than intuition suggests, and in practice partials contribute a large share of the final matrix at nearly no extra sieving cost.

SIQS.NET pushes one step further for large targets: two-large-prime (2LP) relations, where the leftover cofactor is a product of two primes \(\ell_1 \ell_2\), split by the SQUFOF or Pollard-rho algorithm. Combining these is no longer simple pairing — it becomes a graph problem. Build a graph whose vertices are large primes (plus a special vertex “1” for single-prime partials): every partial is an edge between its large primes. Any cycle in this graph is a set of partials in which every large prime is touched an even number of times — multiply the partials along the cycle and all the large primes cancel. The filtering phase hunts these cycles with a union-find forest. For 90-digit-class targets, 2LP mode meaningfully extends the practical reach of the implementation.

7. Filtering: preparing the matrix

Sieving ends with hundreds of thousands of raw relations in relations_*.txt and partials_*.txt. 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 provenance and the cancelled primes — the square-root phase will need them.
• Prune singletons. A prime that appears 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. This typically shrinks the matrix dramatically.
• Trim the surplus. Keep only modestly more rows than columns (each excess row guarantees one more dependency; a few dozen 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, since \(V(x) < 0\) happens on half the interval and \(-1\) must cancel like any prime). Alongside it, relations_filtered.txt preserves the full arithmetic — exponents, \(t\)-values, large-prime cofactors — that the parity-only matrix deliberately forgets.

8. 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 it suffers fill-in — the matrix starts a few dozen nonzeros per row and elimination steadily densifies it until it no longer fits in memory, with \(O(n^3)\) bit-operations on top.

SIQS.NET instead uses Block Lanczos (Montgomery, 1995), the standard workhorse for exactly this job. Two ideas power it:

• The matrix is never modified. Lanczos-type methods touch \(M\) only through matrix–vector products \(v \mapsto Mv\) and \(v \mapsto M^{\mathsf T}v\), which cost 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 vector coordinates, and one hardware XOR processes all of them. The “block” variant makes the word width work as a 64-fold parallel speedup — and has the pleasant side effect of emitting up to 64 nullspace vectors in a single run.

The iteration builds an orthogonal-like sequence of block subspaces and converges after roughly \(n/64\) matrix passes, giving total cost \(O(n \cdot w)\) XOR-word operations for a matrix with \(w\) nonzeros per row — a completely different regime from elimination. Each recovered nullspace vector selects a subset of filtered relations whose parities all cancel: dependencies.txt lists them, typically several dozen per run. Each one is an independent lottery ticket for the final phase.

9. 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\) is the (verified-even) total exponent of factor-base prime \(p_j\), and each combined partial contributes one copy of its cancelled large prime \(\ell\) — the matrix forgot those primes, but arithmetic doesn't. By construction \(X^2 \equiv Y^2 \pmod{\widetilde N}\), hence also mod \(N\), and

are computed against the original \(N\) — the multiplier was scaffolding, and it is dismantled here. Each dependency independently succeeds with probability at least \(1/2\) (it fails only when \(X \equiv \pm Y\), the “trivial” outcome); with sixty-odd dependencies banked, the chance that all of them 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 shape of 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 — bigger factor base versus more required relations, longer intervals versus larger values, more polynomials versus more setup — 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 of \(N\) 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.)

In practice, asymptotics only set the shape; constants are found by measurement. SIQS.NET derives every default from the digit count \(D\) of the target — each calibrated by trial sieving on real targets rather than theory alone:

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, while primes larger than a block — which hit each block at most once per root — are routed through per-block bucket hit lists, and a middle band is rediscovered by resieving, so the per-candidate divisibility gate runs only over genuinely small primes.
• Parallel and distributed sieving. Polynomial families are independent, so they parallelize across cores trivially — and across machines: the Overlord coordinator leases disjoint family ranges to volunteer clients, which rebuild the factor base themselves and verify the job before contributing (see the distributed overview).
• 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 — including re-validating that any supplied parameters match the stored ones.

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.