Pollard's Rho Factorization, in depth: the birthday bound, Floyd and Brent cycle finding, batched gcds, a complete factoring driver and its failure modes

The idea: a birthday collision modulo a prime you cannot see

Pick a polynomial f(x) = x² + c and a start value x₀, and iterate x_i+1 = f(x_i) mod n. Now imagine the same walk reduced modulo an unknown prime factor p of n. Because reduction mod p commutes with squaring and adding, the residues x_i mod p follow the rule x → x² + c mod p on their own. That map lives in a set of only p values, so the walk must eventually repeat, and once one value repeats every later value repeats too. Drawn out, the walk is a tail leading into a loop, the shape of the Greek letter rho.

How soon does it repeat? If the map behaves like a random function, this is the birthday problem: among k random values from p possibilities a collision becomes likely once k is around √p. The expected tail-plus-cycle length is about √(πp/2). For p near 10⁹ that is about 40,000 steps, against a billion for trial division. The randomness is a heuristic, not a proof, but experiments agree with it.

The trick that turns a collision into a factor: if x_i ≡ x_j (mod p) then p divides x_i − x_j, and p also divides n. So gcd(|x_i − x_j|, n) is a multiple of p, computed entirely from values we can see. Unless the walk happens to collide modulo every prime factor at the same step, the gcd is a proper divisor of n. The Euclidean algorithm does that gcd in logarithmic time.

A hand trace on 8051

Take n = 8051, c = 1, x₀ = 2. The walk starts 2, 5, 26, 677, 7474, 2839, 871. Floyd's method keeps two pointers: a tortoise x_i that advances one step and a hare x_2i that advances two, and after each move it computes gcd(|x_i − x_2i|, n).

Three steps, and 8051 = 97 × 83. The diagram shows why. Modulo 97 the walk is 2, 5, 26, 95, 5, 26, 95: a tail of one value and a cycle of length 3. Floyd's pointers meet at the first i that is a multiple of the cycle length and at least the tail length, which here is 3. Modulo 83 the walk is 2, 5, 26, 13, 4, 17, 41, 22, 70, 4, with tail 4 and cycle 5, so it would only have been caught at step 5. Rho returns whichever factor's cycle closes first, and nothing about it favours the smallest prime.

Floyd's version in code

Floyd's version is the one to learn first. It costs three evaluations of f and one gcd per step.

from math import gcd

def rho_floyd(n, c=1, x0=2):
    """Return a divisor of n found by Pollard's rho; may return n itself on failure."""
    f = lambda v: (v * v + c) % n
    x = y = x0
    d = 1
    while d == 1:
        x = f(x)          # tortoise: one step
        y = f(f(y))       # hare: two steps
        d = gcd(abs(x - y), n)
    return d              # 1 < d < n on success, d == n when both cycles closed together

The function can only stop when the gcd exceeds 1, so it always terminates once the walk modulo every prime has cycled. If n is prime the walk modulo n itself cycles and the function returns n, so it never spins forever. It does waste about √n steps, which is why primality is checked first.

Brent's variant and batched gcds

Floyd pays for the hare's extra f evaluations and for a gcd every step, and a gcd costs far more than a modular multiply. Richard Brent's 1980 variant fixes both. It keeps a saved value x and walks y forward over windows of length r = 1, 2, 4, 8 and so on, comparing y with the saved x at each step. Because r doubles, some window is eventually at least as long as the cycle and starts after the tail, and then y meets x. Each step costs one f evaluation instead of three.

The second trick is to multiply the differences together and take one gcd per batch of m steps. If p divides any single difference it divides the product q, so gcd(q, n) still reveals p. The batch can overshoot: when two factors' cycles both close inside one batch, the gcd comes back as n. Brent's answer is to remember the y at the start of the batch and replay that batch one step at a time.

from math import gcd

def rho_brent(n, c, x0=2, m=128):
    """Brent's cycle finding with batched gcds. Returns a divisor of n, possibly n."""
    f = lambda v: (v * v + c) % n
    y, r, q, g = x0, 1, 1, 1
    while g == 1:
        x = y                      # saved point for this window
        for _ in range(r):
            y = f(y)
        k = 0
        while k < r and g == 1:
            ys = y                 # start of the batch, kept for backtracking
            for _ in range(min(m, r - k)):
                y = f(y)
                q = q * abs(x - y) % n
            g = gcd(q, n)
            k += m
        r *= 2
    if g == n:                     # the batch overshot: replay it step by step
        while True:
            ys = f(ys)
            g = gcd(abs(x - ys), n)
            if g > 1:
                break
    return g                       # still n means this c failed; caller retries

Measured with c = 1 and x₀ = 2, counting evaluations of f:

