iv. Primality Testing and Factorization

Task	Method	Cost
one $n\le 10^{12}$	trial division by $2,3$ and $6k\pm 1$	$O(\sqrt{n})$
one $n<2^{64}$	deterministic Miller–Rabin (7 fixed bases)	$O(\log n)$ multiplications per base
factor one $n<2^{64}$	Pollard–Brent rho + Miller–Rabin	expected $\tilde{O}(n^{1/4})$
factor many $n\le 10^7$	smallest-prime-factor sieve (Sieves)	$O(\log n)$ each
all primes in a range	(segmented) sieve	$O(n\log \log n)$
Mersenne numbers $2^p-1$	Lucas–Lehmer	$\tilde{O}(p^2)$
cryptographic sizes, need a certificate	ECPP	heuristic $\tilde{O}({\log }^{4}n)$
theoretical interest only	AKS	$\tilde{O}({\log }^{6}n)$, no unproven hypotheses needed

Smaller base sets suffice for smaller bounds

$\{2,3,5,7\}$ is deterministic for $n<3,215,031,751$; $\{2,7,61\}$ for $n<4.8\cdot 10^9$; $\{2,3,5,7,11,13,17\}$ for $n<3.4\cdot 10^{14}$. Use a shorter set when the input is bounded — fewer modpows per test.

Three recurring bugs

1. $n=1$ is not prime, and $n<2$ must be handled before anything else.
2. i * i <= n overflows when $n$ is near $2^{63}$ — use i <= n / i or __int128.
3. Fermat’s test alone is wrong: Carmichael numbers (561, 1105, 1729, …) satisfy $a^{n-1}\equiv 1$ for every coprime $a$. Korselt: $n$ is Carmichael iff it is squarefree, composite and $p-1\mid n-1$ for every prime $p\mid n$. Miller–Rabin’s square-root check is exactly what fixes this.

using u64 = unsigned long long; using u128 = __uint128_t;
u64 mulm(u64 a, u64 b, u64 m) { return (u128)a * b % m; }
u64 powm(u64 a, u64 e, u64 m) { u64 r = 1 % m; for (a %= m; e; e >>= 1, a = mulm(a, a, m)) if (e & 1) r = mulm(r, a, m); return r; }
bool isPrime(u64 n) {                                     // deterministic for all 64-bit n
    if (n < 2) return false;
    for (u64 p : {2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37}) if (n % p == 0) return n == p;
    int s = __builtin_ctzll(n - 1); u64 d = (n - 1) >> s;
    for (u64 a : {2ULL, 325ULL, 9375ULL, 28178ULL, 450775ULL, 9780504ULL, 1795265022ULL}) {
        if (a % n == 0) continue;
        u64 x = powm(a, d, n);
        if (x == 1 || x == n - 1) continue;
        bool comp = true;
        for (int i = 1; i < s && comp; i++) { x = mulm(x, x, n); if (x == n - 1) comp = false; }
        if (comp) return false;                            // a witnesses that n is composite
    }
    return true;
}
u64 pollard(u64 n) {                                      // a nontrivial factor of an odd composite n (Brent, batched gcd)
    static mt19937_64 rng(chrono::steady_clock::now().time_since_epoch().count());
    if (n % 2 == 0) return 2;
    for (;;) {
        u64 c = rng() % (n - 1) + 1, y = rng() % n, x = y, ys = y, g = 1, q = 1;
        auto f = [&](u64 v) { return (mulm(v, v, n) + c) % n; };
        for (u64 r = 1; g == 1; r <<= 1) {
            x = y; for (u64 i = 0; i < r; i++) y = f(y);
            for (u64 k = 0; k < r && g == 1; k += 128) {
                ys = y;
                for (u64 i = 0; i < min<u64>(128, r - k); i++) y = f(y), q = mulm(q, x > y ? x - y : y - x, n);
                g = std::gcd(q, n);
            }
        }
        if (g == n) do ys = f(ys), g = std::gcd(x > ys ? x - ys : ys - x, n); while (g == 1);   // replay the last block
        if (g != n) return g;                              // otherwise the walk closed without splitting: new c
    }
}
void factorRec(u64 n, vector<u64>& out) {
    if (n == 1) return;
    if (isPrime(n)) { out.push_back(n); return; }
    u64 d = pollard(n); factorRec(d, out); factorRec(n / d, out);
}
vector<pair<u64, int>> factorize(u64 n) {                 // sorted (prime, exponent)
    vector<u64> ps;
    for (u64 p : {2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37}) while (n % p == 0) ps.push_back(p), n /= p;
    factorRec(n, ps); sort(all(ps));
    vector<pair<u64, int>> res; for (u64 p : ps) { if (sz(res) && res.back().first == p) res.back().second++; else res.push_back({p, 1}); }
    return res;
}
vector<u64> divisorsFrom(const vector<pair<u64, int>>& f) {   // all divisors (unsorted), d(n) of them
    vector<u64> ds{1};
    for (auto [p, e] : f) { int m = sz(ds); u64 pk = 1; for (int k = 1; k <= e; k++) { pk *= p; for (int i = 0; i < m; i++) ds.push_back(ds[i] * pk); } }
    return ds;
}

