Lab 06 — Miller-Rabin & Pollard’s Rho

Goal

Deterministic primality testing for 64-bit integers, and full factorization of numbers up to ~10^18 in milliseconds.

Background

Fermat test: if n is prime, a^(n−1) ≡ 1 (mod n). Necessary, not sufficient — Carmichael numbers (561 = 3×11×17 is the smallest) pass for every a coprime to n. Infinitely many exist.

Miller-Rabin (strong pseudoprime test): write n−1 = d·2^s with d odd. For prime n, the sequence a^d, a^(2d), …, a^(2^s·d) must reach 1, and the only square roots of 1 mod a prime are ±1 — so either a^d ≡ 1 or some a^(2^r·d) ≡ −1. An a violating this is a witness: a certificate that n is composite. A random a is a witness for composite n with probability ≥ 3/4 — no Carmichael-style escape.

Deterministic for 64-bit: testing the first 12 primes {2,3,5,7,11,13,17,19,23,29,31,37} is correct for all n < 3.317×10^24 (Sorenson & Webster 2015) — covers all of uint64.

Pollard’s rho (1975): iterate x → x² + c mod n. Viewed mod p (smallest prime factor), the sequence enters a cycle after expected O(√p) steps — birthday paradox. A collision mod p but not mod n means gcd(|x−y|, n) is a nontrivial factor. Floyd (tortoise/hare) or Brent (backtracking powers of 2, ~24% fewer evals, batched gcds) detects it.

Interview Context

  • Codeforces / ICPC: standard tool; Project Euler: constant
  • Cryptography engineering: primality generation for RSA/DH keys is exactly this
  • CTF / security: factoring weak keys is a staple
  • Quant: occasional number-theory screens
  • Standard FAANG: near zero — at most “check if n is prime” (trial division suffices)

When to Skip This Topic

Skip if any of these are true:

  • You’re not targeting competitive programming, crypto, or security roles
  • You haven’t internalized modular exponentiation and gcd
  • Your pipeline never sees integers past 10^12 (trial division is fine there)

If FAANG-only: skip the whole lab, keep the trial-division one-liner.

Problem Statement

Fast factorization.

Given up to 100 integers, each up to 10^18, output the full prime factorization of each (primes with multiplicity, sorted). Products of two ~10^9 primes must work — trial division dies there.

Constraints

  • ≤ 100 queries, each 1 ≤ n ≤ 10^18
  • Total wall-clock: < 1 sec
  • Deterministic output required (no failure probability)

Clarifying Questions

  1. How large can n be? (10^18 — fits uint64; √n loop does not fit the time budget.)
  2. Is a probabilistic answer acceptable? (No — and it needn’t be: fixed witnesses are deterministic for 64-bit.)
  3. Output format? (Sorted prime list with multiplicity, e.g. 12 → [2, 2, 3].)
  4. Edge inputs? (1 → empty factorization; assume n ≥ 1, no negatives.)

Examples

Miller-Rabin on n = 221 (= 13×17). n−1 = 220 = 55·2², so d = 55, s = 2.

a = 174: 174^55 mod 221 = 47;  47² mod 221 = 220 = n−1 → passes (174 is a strong liar)
a = 137: 137^55 mod 221 = 188; 188² mod 221 = 205 ≠ ±1 → 137 is a witness: composite

One liar survives; the fixed 12-prime battery does not lie for any 64-bit n.

Pollard’s rho on 8051 = 83×97, f(x) = x²+1, x = y = 2:

step 1: x=5,   y=26   gcd(21, 8051)  = 1
step 2: x=26,  y=7474 gcd(7448, 8051) = 1
step 3: x=677, y=871  gcd(194, 8051)  = 97   → factor found

Brute Force

Trial division by 2, 3, 5, … up to √n. Fine to n ≈ 10^12 (10^6 iterations). For a semiprime near 10^18: 10^9 iterations per query — dead.

Brute Force Complexity

  • Time: O(√n) per query — 10^9 ops worst case
  • Space: O(1)