The last row matches the theory: √(π · 10⁹ / 2) is about 39,600. Batching barely changes the f calls but removes 127 of every 128 gcds, which is most of the run time in big-integer code. On 10403 the batch overshot to n and the replay recovered 101.

From one divisor to a full factorization

Rho finds one divisor. The driver around it strips small primes, proves pieces prime with Miller-Rabin, and recurses on whatever rho splits off.

import math, random

SMALL = [p for p in range(2, 1000) if all(p % q for q in range(2, math.isqrt(p) + 1))]
MR_BASES = (2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41)   # deterministic below 3.3e24

def is_prime(n):
    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 //= 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 find_divisor(m, rng, attempts=64):
    for _ in range(attempts):
        c = rng.randrange(1, m - 1)
        if c == m - 2:             # c = -2 gives a degenerate map; skip it
            continue
        d = rho_brent(m, c, x0=rng.randrange(m))
        if d != m:
            return d
    raise RuntimeError(f"rho gave up on {m}; hand it to ECM")

def factor(n, seed=0):
    rng, out = random.Random(seed), []
    for p in SMALL:
        while n % p == 0:
            out.append(p)
            n //= p
    stack = [n] if n > 1 else []
    while stack:
        m = stack.pop()
        if is_prime(m):
            out.append(m)
            continue
        r = math.isqrt(m)
        if r * r == m:
            stack += [r, r]
            continue
        d = find_divisor(m, rng)
        stack += [d, m // d]
    return sorted(out)

assert factor(8051) == [83, 97]
assert factor(600851475143) == [71, 839, 1471, 6857]
assert factor(2**64 + 1) == [274177, 67280421310721]

One run shows why the recursion matters. Brent with m = 128 on 600851475143 returned 59569, which is 71 × 839 and not prime, so every piece must go back through the primality gate. The Miller-Rabin bases above are proven deterministic for every n below 3.3 × 10²⁴. Beyond that the test is probabilistic, so add random bases or a proving step.

Failure modes

These are the failures seen in real implementations, roughly in order of how often they cause trouble.

• Calling rho on a prime. It returns n after about √n steps, which for a 60-bit prime is a billion iterations. Always test primality first.
• Treating the returned divisor as prime. As the 59569 example shows, it often is not. Recurse.
• gcd equals n. Both factors' cycles closed together, or the batch overshot. Backtrack the batch. If it still equals n, change c and x₀. Retrying with the same parameters repeats the same walk exactly.
• Bad constants. c = 0 gives x_i = x₀²ⁱ, and c = −2 is linked to Chebyshev polynomials. Both walks have structure that breaks the random-map assumption. Pick c at random and avoid those two.
• Even n and prime powers. Rho modulo 2 has a tiny cycle and gains nothing. Strip small primes by trial division. For n = 25 with c = 1 the run returned 25 itself, while c = 2 found 5. A cheap square check avoids the round trip.
• Overflow in fixed-width languages. In C, C++, Java or Rust, x * x overflows 64 bits once n passes 2³². Use a 128-bit product, or Montgomery multiplication as described in the modular exponentiation article. Overflow does not crash; it silently produces wrong residues.
• No time budget. With two 30-digit prime factors rho will not finish. Bound the iterations and escalate.

Where rho stops being the right tool

Rho's cost depends on the smallest factor, not on n. By the square-root estimate a 10-digit factor takes about 10⁵ steps, a 15-digit one about 4 × 10⁷ and a 20-digit one about 10¹⁰. Each extra two digits in p multiplies the work by ten, so about 20 digits is the practical limit.

The next tool is Lenstra's elliptic curve method (ECM), which also depends on the size of p but grows sub-exponentially, and takes over from about 20 digits up to 60 or more. When n is a product of two primes of similar size, as an RSA modulus is, all such methods lose. They are replaced by the quadratic sieve and the general number field sieve, whose cost depends on the size of n itself. That is why rho is no threat to a 2048-bit RSA key: its expected work would be on the order of 2⁵¹² steps.

Trade-offs

What to do next

1. Type in rho_floyd and reproduce the 8051 trace by hand, then confirm the cycle lengths 3 and 5 by reducing the sequence mod 97 and mod 83.
2. Implement rho_brent with batching and check it on 10403, where the backtrack path fires with m = 128.
3. Wrap it in the factor driver and run the three asserts, then a random test that multiplies factor(n) back together for a million random n below 2⁶².
4. If you work in a fixed-width language, add a 128-bit or Montgomery mulmod and test it against Python on 62-bit semiprimes.
5. If you meet factors beyond about 20 digits, use an existing ECM implementation such as GMP-ECM or PARI/GP rather than tuning rho further.