From $n=\prod p_i^{e_i}$	Value
number of divisors $d(n)$	$\prod (e_i+1)$
sum of divisors $\sigma (n)$	$\prod \frac{p_i^{e_i+1}-1}{p_i-1}$
Euler $\phi (n)$	$n\prod (1-1/p_i)$
Carmichael $\lambda (n)$ (exponent of ${\mathbb{Z}}_n^{\ast }$)	$\text{lcm}\lambda (p_i^{e_i})$, $\lambda (p^e)=p^{e-1}(p-1)$ except $\lambda (2^e)=2^{e-2}$ for $e\ge 3$
perfect square / $k$-th power	all $e_i$ even / divisible by $k$
primitive root exists	$n\in \{1,2,4,p^k,2p^k\}$ (Primitive Roots)

Proof: Miller–Rabin

Write $n-1=2^{s}d$ with $d$ odd. If $n$ is an odd prime, ${\mathbb{Z}}_n$ is a field, so $x^2=1$ has only the roots $\pm 1$. By Fermat $a^{2^{s}d}=1$, and walking the sequence $a^d,a^{2d},\dots ,a^{2^{s}d}$ backwards from the final 1, each earlier term is a square root of the next: either the sequence is all ones from the start ($a^d=1$), or the term just before the first 1 is $-1$. So a prime never produces the “composite” verdict, and any base that does is a certificate of compositeness. For odd composite $n$, at most a quarter of the bases in $[1,n-1]$ fail to witness (Rabin–Monier), and the seven bases used here were verified by exhaustive search to have no common strong pseudoprime below $2^{64}$; the trial division by primes $\le 37$ handles small $n$ and bases divisible by $n$ are skipped (they carry no information).

Proof sketch: Pollard's rho runs in about $n^{1/4}$

Let $p$ be the smallest prime factor of $n$ ($p\le \sqrt{n}$). Modulo $p$, the map $x\mapsto x^2+c$ behaves like a random function, so the sequence repeats modulo $p$ after $O(\sqrt{p})$ steps (birthday paradox, Probability), i.e. $x_i\equiv x_j(\mathrm{mod}p)$ and so $p\mid \gcd (x_i-x_j,n)$. Brent’s cycle finding compares $x_{2^k}$ with the following $2^k$ values, which detects any cycle within a constant factor of its length. The gcd is $n$ only if the walk also closed modulo every other prime factor at the same moment, which happens with small probability; then a new $c$ is drawn. Batching 128 differences into one product loses nothing: if the product’s gcd with $n$ is $n$, the block is replayed one step at a time. So a factor appears after expected $O(\sqrt{p})=O(n^{1/4})$ multiplications (heuristically; no proof for this specific map is known).