Optimization Path

  1. Strip primes < 100 by trial division (kills most mass cheaply).
  2. Miller-Rabin with the fixed 12-witness set → deterministic is_prime for uint64 in O(12 log³ n).
  3. If composite: Pollard’s rho finds a factor in expected O(n^{1/4}); recurse on both halves.
  4. Constant factors: Brent’s cycle detection + batching gcds (multiply ~128 differences mod n, one gcd).

Final Expected Approach

import math

WITNESSES = (2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37)

def is_prime(n):
    if n < 2:
        return False
    for p in WITNESSES:
        if n % p == 0:
            return n == p
    d, s = n - 1, 0
    while d % 2 == 0:
        d //= 2
        s += 1
    for a in WITNESSES:
        x = pow(a, d, n)
        if x == 1 or x == n - 1:
            continue
        for _ in range(s - 1):
            x = x * x % n
            if x == n - 1:
                break
        else:
            return False
    return True

def pollard_rho(n):
    if n % 2 == 0:
        return 2
    c = 1
    while True:
        x = y = 2
        d = 1
        while d == 1:
            x = (x * x + c) % n
            y = (y * y + c) % n
            y = (y * y + c) % n
            d = math.gcd(abs(x - y), n)
        if d != n:
            return d
        c += 1

def factor(n):
    if n == 1:
        return []
    if is_prime(n):
        return [n]
    d = pollard_rho(n)
    return sorted(factor(d) + factor(n // d))

is_prime must run before pollard_rho — rho never terminates on a prime. The small-prime loop doubles as the n ≤ 37 guard.

Data Structures

  • None beyond integers — the whole lab is arithmetic
  • Recursion stack for factor(): depth ≤ log₂ n ≈ 60

Correctness Argument

  • Miller-Rabin: for prime n, a^(n−1) ≡ 1 and repeated square-rooting from 1 can only pass through ±1; so the d·2^r ladder must show a^d ≡ 1 or hit −1. Any a breaking this proves n composite — no false “composite” verdicts, ever.
  • Determinism: the 12-prime witness set has been exhaustively verified correct below 3.317×10^24 > 2^64. No probability involved.
  • Rho: x_i mod p cycles with expected tail+cycle O(√p); Floyd guarantees some i with x_i ≡ x_{2i} (mod p), giving p | gcd(|x_i − x_{2i}|, n). gcd = n only if all prime factors collide simultaneously — then retry with the next c changes the sequence.
  • factor(): splits n into two smaller cofactors and recurses; termination since each part shrinks and primes are leaves.

Complexity

  • is_prime: O(12 log³ n) bit ops — microseconds via three-arg pow
  • pollard_rho: expected O(n^{1/4}) iterations — ≈ 31,600 for the worst 64-bit semiprime (p ≈ 10^9)
  • factor: O(n^{1/4} log n) expected; 100 worst-case queries well under 1 sec
  • Space: O(log n)

Implementation Requirements

  • Three-arg pow(a, d, n) only — never a**d % n
  • Inner squaring loop runs s−1 times, not s (a^d already checked)
  • Even n handled before rho; is_prime before rho, always
  • c advances deterministically 1, 2, 3, … on gcd = n — reproducible runs, no RNG
  • Return factors sorted with multiplicity

Tests

  • Primes: 2, 3, 5, 2^61−1, 999999999999999989 (largest prime < 10^18)
  • Carmichael: 561 → [3, 11, 17]; 41041 → [7, 11, 13, 41]
  • Strong-liar case: 221 with witness list must say composite
  • Semiprime: 10^18-scale product of two ~10^9 primes, e.g. 999999937 × 999999893
  • Edge: 1 → [], 4 → [2, 2], perfect powers: 2^60, 3^37
  • Stress: factor(n) product check + all-parts-prime check for 10^4 random n < 10^12 vs trial division
  • Performance: 100 random 18-digit numbers in < 1 sec

Follow-up Questions

  • Why is the fixed witness set deterministic — isn’t Miller-Rabin probabilistic? → Probabilistic only for random witnesses. Explicit computation found every strong pseudoprime to small prime bases: Jaeschke (1993) pushed the first 8 primes to 3.4×10^14; Sorenson–Webster (2015) verified the first 12 primes below 3.317×10^24. Jim Sinclair’s set {2, 325, 9375, 28178, 450775, 9780504, 1795265022} covers uint64 with only 7 rounds.
  • Why is rho expected O(n^{1/4})? → Let p be the smallest prime factor, p ≤ √n. The pseudorandom sequence mod p collides after expected O(√p) steps (birthday bound over p values). O(√p) ≤ O(n^{1/4}).
  • Why does rho fail on primes, and on n = 4? → On prime n the only “collision” is mod n itself, so gcd is always 1 or n — infinite loop; call is_prime first. n = 4: √p = √2 gives a degenerate 1-element cycle that x²+c may never expose — the even-n guard (and small-prime stripping) removes it.
  • Brent vs Floyd? → Brent keeps one saved point and checks x_i against x_{2^k} — ~24% fewer function evaluations. The real win is batching: accumulate ∏|x−y| mod n over ~128 steps, take one gcd; a 30–50% wall-clock cut since gcd is the expensive step.
  • Beyond 64 bits? → ECM digs out factors up to ~10^20 from arbitrarily large n (cost scales with the factor, not n); quadratic sieve is fastest for n up to ~100 digits; GNFS beyond that — RSA-250 (829 bits) fell to GNFS in 2020 after ~2700 core-years.
  • Relation to RSA? → RSA keygen is Miller-Rabin (find two ~1024-bit primes); RSA security is the absence of anything like rho at that scale — rho is n^{1/4}, i.e. 2^512 steps for a 2048-bit modulus. These two algorithms are both sides of the asymmetry.

Product Extension

  • Key-generation services: primality testing at scale (HSMs, TLS cert tooling)
  • Weak-key auditing: batch-factoring RSA moduli harvested from devices (Mining-your-Ps-and-Qs style)
  • Hash-table / hashing libraries: next-prime sizing
  • Number-theory backends: SymPy’s factorint, Pari/GP — same rho + MR core

Language/Runtime Follow-ups

  • Python: three-arg pow is C-speed; native bigints make this lab genuinely fast in pure Python — a rare case.
  • C++: mulmod of two uint64 overflows — use __int128; witnesses reduced mod n need the a % n == 0 guard.
  • Java: Math.multiplyHigh for mulmod, or BigInteger.isProbablePrime(certainty) as the library answer.
  • Rust: u128 intermediates; num-prime crate mirrors this exact design.

Common Bugs

  1. s vs s−1: squaring loop after the initial a^d check must run s−1 times; running s wrongly passes some composites through the final compare.
  2. Missing x == n−1 pre-check: if a^d ≡ n−1 immediately, the round passes; forgetting it yields false “composite” on primes.
  3. Witness divides n: for small n, a ∈ WITNESSES may equal or divide n — the small-prime loop must run first or 3 gets declared composite.
  4. Rho on a prime: infinite loop. is_prime gate is mandatory, not an optimization.
  5. gcd = n treated as failure: it means all factors collided at once — retry with next c, don’t return n.
  6. a**d % n: computes a full 10^17-bit power first; always three-arg pow.
  7. C++ mulmod overflow: (x*x) % n with x near 10^18 silently wraps; everything downstream is garbage.

Debugging Strategy

  • Cross-check is_prime against a sieve for all n < 10^6
  • factor(n): assert product of parts == n and is_prime on every part — self-verifying by construction
  • Deterministic c sequence means failures reproduce exactly; print (c, iteration count) on retry
  • Hand-trace 221 and 8051 (tables above) against your implementation’s intermediate values

Mastery Criteria

  • Deterministic is_prime from memory in ≤ 15 min
  • Full factor() with rho in ≤ 35 min
  • Reproduce the n = 221 witness/liar trace by hand
  • Derive the O(n^{1/4}) birthday bound in two sentences
  • State the deterministic witness set and its verified bound
  • Factor 100 random 10^18 integers in < 1 sec, first run