Prime Factorization, in depth: a tiered pipeline for every 64-bit integer, from trial division and deterministic Miller-Rabin to Montgomery arithmetic and the hand-off to ECM

The contract and what makes an input hard

Write the contract down first. The input is an integer n with 1 <= n < 2⁶⁴. The output maps each prime to its exponent, so 360 becomes {2: 3, 3: 2, 5: 1} and 1 becomes the empty map. The product of p^e over the map must equal n, and every key must be prime. Both properties are cheap to check, so check them.

It also helps to know what makes an input hard. The largest prime factor is never the problem, because it ends up as a cofactor that a primality test settles. The cost is set by the second-largest prime factor: trial division needs work proportional to it, and rho needs roughly its square root in steps. The worst 64-bit inputs are semiprimes with two factors near 2³².

The pipeline

The tiers are ordered by cost per factor found.

Tier 1, trial division by the 168 primes below 1000. About 92% of random integers have a prime factor below 1000. These small factors are exactly the ones where rho's fixed overhead per call is worst relative to their size.

Tier 2, a smallest-prime-factor table for any piece below 2²⁰. It factors such a piece in at most 20 lookups, and it is the right tool for bulk work such as computing Euler's totient over a range.

Tier 3, deterministic Miller-Rabin. A piece that passes is prime. Without this tier, rho would spin forever on a prime cofactor.

Tier 4, splitting. Check for a perfect power, then run Pollard-Brent with increasing constants until it returns a proper divisor d. Push d and m/d back on the stack.

Proving primality with a fixed base set

Miller-Rabin writes n - 1 = d · 2^s with d odd, computes a^d mod n for a base a, and squares it up to s - 1 times. For prime n, the sequence starts at 1 or reaches n - 1. A composite that passes for base a is a strong pseudoprime to that base. The smallest for base 2 is 2047 = 23 · 89, and 3215031751 fools bases 2, 3, 5 and 7 together. Background is in Miller-Rabin primality, in depth.

For bounded n, a fixed set of bases is a proof rather than a probability. The smallest strong pseudoprime to all of the first 12 primes, 2 through 37, is ψ₁₂ = 318,665,857,834,031,151,167,461, about 3.18 × 10²³, as computed by Sorenson and Webster. That is far above 2⁶⁴, about 1.8 × 10¹⁹. So those 12 bases decide primality for every unsigned 64-bit integer. Shorter 64-bit base sets exist, but they use large bases that must be reduced mod n, with the zero case skipped. The 12 small primes avoid that trap and are easy to audit.

Splitting: rho, perfect powers and overshoot

Rho iterates x → x² + c mod m and detects, with a gcd, when two values collide modulo an unknown factor p. Brent's version batches the gcds, so each step costs one or two modular multiplications. A batch can overshoot. The gcd jumps from 1 straight to m because every factor collided inside the same batch. The code then replays the batch one step at a time from a saved point, and if that also returns m, it moves on to the next constant c.

Perfect powers get their own check, because an integer k-th root is cheap. The float root round(m ** (1/k)) is approximate above 2⁵³, so test r - 1, r and r + 1 with exact integer powers.

The complete routine

The complete routine relies on Python's big integers, so (y * y + c) % n cannot overflow. The C version that needs care follows.

import math
from collections import Counter

SPF_LIMIT = 1 << 20
SPF = list(range(SPF_LIMIT))               # smallest prime factor table
for i in range(2, math.isqrt(SPF_LIMIT - 1) + 1):
    if SPF[i] == i:
        for j in range(i * i, SPF_LIMIT, i):
            if SPF[j] == j:
                SPF[j] = i
SMALL_PRIMES = [p for p in range(2, 1000) if SPF[p] == p]
MR_BASES = (2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37)

def is_prime(n):
    """Deterministic Miller-Rabin below 3.18e23, which covers all of uint64."""
    if n < 2:
        return False
    for p in MR_BASES:
        if n % p == 0:
            return n == p
    d, s = n - 1, 0
    while d % 2 == 0:
        d, s = d // 2, s + 1
    for a in MR_BASES:
        x = pow(a, d, n)
        if x in (1, n - 1):
            continue
        for _ in range(s - 1):
            x = x * x % n
            if x == n - 1:
                break
        else:
            return False
    return True

def perfect_power(n):
    for k in range(n.bit_length(), 1, -1):
        r = round(n ** (1.0 / k))
        for cand in (r - 1, r, r + 1):
            if cand > 1 and cand ** k == n:
                return cand
    return None

def brent(n, c, batch=128):
    y, r, q, g = 2, 1, 1, 1
    while g == 1:
        x = y
        for _ in range(r):
            y = (y * y + c) % n
        k = 0
        while k < r and g == 1:
            ys = y
            for _ in range(min(batch, r - k)):
                y = (y * y + c) % n
                q = q * abs(x - y) % n
            g = math.gcd(q, n)
            k += batch
        r *= 2
    if g == n:                             # overshot: replay one step at a time
        g = 1
        while g == 1:
            ys = (ys * ys + c) % n
            g = math.gcd(abs(x - ys), n)
    return g

def split(n):
    root = perfect_power(n)
    if root:
        return root
    for c in range(1, 64):
        d = brent(n, c)
        if 1 < d < n:
            return d
    raise RuntimeError(f"no split found for {n}")

