Python / Number Theory / Cryptography
The algorithm begins by determining a smoothness bound B using the heuristic exp(sqrt(ln(n)*ln(ln(n)))/2). We construct a factor base of prime numbers p ≤ B where the Legendre symbol (n/p) equals 1. To efficiently identify B-smooth numbers, we solve the congruence x2 ≡ n (mod p) for each prime in the base using the Tonelli-Shanks algorithm. These roots define the starting offsets for our sieve.
We initialize a sieve array over an interval and add the logarithm of p to indices corresponding to the arithmetic progressions of the roots. Indices where the accumulated log-sum exceeds a threshold (derived from the log of the target polynomial) are flagged as candidate B-smooth numbers. We verify these candidates via trial division and collect their exponent vectors. Once the number of relations exceeds the size of the factor base, we form a matrix of exponent parities. We perform Gaussian elimination over the field GF(2) to find the null space. Each vector in the kernel corresponds to a subset of relations whose product is a perfect square. This yields the congruence X2 ≡ Y2 (mod n), allowing us to compute gcd(X - Y, n) to reveal the non-trivial prime factors.