Why this matters when you build with AI
An AI coding assistant can write a primality test in ten seconds. The code compiles, the first three examples pass, and it still returns the wrong answer for numbers near 10^18 because the multiplication inside it overflows. That kind of bug does not crash. It returns a plausible true or false, and you only find it when a number that should be prime is rejected.
This lesson gives you two tools for numbers up to 2^64 (about 1.8 × 10^19). The first, Miller-Rabin, answers "is this number prime?" fast. The second, Pollard rho, finds a factor of a number that is not prime. Together they let you factorise numbers that are far too large for trial division. Trial division means dividing by every number up to the square root, which takes about 10^9 steps for a number near 10^18.
You will also learn the checks that catch AI-written versions of these algorithms. Those checks matter more than the algorithms, because the bugs look correct on small inputs.
The problem in plain words
A number is prime if its only divisors are 1 and itself. Two questions are different:
- Primality test: is
nprime? Answer yes or no. - Factorisation: what are the prime factors of
n?
For a number with 18 digits, a primality test can be fast even when factorisation is hard. For example, 1000000016000000063 equals 1000000007 times 1000000009. Both factors have 10 digits. Trial division would have to try about 10^9 divisors to find 1000000007. Pollard rho finds it in roughly 3 × 10^4 steps, tens of thousands of times less work.
A fast test that is wrong: the Fermat test
The simplest fast test uses Fermat's little theorem. If p is prime and a is not a multiple of p, then a^(p-1) mod p equals 1. Check this for a random a. If the answer is not 1, n is composite for sure. If it is 1, n is probably prime.
The problem is that some composite numbers pass for many bases. The number 561 = 3 × 11 × 17 passes for every base that is coprime to it. Such numbers are called Carmichael numbers. A test that stops at base 2 has a second problem: 3215031751 = 151 × 751 × 28351 is a strong pseudoprime to bases 2, 3, 5 and 7, which means it passes the stronger test that Miller-Rabin uses. We checked it: our base-2-only check returns "probably prime" for 3215031751, and the full test rejects it.
So a single base is not enough. Miller-Rabin adds a stronger check and a set of fixed bases that make it deterministic for a known range.
Miller-Rabin: the idea with a worked example
Write n - 1 as d * 2^s, where d is odd. For n = 561, n - 1 = 560 = 35 * 2^4, so d = 35 and s = 4.
Take a base a, here a = 2. Compute x = a^d mod n. Then square x up to s - 1 times, checking after each step.
- If
xis 1 orn - 1at the start, the base passes. - If some square equals
n - 1, the base passes. - If neither happens,
ais a witness: it provesnis composite.
For 561 with base 2, the values are:
2^35 mod 561 = 263
squaring: 263 -> 166 -> 67 -> 1The value never equals 560 (n - 1), and it reaches 1 after 67. A square root of 1 that is not 1 or n - 1 proves that n is composite, so base 2 is a witness. We ran this chain and it matches the output of the code below.
The key fact is this: for a prime n, the only square roots of 1 modulo n are 1 and n - 1. So a chain that reaches 1 without passing through n - 1 cannot come from a prime.
Which bases to use for 64-bit numbers
A test with random bases is correct with high probability, but this lesson uses fixed bases so the answer is always the same. The set 2, 325, 9375, 28178, 450775, 9780504, 1795265022 is attributed to Jim Sinclair. We found it in several library sources (a Rust crate, a Rust factorisation module and a MATLAB paper that cites it). Those sources are secondary, so treat the claim as the published result, not as something we proved.
What we did check: the code agrees with a sieve for every n up to 10^7 (664579 primes, 0 mismatches), and it agrees with trial division on 200000 random numbers up to 10^12. We did not check every 64-bit number, because that would take far longer than this lesson allows.
The skip rule
A base a that is a multiple of n tells you nothing. For example, 28178 = 2 × 73 × 193. For n = 73, 28178 mod 73 is 0, so the base is useless and the test must skip it.
We wrote a version without the skip rule to see what happens. It declares the prime 73 composite. The first wrong answer below 10^6 is 73, and the version with the skip rule makes no mistakes in that range. Always write if (a % n == 0) continue;.
Overflow in the multiplication
Every step multiplies two numbers below n. For n near 10^18, the product can be near 10^36, which does not fit in 64 bits. The product wraps around and the result is garbage. Use a 128-bit product, (unsigned __int128)a * b % m, which is supported by g++ and clang on 64-bit targets.
We ran the same test with the plain 64-bit product a * b % m on 20000 odd numbers near 10^18. The correct test found 956 primes. The plain version found 0. Every answer was wrong, and nothing crashed.
The standard code: Miller-Rabin and Pollard-Brent
Pollard rho is the second half of this lesson. Its idea is covered in the next section, so for now read the whole file once. It is a complete program with a main that factorises the numbers you type in. It compiles with g++ -O2 -std=c++17 -Wall.
// Deterministic Miller-Rabin and Pollard-Brent rho for 64-bit n.
#include <bits/stdc++.h>
using namespace std;
typedef unsigned long long u64;
typedef unsigned __int128 u128;
// (a * b) mod m without overflow. a and b must be below m.
u64 mulmod(u64 a, u64 b, u64 m) {
return (u64)((u128)a * b % m);
}
u64 powmod(u64 a, u64 e, u64 m) {
u64 r = 1 % m;
a %= m;
while (e) {
if (e & 1) r = mulmod(r, a, m);
a = mulmod(a, a, m);
e >>= 1;
}
return r;
}
// One Miller-Rabin round with base a. n is odd, n-1 = d * 2^s with d odd.
// Returns true when n passes the round (n is probably prime).
bool mr_round(u64 n, u64 a, u64 d, int s) {
a %= n;
if (a == 0) return true; // base is a multiple of n: no information
u64 x = powmod(a, d, n);
if (x == 1 || x == n - 1) return true;
for (int i = 1; i < s; i++) {
x = mulmod(x, x, n);
if (x == n - 1) return true;
}
return false; // a is a witness: n is composite
}
// Deterministic for every n < 2^64 with these 7 bases (set credited to Jim Sinclair).
bool is_prime(u64 n) {
if (n < 2) return false;
static const u64 small[] = {2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37};
for (u64 p : small) {
if (n % p == 0) return n == p;
}
u64 d = n - 1;
int s = 0;
while ((d & 1) == 0) { d >>= 1; s++; }
static const u64 bases[] = {2, 325, 9375, 28178, 450775, 9780504, 1795265022};
for (u64 a : bases) {
if (!mr_round(n, a, d, s)) return false;
}
return true;
}
// Pollard rho with Brent's cycle finding. Returns a non-trivial factor of composite n.
u64 pollard_brent(u64 n) {
if (n % 2 == 0) return 2;
static mt19937_64 rng(20261108);
while (true) {
u64 y = rng() % (n - 1) + 1;
u64 c = rng() % (n - 1) + 1;
const u64 m = 128;
u64 g = 1, r = 1, q = 1, x = 0, ys = 0;
auto f = [&](u64 v) -> u64 {
u128 t = (u128)mulmod(v, v, n) + c; // the sum can exceed 2^64 when n is large
return (u64)(t % n);
};
do {
x = y;
for (u64 i = 0; i < r; i++) y = f(y);
u64 k = 0;
while (k < r && g == 1) {
ys = y;
u64 lim = min(m, r - k);
for (u64 i = 0; i < lim; i++) {
y = f(y);
u64 diff = x > y ? x - y : y - x;
q = mulmod(q, diff, n);
}
g = std::gcd(q, n);
k += m;
}
r *= 2;
} while (g == 1);
if (g == n) {
// The batch product hit n. Redo the last batch one step at a time.
do {
ys = f(ys);
u64 diff = x > ys ? x - ys : ys - x;
g = std::gcd(diff, n);
} while (g == 1);
}
if (g != n) return g; // otherwise retry with another c
}
}
void factor_rec(u64 n, vector<u64>& out) {
if (n == 1) return;
if (is_prime(n)) { out.push_back(n); return; }
u64 d = pollard_brent(n);
factor_rec(d, out);
factor_rec(n / d, out);
}
// Full factorisation: trial division by small primes, then MR and rho.
vector<pair<u64, int>> factorise(u64 n) {
vector<u64> primes;
for (u64 p = 2; p < 1000 && p * p <= n; p++) {
while (n % p == 0) { primes.push_back(p); n /= p; }
}
if (n > 1) factor_rec(n, primes);
sort(primes.begin(), primes.end());
vector<pair<u64, int>> res;
for (u64 p : primes) {
if (!res.empty() && res.back().first == p) res.back().second++;
else res.push_back({p, 1});
}
return res;
}
int main() {
int t;
if (!(cin >> t)) return 0;
while (t--) {
u64 n;
cin >> n;
auto f = factorise(n);
cout << n << ":";
for (auto [p, e] : f) {
cout << " " << p;
if (e > 1) cout << "^" << e;
}
cout << "\n";
}
return 0;
}Pollard rho: the idea with a worked example
Pick a function that looks random, here f(x) = x^2 + c mod n. Start at a value and repeat: x, f(x), f(f(x)), .... Modulo any prime factor p of n, the sequence eventually repeats, because there are only p possible values. When it repeats, the difference between two positions in the sequence is a multiple of p. Taking gcd(difference, n) then reveals p. The gcd is the largest number that divides both values, and the Euclidean algorithm computes it fast.
Here is a trace with n = 8051 = 83 × 97, f(x) = x^2 + 1, and x starting at 2. The slow pointer x moves one step per round. The fast pointer y moves two steps. We print gcd(|x - y|, n).
round x y gcd(|x-y|, 8051)
1 5 26 1
2 26 7474 1
3 677 871 97At round 3 the gcd is 97, which splits 8051 into 97 and 83. The sequence modulo 97 repeats first, so the difference becomes a multiple of 97 before it becomes a multiple of 83.
Brent's improvement keeps one saved value x and moves a second value y forward. It checks for a repeat only at powers of two, and it multiplies the differences together so there is one gcd per batch of 128 steps rather than one per step. R. P. Brent, "An improved Monte Carlo factorization algorithm", BIT Numerical Mathematics, volume 20, 1980, describes the method. We confirmed that citation on Brent's publication page and a bibliographic listing, not by reading the full paper. The batching is why the code above is fast enough for 10^18.
Complexity
Miller-Rabin costs about one modular exponentiation per base. Each exponentiation has about 64 squarings, so one primality test is about 7 × 64 = 448 multiplications of 128-bit products, plus a few more for the squaring chain. Each step is a constant number of operations, so the test is cheap.
Pollard rho is expected to find a factor p in about the square root of p steps. For a semiprime near 10^18 whose smallest factor is near 10^9, that is about 3 × 10^4 steps. We factorised 100 random products of two primes, each between 10^9 and 2 × 10^9, with the code above. It took about 0.1 seconds on the WSL2 machine used to write this lesson. Your times will differ with your machine and compiler.
Worked problem: factor the numbers in a batch
You receive t numbers, each at most 10^18. For each one, print its prime factorisation, with the primes in increasing order and repeats shown as powers.
Input and output with the program above, using the file we ran:
5
360
3215031751
1000000000000000000
1000000016000000063
999999999989360: 2^3 3^2 5
3215031751: 151 751 28351
1000000000000000000: 2^18 5^18
1000000016000000063: 1000000007 1000000009
999999999989: 999999999989Check the output yourself. Multiply the factors back together for each line: 151 × 751 × 28351 = 3215031751. Also check that each listed factor is prime. The program does not print the check; you should write one.
The program also passed two brute-force tests. It agreed with a sieve for every n up to 10^7, and its factorisations of 200000 random numbers up to 10^12 multiplied back to the original number, with all factors prime.
Common mistakes
Plain 64-bit multiplication. The symptom is a test that declares every number near 10^18 composite. In our run the wrong version found 0 primes among 20000 odd numbers near 10^18, while the correct test found 956. Use unsigned __int128 for the product.
No skip for bases that are multiples of n. The symptom is a prime reported as composite. The version without the skip rule declared 73 composite. Add if (a % n == 0) continue;.
Using only base 2. The symptom is a composite reported as prime. The number 3215031751 passes the base-2 test. Use the full base set.
Calling Pollard rho on a prime. Pollard rho never finds a factor of a prime. The program then loops forever. We ran pollard_brent(1000000007) with a 5-second timeout and it did not return. Always call is_prime first, as factor_rec does.
Trusting AI output without a test. AI assistants often give Miller-Rabin code that has one of the bugs above. They also often claim that the first few primes as bases are enough for 64-bit numbers. That claim is false for base 2 alone, as the example of 3215031751 shows. Test the code against a sieve before you trust it.
Checklist to keep
- Write
n - 1asd * 2^swithdodd before the test starts. - Multiply with
unsigned __int128for any modulus above 2^32. - Skip every base that is a multiple of
n. - Use the seven fixed bases for 64-bit inputs.
- Run
is_primebefore calling Pollard rho. - Check each factor with
is_primeand multiply the factors back ton. - Compare with a sieve on all numbers up to 10^6 or 10^7 before you trust it.
Summary
Miller-Rabin tests primality with a strong check that a single Fermat test does not give. A fixed set of seven bases makes the result deterministic for 64-bit inputs, as the published result states. Pollard rho with Brent's cycle finding finds a factor of a composite number in about the square root of that factor's size. Both need 128-bit multiplication, and both fail quietly when written carelessly.
Practice
Use this prompt when you ask an AI assistant to write the code yourself. Check the answer with the tests in this lesson.
Write a C++17 function bool is_prime(unsigned long long n) that is exact for every n below 2^64.
Use Miller-Rabin with a fixed set of bases. Multiply with unsigned __int128.
Skip any base that is a multiple of n. Do not use any library primality test.
Show the bases you used and explain why the test is exact for 64-bit inputs.A good answer uses the seven bases above, has the skip rule, and uses a 128-bit product. Check it against a sieve up to 10^6 and against 3215031751, which must be composite.
Practice 1: Trace 561 by hand
Compute 2^35 mod 561 by repeated squaring, then check the squaring chain for base 2. Is 561 composite or probably prime according to this base alone?
Show the answer
2^35 mod 561 = 263. Squaring gives 166, then 67, then 1. The chain never reaches 560, so base 2 is a witness and 561 is composite. (Use the square-and-multiply steps from the powmod function to get the same numbers.)
Practice 2: Find the first failure of the skip-free test
Write is_prime without the skip rule and run it for every n from 2 to 10^6 against a sieve. Report the first wrong answer.
Show the answer
The first wrong answer is 73. The base 28178 equals 2 × 73 × 193, so it is a multiple of 73 and gives no information. Adding if (a % n == 0) continue; removes the error.
Practice 3: Largest prime factor for many numbers
You receive up to 1000 numbers, each up to 10^18. Print the largest prime factor of each. Reuse factorise from the code above.
Hint: factorise returns the primes in increasing order, so the last entry of the result is the largest prime.
Show the answer
Replace the output loop in main with:
auto f = factorise(n);
cout << f.back().first << "\n";For n = 1, factorise returns an empty list, so guard that case with if (n == 1) { cout << 1 << "\n"; continue; } before the call. Test it on 1000000016000000063, which must print 1000000009.
Practice 4: Why rho must not run on a prime
Explain, in two sentences, why calling pollard_brent on a prime loops forever. Then describe the one check that prevents it.
Hint: think about what gcd(difference, p) can return when p is prime.
Show the answer
For a prime p, every difference between two sequence values is either 0 or a number below p, so gcd(difference, p) is 1 or p, and the loop never finds a non-trivial divisor. The check is is_prime(n) before calling pollard_brent, as factor_rec does.
Practice 5: Bound the 128-bit sum
In pollard_brent, the function f adds c to a value that is already below n. Explain why the sum must be computed in 128 bits when n is above 2^63.
Hint: write the largest possible value of mulmod(v, v, n) plus c.
Show the answer
Both mulmod(v, v, n) and c are below n, so their sum is below 2n. When n is above 2^63, 2n is above 2^64, so a 64-bit sum can wrap around. The code computes the sum as u128 and then takes it modulo n.