def factorize(n):
    if not 1 <= n < 1 << 64:
        raise ValueError("n must be in [1, 2**64)")
    out = Counter()
    for p in SMALL_PRIMES:                 # tier 1
        if p * p > n:
            break
        while n % p == 0:
            out[p] += 1
            n //= p
    stack = [n] if n > 1 else []
    while stack:
        m = stack.pop()
        if m < SPF_LIMIT:                  # tier 2
            while m > 1:
                out[SPF[m]] += 1
                m //= SPF[m]
        elif is_prime(m):                  # tier 3
            out[m] += 1
        else:                              # tier 4
            d = split(m)
            stack += [d, m // d]
    return out

Test it against a slow oracle: naive trial division on random n below 10⁹, products of known 20- to 40-bit primes, the product check, and is_prime against a sieve and the published pseudoprimes. This exact code passed all of those.

Where the time goes in C: Montgomery multiplication

In C, a 64-bit modular multiplication needs a 128-bit product and a 128-by-64-bit remainder. (unsigned __int128)a * b % n compiles to a runtime division routine costing tens of cycles, and rho does one or two per step. Montgomery multiplication removes the division. Keep every value as aR mod n with R = 2⁶⁴. The product of two such values, divided by R, is again in that form, and dividing by R becomes a shift. The variant below subtracts high halves, which avoids a 129-bit intermediate. It needs n odd, which tier 1 guarantees.

#include <stdint.h>
typedef unsigned __int128 u128;

/* n^-1 mod 2^64 by Newton iteration: each step doubles the correct low bits. */
static uint64_t inv64(uint64_t n) {
    uint64_t x = n;                 /* correct to 3 bits for odd n */
    for (int i = 0; i < 5; i++) x *= 2 - n * x;
    return x;
}

/* Returns T / 2^64 mod n for T < n * 2^64; ninv = inv64(n). */
static inline uint64_t redc(u128 T, uint64_t n, uint64_t ninv) {
    uint64_t m    = (uint64_t)T * ninv;               /* low halves of T and m*n match */
    uint64_t t_hi = (uint64_t)(T >> 64);
    uint64_t mn_hi = (uint64_t)(((u128)m * n) >> 64);
    uint64_t t = t_hi - mn_hi;
    return t_hi < mn_hi ? t + n : t;
}

static inline uint64_t mont_mul(uint64_t a, uint64_t b, uint64_t n, uint64_t ninv) {
    return redc((u128)a * b, n, ninv);                /* a, b in Montgomery form */
}

Convert c into Montgomery form once, and y*y + c stays in the domain. The gcd works on Montgomery values directly, because gcd(xR - yR, n) = gcd(x - y, n). The same multiplication speeds up Miller-Rabin's exponentiation; see modular exponentiation, in depth. This REDC was checked against exact arithmetic on 200,000 random odd moduli.

Worked examples

Input 1: 2⁶⁴ - 1 = 18,446,744,073,709,551,615. Tier 1 strips 3, 5, 17, 257 and 641. The cofactor 439,125,228,929 fails Miller-Rabin and is not a perfect power. Brent with c = 1 returns 65,537 after 510 steps, a small multiple of √p = 256, as the birthday argument predicts. The table confirms 65,537 is prime, and 6,700,417 passes Miller-Rabin. The result is seven primes with exponent 1, and the product checks.

Input 2: 4,294,967,291 × 4,294,967,279 = 18,446,743,979,220,271,189. These are the two largest primes below 2³², the hardest kind of 64-bit input. Tier 1 does 168 useless divisions, Miller-Rabin rejects the number, and Brent needs 119,038 steps to return 4,294,967,291. That took about 70 ms in CPython. A Montgomery C version should need around a millisecond. That is an estimate, not a measurement.

A factor twice as long in bits costs about 2¹⁶ times as many trial divisions but only 2⁸ times as many rho steps. That square-root scaling is why rho is the default splitter at this size.

Beyond 64 bits: choosing a method

Above 64 bits, rho's square-root cost bites: a 40-digit factor would need about 10²⁰ steps. The standard ladder of methods is below.

Do not write these yourself. PARI/GP's factor, GMP-ECM, YAFU, msieve and CADO-NFS cover the ladder. The skill is routing: cheap methods first, ECM with growing bounds while the cofactor is large, and a sieve only when ECM has probably found everything small. Why large semiprimes matter is covered in RSA public-key cryptography.

Failure modes

• Probabilistic bases as a proof. One random base, or base 2 alone, accepts 2047. Use the fixed 12-base set and keep its bound in a test.
• Rho on a prime. It never returns. Test primality before splitting, and cap attempts so a bug ends in an error rather than a hang.
• Batch overshoot. Without the replay step, a batch whose gcd is n makes the routine retry forever.
• Overflow in C. y * y + c in uint64_t silently wraps, and adding c can overflow when n is near 2⁶⁴.
• Float roots. Trusting n ** 0.5 misses squares above 2⁵³.
• Unbounded inputs. A 200-bit number sent to a 64-bit routine runs for days. Validate the range and set a time budget.
• Unverified output. Pushing d but forgetting m/d gives a plausible wrong map. The product check catches it.

What to do next

1. Copy the routine and add the contract tests: oracle comparison, known semiprimes and the product check.
2. Assert that ψ₁₂ passes all 12 bases, so the edge of the guarantee is visible.
3. For speed, port tiers 3 and 4 to C with Montgomery arithmetic and fuzz them against the Python version.
4. Measure your input mix, then tune the trial-division bound and decide whether the 4 MB table earns its memory.
5. For inputs above 64 bits, call PARI/GP, YAFU or GMP-ECM instead of extending rho.
6. Put a size limit and a timeout on every factoring endpoint that untrusted callers reach.