Beyond Pollard rho: Baillie–PSW (one Miller–Rabin base-2 test plus a Lucas test, no known counterexample) is a fast primality heuristic; Pollard’s Kangaroo adapts the same random-walk idea to discrete logarithms; Lenstra’s ECM factors numbers with a medium-sized factor (up to ~60 digits) faster than rho; the quadratic sieve / GNFS are what’s needed for genuinely large semiprimes, where rho is hopeless — none of these are typically hand-implemented in a contest, but recognise the names.

Proof: divisor formulas and Korselt

By unique factorisation a divisor is $\prod p_i^{a_i}$ with $0\le a_i\le e_i$, giving $d(n)$; summing, $\sigma (n)=\prod (1+p_i+\cdots +p_i^{e_i})$, and each factor is a geometric series. $\phi$ is inclusion–exclusion over the primes dividing $n$ (Arithmetic Functions). Korselt (if): for squarefree $n$ with $p-1\mid n-1$, $a^{n-1}\equiv 1(\mathrm{mod}p)$ for every $p\mid n$ coprime to $a$, hence modulo $n$ by CRT. (Only if): if $p^2\mid n$, take $a$ of order $p$ modulo $p^2$ ($a=1+p$); $a^{n-1}\equiv 1$ would need $p\mid n-1$, impossible as $p\mid n$. If $p-1\nmid n-1$, a primitive root $g$ modulo $p$ (lifted by CRT to be coprime to $n$) has $g^{n-1}\not\equiv 1(\mathrm{mod}p)$.

Wilson’s theorem and factorials modulo $n$

$(n-1)!\equiv -1(\mathrm{mod}n)$ iff $n$ is prime; for composite $n>4$, $(n-1)!\equiv 0$, and $3!\equiv 2(\mathrm{mod}4)$. Too slow as a test, but it evaluates products like $(p-1)!$, ${\prod }_{i<p,p\nmid i}i$ and $\left(\genfrac{}{}{0ex}{}{p-1}{k}\right)\equiv (-1{)}^k$ in factorial-mod-$p$ problems (Factorials mod p).

Proof: Wilson

For prime $p$, pair each $a\in \{1,\dots ,p-1\}$ with its inverse; only $a=\pm 1$ are self-inverse ($a^2\equiv 1⟺p\mid (a-1)(a+1)$), so the product is $1\cdot (p-1)\equiv -1$. For composite $n=ab$ with $1<a\le b<n$: if $a<b$ both appear in $(n-1)!$, so $n\mid (n-1)!$; if $n=a^2$ with $a>2$, then $a$ and $2a$ are distinct and $<n$, so $a^2\mid (n-1)!$; $n=4$ gives $6\equiv 2$. In both composite cases the product is not $-1$.

Useful consequences: $(p-2)!\equiv 1(\mathrm{mod}p)$ (drop the $-1$ factor from Wilson), and more generally a factorial near $p$ converts to a small one via $(p-1-k)!\equiv \frac{-1}{(-1{)}^{k+1}k!}(\mathrm{mod}p)$ — handy when a formula needs $(p-k)!$ for small $k$. $(p-1)!\equiv -1(\mathrm{mod}{p}^2)$ holds only for the three known Wilson primes 5, 13, 563 (none other known below $2\cdot 10^{13}$). Defining $(n!{)}_p$ as the product of integers $\le n$ not divisible by $p$, Wilson gives $(n!{)}_p\equiv (-1{)}^{\lfloor n/p\rfloor }(n\mathrm{mod}p)!(\mathrm{mod}p)$ — the engine behind the generalised Lucas theorem (Factorials mod p).

Special forms worth checking before general factorization

Form	Approach
$n=a^b$ (perfect power)	check each $b$ from 2 to ${\log }_{2}n$, taking an integer $b$-th root
$n=x^2-y^2$ (difference of squares)	Fermat’s method: fast when $n$‘s two factors are close to $\sqrt{n}$, useless otherwise
$n=2^{k}m$	strip powers of two with __builtin_ctzll first
factorials and binomials	never factor the number itself — use Legendre’s formula for the exponent of $p$ in $n!$