Primality testing and factoring

Gauss wrote in the Disquisitiones that the problem of distinguishing prime numbers from composite numbers, and of resolving the latter into their prime factors, is one of the most important and useful in arithmetic. Two centuries later the two halves of his problem have turned out to be of very different difficulty. Deciding whether a 600-digit number is prime takes a laptop a fraction of a second. Finding the factors of a 600-digit product of two primes is, as far as anyone knows, beyond all the computers on Earth. That gap is exactly what public-key cryptography exploits. This lesson explains both halves: fast probabilistic primality tests built from Fermat's theorem, and the cleverest elementary factoring methods.

Trial division and its limits

The simplest test divides \(n\) by every prime up to \(\sqrt n\). For a 20-digit number that is about \(4 \times 10^8\) divisions, quite feasible. For a 100-digit number it is about \(10^{48}\), and even at a billion divisions per second this would take far longer than the age of the universe. The cost grows exponentially in the number of digits. A useful algorithm must run in time polynomial in the number of digits, \(\log n\), not in \(n\) or \(\sqrt n\).

The Fermat test and its failure

Fermat's theorem gives a fast test: pick a base \(a\) and compute \(a^{n-1} \bmod n\) by square-and-multiply, which takes about \(2\log_2 n\) multiplications. If the result is not 1, then \(n\) is certainly composite, and \(a\) is called a Fermat witness. If the result is 1, \(n\) is probably prime, but not certainly: a composite \(n\) with \(a^{n-1} \equiv 1 \pmod n\) is a Fermat pseudoprime to base \(a\). The smallest pseudoprime to base 2 is \(341 = 11 \cdot 31\).

Pseudoprimes to base 2 are rare. There are 78,498 primes below a million but only 245 base-2 pseudoprimes. One might hope that testing several bases removes them all. Carmichael numbers spoil that hope: they are composites \(n\) with \(a^{n-1} \equiv 1 \pmod n\) for every \(a\) coprime to \(n\).

