Quadratic Sieve

Python / Number Theory / Cryptography

An implementation of the quadratic sieve, the sub-exponential integer factoring algorithm behind practical attacks on small RSA keys. Written in Python as a notebook following Carl Pomerance's 2008 paper directly, then with vectorized performance optimizations.

Implementation Details

The algorithm selects a smoothness bound B via the heuristic exp(0.5*sqrt(ln(n)*ln(ln(n)))) and builds a factor base from the primes p below B where the Legendre symbol (n/p) equals 1, the only primes that can divide \(Q(x) = x^2 - n\). For each such prime we solve the modular square root congruence \(x^2 \equiv n \pmod p\) with the Tonelli-Shanks algorithm; the two roots define the arithmetic progressions along which the sieve accumulates.

Algorithm Pseudocode — Fig 1. Formal definition of the Quadratic Sieve factorization implemented in Python. \(S[i]\) denotes the accumulated log-sum at index \(i\) of the sieve array, used to identify smooth numbers. The final step computes the GCD using the constructed congruence \(X^2 \equiv Y^2 \pmod n\), where \(sqrt(\prod x_i^2)\) is the product of the collected roots (\(X\)) and \(sqrt(\prod Q(x_i))\) is the square root of the product of the B-smooth numbers (\(Y\)).

Sieving proceeds in blocks over intervals starting at \(\lceil\sqrt{n}\rceil\). Each block is a NumPy float32 array, and for every factor-base prime we add \(\log_2 p\) at the root offsets using vectorized slice arithmetic (sieve[offset::p] += log(p)) instead of per-index loops. Indices whose accumulated log-sum exceeds a threshold are flagged as candidate B-smooth values; the threshold sits at 75% of the expected magnitude \(\log_2 Q(x)\) to tolerate rounding error and small unsieved factors.

Log-Sum Visualization — Fig 2. Visualization of the sieving interval. The red line represents the target threshold derived from the polynomial. Indices where the accumulated log sum exceeds this threshold (blue bar) are flagged as candidate smooth numbers.

Candidates are confirmed via trial division over the factor base, and their exponent vectors are collected until the relation count exceeds the factor base size by a safety margin. We reduce the matrix of exponent parities over GF(2) with Gaussian elimination and extract a basis for its null space; each kernel vector selects a subset of relations whose product is a perfect square. This yields a congruence of squares \(X^2 \equiv Y^2 \pmod n\), and gcd(X - Y, n) reveals a non-trivial factor, retrying with the next kernel vector whenever the split is trivial.

GF(2) Matrix Reduction — Fig 3. The linear algebra part, solving for the null space of the exponent matrix (mod 2) to find a linear dependency, identifying the perfect square required to factor the composite number.