Theorem (Korselt's criterion, 1899). A composite \(n\) is a Carmichael number if and only if \(n\) is squarefree and \(p - 1 \mid n - 1\) for every prime \(p\) dividing \(n\).

Proof of sufficiency. If \(p \mid n\) and \(p - 1 \mid n - 1\), then for \(a\) coprime to \(n\), Fermat gives \(a^{n-1} = (a^{p-1})^{(n-1)/(p-1)} \equiv 1 \pmod p\). Since \(n\) is a product of distinct primes, each dividing \(a^{n-1} - 1\), so does \(n\). (Necessity uses primitive roots modulo \(p\) and modulo \(p^2\).) Example: \(561 = 3 \cdot 11 \cdot 17\), and 2, 10 and 16 all divide 560. Other examples are \(1105 = 5 \cdot 13 \cdot 17\), \(1729 = 7 \cdot 13 \cdot 19\) and \(2465 = 5 \cdot 17 \cdot 29\). Every Carmichael number has at least three prime factors. In 1994 Alford, Granville and Pomerance proved that there are infinitely many.

The Miller–Rabin test

The repair uses one more fact about primes: modulo a prime, the only square roots of 1 are \(\pm 1\), by Lagrange's theorem. For odd \(n\), write \(n - 1 = 2^sd\) with \(d\) odd. If \(n\) is prime and \(p \nmid a\), consider the sequence \[a^d,\ a^{2d},\ a^{4d},\ \ldots,\ a^{2^sd} = a^{n-1} \pmod n.\] Each term is the square of the one before, and the last is 1 by Fermat. Going backwards from the last term, each earlier term is a square root of the next, so as long as we are at 1 the previous term is \(\pm 1\). Hence either \(a^d \equiv 1\), or some term \(a^{2^rd} \equiv -1\) with \(0 \le r < s\).

Definition. An odd composite \(n\) is a strong pseudoprime to base \(a\) if \(a^d \equiv 1\) or \(a^{2^rd} \equiv -1 \pmod n\) for some \(0 \le r < s\). Otherwise \(a\) is a strong witness that \(n\) is composite.

Example: \(n = 561\), \(n - 1 = 560 = 2^4 \cdot 35\). With \(a = 2\): \(2^{35} \equiv 263\), \(2^{70} \equiv 166\), \(2^{140} \equiv 67\), \(2^{280} \equiv 1 \pmod{561}\). The sequence reaches 1 from 67, which is a square root of 1 other than \(\pm 1\). That is impossible modulo a prime, so 561 is composite. In fact it reveals a factor: \(67^2 \equiv 1\) means \(561 \mid (67 - 1)(67 + 1)\), and \(\gcd(66, 561) = 33\).

Theorem (Rabin, Monier, 1980). If \(n\) is an odd composite greater than 9, at most one quarter of the bases \(a\) in \(1, \ldots, n - 1\) are strong liars, bases to which \(n\) is a strong pseudoprime.

So there are no strong analogues of Carmichael numbers. The Miller–Rabin test picks \(k\) random bases; if \(n\) is composite, the chance that all \(k\) are liars is at most \(4^{-k}\). With \(k = 40\) that is less than \(10^{-24}\), far smaller than the chance of a hardware error during the computation. Each round costs about \(\log_2 n\) modular multiplications, so testing a 600-digit number takes milliseconds. This is the test that cryptographic software uses to generate primes. For numbers below \(2^{64}\), testing the twelve bases 2, 3, 5, …, 37 is known to give a correct answer every time, so the test becomes deterministic.

Gary Miller showed in 1976 that if the generalized Riemann hypothesis is true, testing all bases up to \(2(\ln n)^2\) suffices, giving a deterministic polynomial-time test. In 2002 Manindra Agrawal and his students Neeraj Kayal and Nitin Saxena, at the Indian Institute of Technology Kanpur, found a deterministic polynomial-time test needing no unproved hypothesis, the AKS test, based on the polynomial congruence \((x + a)^n \equiv x^n + a \pmod n\), which holds for all \(a\) exactly when \(n\) is prime. It settled a famous theoretical question, although Miller–Rabin remains faster in practice.

A probable prime can also be certified. Lucas observed that \(n\) is prime if some \(a\) has order exactly \(n - 1\) modulo \(n\), that is, \(a^{n-1} \equiv 1\) but \(a^{(n-1)/q} \not\equiv 1\) for every prime \(q \mid n - 1\): then \(\varphi(n) \ge n - 1\), forcing \(n\) prime. This is the primitive root test of lesson 8 run in reverse. Given the factorization of \(n - 1\), anyone can check such a certificate quickly.

Factoring by differences of squares

Now the hard half. Fermat proposed factoring \(n\) by writing it as a difference of squares: if \(n = x^2 - y^2\), then \(n = (x - y)(x + y)\). Start at \(x = \lceil\sqrt n\rceil\) and increase \(x\) until \(x^2 - n\) is a perfect square. For \(n = 8051\): \(\lceil\sqrt{8051}\rceil = 90\) and \(90^2 - 8051 = 49 = 7^2\), so \(8051 = (90 - 7)(90 + 7) = 83 \cdot 97\). Fermat's method is fast when the two factors are close together, since then \(x\) is close to \(\sqrt n\), and slow otherwise. That is why RSA keys use primes that are not too close. Every odd \(n = ab\) is such a difference, with \(x = (a + b)/2\) and \(y = (b - a)/2\).

The modern methods generalize this idea. If we find \(x, y\) with \(x^2 \equiv y^2 \pmod n\) but \(x \not\equiv \pm y\), then \(n\) divides \((x - y)(x + y)\) without dividing either factor, so \(\gcd(x - y, n)\) is a proper factor. The quadratic sieve (Pomerance, 1981) and the number field sieve (Pollard and others, around 1990) build such congruences by multiplying together many numbers whose factorizations involve only small primes, using linear algebra modulo 2 to choose which to combine. The number field sieve is the fastest known general method; in 2020 it factored RSA-250, a 250-digit number, using about 2,700 core-years of computing time. Its running time grows faster than any power of the number of digits, though slower than exponentially.

Pollard's rho method

In 1975 John Pollard found a method of startling simplicity that finds a prime factor \(p\) of \(n\) in roughly \(\sqrt p\) steps rather than \(p\). Choose the map \(f(x) = x^2 + 1 \bmod n\) and iterate it from \(x_0 = 2\). The sequence modulo the unknown prime \(p\) lives in a set of only \(p\) values, so it must eventually repeat, and from then on it cycles. Its shape, a tail leading into a loop, looks like the Greek letter rho.

The birthday paradox explains why the repetition comes soon. Among about \(1.2\sqrt p\) random values from a set of size \(p\), the chance that two coincide is about one half, just as 23 people are likely to include two with the same birthday. The map \(x^2 + 1\) behaves enough like a random map that the tail and cycle together have length on the order of \(\sqrt p\). The figure shows the sequence modulo 61: a tail of 3 and a cycle of 10, compared with \(\sqrt{61} \approx 7.8\).

We do not know \(p\), but we can detect the repetition anyway. If \(x_i \equiv x_j \pmod p\) with \(x_i \ne x_j\), then \(\gcd(x_i - x_j, n)\) is divisible by \(p\) and is usually a proper factor of \(n\). Floyd's cycle-finding trick avoids storing the sequence: run two copies, a tortoise \(x_i\) and a hare \(x_{2i}\) moving twice as fast, and compute \(\gcd(x_{2i} - x_i, n)\) at each step. Once both are inside the cycle the hare gains one step per iteration and catches the tortoise after at most one lap. For \(n = 8051\): \(x_1 = 5\), \(x_2 = 26\), \(\gcd(21, 8051) = 1\); \(x_2 = 26\), \(x_4 = 7474\), \(\gcd(7448, 8051) = 1\); \(x_3 = 677\), \(x_6 = 871\), \(\gcd(194, 8051) = 97\). Three steps factor \(8051 = 83 \cdot 97\).

Pollard's rho finds factors of 10 to 20 digits easily and famously factored the Fermat number \(F_8\) in 1980, revealing a 16-digit prime factor. Pollard's other method, the \(p - 1\) method, computes \(\gcd(a^M - 1, n)\) for \(M\) a product of many small prime powers; it succeeds quickly when \(p - 1\) has only small prime factors, because then \(p - 1 \mid M\) and \(a^M \equiv 1 \pmod p\) by Fermat. Lenstra's elliptic curve method (1985) generalizes it by replacing the group of units modulo \(p\), of order \(p - 1\), with the group of points on an elliptic curve, whose order varies as the curve changes.

Why the gap matters

Multiplying two 300-digit primes takes microseconds; recovering them from the product would, with the best known algorithm, take thousands of years of the world's computing power. No one has proved that factoring is hard, and a large enough quantum computer running Shor's algorithm of 1994 could factor quickly; no such machine yet exists, and new cryptographic standards are already being adopted in preparation. For now the asymmetry between primality testing and factoring is the foundation on which the RSA system of the next lesson rests.

Testing and factoring by computer

The program counts base-2 pseudoprimes below a million, finds Carmichael numbers with Korselt's criterion, implements Miller–Rabin and counts strong liars, and factors numbers with Fermat's method and Pollard's rho.

import math, random

N = 10**6
sieve = bytearray([1]) * (N + 1); sieve[0:2] = b"\x00\x00"
for p in range(2, math.isqrt(N) + 1):
    if sieve[p]:
        sieve[p*p::p] = bytearray(len(range(p*p, N + 1, p)))
psp = [n for n in range(3, N, 2) if not sieve[n] and pow(2, n - 1, n) == 1]
print("primes below 10^6:", sum(sieve), " base-2 pseudoprimes:", len(psp), " first five:", psp[:5])

def factor(n):
    f, d = {}, 2
    while d * d <= n:
        while n % d == 0:
            f[d] = f.get(d, 0) + 1; n //= d
        d += 1
    if n > 1:
        f[n] = f.get(n, 0) + 1
    return f

def korselt(n):
    f = factor(n)
    return len(f) > 1 and all(e == 1 for e in f.values()) and all((n - 1) % (p - 1) == 0 for p in f)
print("Carmichael numbers below 10^4:", [n for n in range(3, 10**4, 2) if korselt(n)])

def strong_probable_prime(n, a):
    d, s = n - 1, 0
    while d % 2 == 0:
        d //= 2; s += 1
    x = pow(a, d, n)
    if x in (1, n - 1):
        return True
    for _ in range(s - 1):
        x = x * x % n
        if x == n - 1:
            return True
    return False

def miller_rabin(n, rounds=40):
    if n < 4:
        return n in (2, 3)
    if n % 2 == 0:
        return False
    return all(strong_probable_prime(n, random.randrange(2, n - 1)) for _ in range(rounds))

random.seed(7)
print("561 strong pseudoprime to base 2:", strong_probable_prime(561, 2))
for n in [561, 1105, 1729, 2047, 8911]:
    liars = sum(strong_probable_prime(n, a) for a in range(1, min(n, 20000)))
    print(f"n={n}: strong liars among first {min(n, 20000) - 1} bases: {liars}")
big = 2**521 - 1
print("2^521 - 1 probably prime:", miller_rabin(big), " 2^523 - 1:", miller_rabin(2**523 - 1))

def fermat_factor(n):
    x = math.isqrt(n)
    if x * x < n:
        x += 1
    while True:
        y2 = x * x - n
        y = math.isqrt(y2)
        if y * y == y2:
            return x - y, x + y
        x += 1

def pollard_rho(n, c=1):
    f = lambda x: (x * x + c) % n
    x = y = 2; d = 1; steps = 0
    while d == 1:
        x, y = f(x), f(f(y)); steps += 1
        d = math.gcd(abs(x - y), n)
    return d, steps

print("Fermat: 8051 =", fermat_factor(8051), " 1000009 =", fermat_factor(1000009))
print("rho on 8051:", pollard_rho(8051))
n = 1000000007 * 998244353
print("rho on a product of two 10-digit primes:", pollard_rho(n))
F6 = 2**64 + 1
d, steps = pollard_rho(F6)
print("rho on F6 = 2^64 + 1 found", d, "after", steps, "steps; cofactor", F6 // d, "prime:", miller_rabin(F6 // d))

primes below 10^6: 78498  base-2 pseudoprimes: 245  first five: [341, 561, 645, 1105, 1387]
Carmichael numbers below 10^4: [561, 1105, 1729, 2465, 2821, 6601, 8911]
561 strong pseudoprime to base 2: False
n=561: strong liars among first 560 bases: 10
n=1105: strong liars among first 1104 bases: 30
n=1729: strong liars among first 1728 bases: 162
n=2047: strong liars among first 2046 bases: 242
n=8911: strong liars among first 8910 bases: 1782
2^521 - 1 probably prime: True  2^523 - 1: False
Fermat: 8051 = (83, 97)  1000009 = (293, 3413)
rho on 8051: (97, 3)
rho on a product of two 10-digit primes: (998244353, 13433)
rho on F6 = 2^64 + 1 found 274177 after 808 steps; cofactor 67280421310721 prime: True

The counts of strong liars are small fractions of the bases, as the Rabin–Monier theorem promises, while every base coprime to a Carmichael number is a Fermat liar. The Mersenne number \(2^{521} - 1\), proved prime by Raphael Robinson in 1952 on one of the first electronic computers, passes 40 rounds in a moment. Pollard's rho splits a product of two primes near a billion in about 13,000 steps, fewer even than the square root of the smaller factor, about 31,600, and it finds the factor 274177 of the Fermat number \(F_6 = 2^{64} + 1\), first found by Landry in 1880 after months of hand calculation, in a few hundred steps.