Pollard's Rho with Brent Cycle Detection

Intuition, Vocabulary and Foundations

Pollard’s rho searches for a factor without trying every divisor. It repeatedly applies a polynomial such as $f(x)=x^2+c\pmod n$ and watches the orbit. Because only finitely many residues exist, the walk eventually repeats and forms a tail followed by a cycle—the shape that gives the method its “rho” name.

Suppose the composite candidate has an unknown prime factor $p$. Although the implementation computes modulo $n$, the same orbit exists modulo $p$. By the birthday-paradox heuristic, two positions likely collide modulo $p$ after roughly $O(\sqrt p)$ steps. If $x\equiv y\pmod p$ but $x\not\equiv y\pmod n$, then $p$ divides $\lvert x-y\rvert$, and $1\lt\gcd(\lvert x-y\rvert,n)\lt n$ reveals a nontrivial factor.

Brent’s cycle-detection strategy compares blocks of orbit values and increases the checkpoint distance geometrically. Instead of computing a gcd for every difference, it multiplies several differences modulo $n$ and performs one gcd on the batch product. If the gcd becomes $n$, the batch contained too much information; replaying it one difference at a time can recover a factor.

The polynomial, starting value, and batch behavior affect runtime. An unlucky walk can collapse into an unhelpful cycle. Randomized restarts are therefore part of the algorithm, not merely a performance embellishment.

Pollard’s rho does not prove that its current walk will succeed; it relies on collision probability modulo an unknown factor. Correct implementations detect $g=n$, recover within the failed batch when possible, and restart with fresh polynomial parameters when necessary.

Concept	Definition or invariant	Concrete example	Non-example or boundary	Question answered
Polynomial orbit	$x_{k+1}=f(x_k)\bmod n$	$f(x)=x^2+c\bmod n$	A constant map creates useless cycles	How are pseudorandom residues generated?
Birthday collision	Repetition expected after about the square root of the state-space size	Collision modulo $p$ after roughly $\sqrt p$ steps	It is an expectation, not a deadline	Why can a factor appear early?
Factor signal	$g=\gcd(\lvert x-y\rvert,n)$	$1\lt g\lt n$ yields a factor	$g=1$ gives no factor; $g=n$ gives no separation	How does collision expose divisibility?
Brent checkpoint	Hold one value while advancing another through a power-of-two block	Block lengths $1,2,4,8,\ldots$	Advancing both identically keeps difference zero	How is cycle detection organized?
Gcd batching	$q=\prod\lvert x-y_j\rvert\bmod n$	One gcd checks many differences	Unreduced multiplication can overflow	How are gcd calls amortized?
Recovery phase	Replay a failed batch one difference at a time	Used after $\gcd(q,n)=n$	Returning $n$ makes no progress	How is an overfull batch repaired?
Restart	Choose a new seed and polynomial constant	Replace $(x_0,c)$ after failure	Repeating identical parameters repeats failure	How does the search escape bad walks?
Recursive splitting	Split composites until every leaf is prime	Split $n$ into $d$ and $n/d$	A rho divisor is not automatically prime	How is full factorization produced?

Step-by-Step Derivations and Algorithms

Begin by removing tiny prime factors and testing whether $n$ is prime. Pollard’s rho should receive an odd composite; otherwise it performs unnecessary work or may return the input itself.

For one Brent attempt, choose a starting value $y$, polynomial constant $c$, and batch size $m$. Set the block length to $r=1$. At each outer round, save $x=y$, advance $y$ by $r$ polynomial steps, then compare $x$ with successive values in the new block.

For each sub-block, maintain $q\equiv\prod_j\lvert x-y_j\rvert\pmod n$. Every prime factor dividing an accumulated difference also divides $q$, so $\gcd(q,n)$ detects a collision modulo a factor. If the gcd is $1$, continue. If it is proper, return it. If it equals $n$, restore the most recent sub-block start and advance one step at a time until a proper factor appears or recovery fails.

After finding $d$, split the value into $d$ and $n/d$. Test each piece for primality. Prime pieces enter the output; composite pieces return to rho. Sorting at the end gives a canonical factor list.

Expected work to find the smallest prime factor $p$ is $O(\sqrt p)$ polynomial steps under the random-map heuristic. For a balanced semiprime $n=pq$ with $p\approx q\approx\sqrt n$, this is approximately $O(n^{1/4})$. It is not a worst-case guarantee.

Pollard's Rho with Brent Cycle Detection algorithm map

flowchart TD
    A["Receive composite n"] --> B["Choose random seed and c"]
    B --> C["Advance through power-of-two block"]
    C --> D["Batch products of differences"]
    D --> E["Compute gcd with n"]
    E -->|g equals 1| F{"Block complete?"}
    F -->|No| D
    F -->|Yes| G["Double block length"]
    G --> C
    E -->|proper factor| H["Return factor d"]
    E -->|g equals n| I["Replay last batch individually"]
    I --> J{"Factor recovered?"}
    J -->|Yes| H
    J -->|No| K["Restart with new parameters"]
    K --> B
    H --> L["Split d and n divided by d"]
    L --> M{"Each piece prime?"}
    M -->|No| A
    M -->|Yes| N["Collect and sort factors"]

Within a Brent sub-block, $q$ is the product modulo $n$ of all tested differences, so any factor detected by an individual difference also divides $\gcd(q,n)$. Batching reduces gcd frequency, while expected search work remains about $O(\sqrt p)$ modular steps for the smallest prime factor $p$; recursive splitting adds primality tests and further searches.

Implementation, Complexity and Engineering

The reference implementation combines deterministic Miller–Rabin for unsigned 64-bit primality with Brent’s rho search. A seeded random generator supports reproducible tests; production callers should use securely seeded randomness when inputs or timing behavior are sensitive.

Python’s arbitrary-precision multiplication safely evaluates $x^2+c$ and batched products. Fixed-width implementations need unsigned 128-bit intermediates or another overflow-safe modular multiplication method. The factorization routine uses an explicit stack and sorts its final factors.

from math import gcd
from random import Random, SystemRandom


UINT64_MAX = (1 << 64) - 1
MR_BASES_64 = (2, 325, 9_375, 28_178, 450_775, 9_780_504, 1_795_265_022)
SMALL_PRIMES = (2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37)


def is_prime_u64(n):
    if not isinstance(n, int):
        raise TypeError("n must be an integer")
    if not 0 <= n <= UINT64_MAX:
        raise ValueError("n must be an unsigned 64-bit integer")
    if n < 2:
        return False
    for prime in SMALL_PRIMES:
        if n == prime:
            return True
        if n % prime == 0:
            return False

    d = n - 1
    s = 0
    while d % 2 == 0:
        s += 1
        d //= 2

    for base in MR_BASES_64:
        base %= n
        if base == 0:
            continue
        value = pow(base, d, n)
        if value == 1 or value == n - 1:
            continue
        for _ in range(s - 1):
            value = value * value % n
            if value == n - 1:
                break
            if value == 1:
                return False
        else:
            return False
    return True


def brent_factor(n, rng=None, batch_size=128, max_restarts=64):
    if not isinstance(n, int):
        raise TypeError("n must be an integer")
    if not 4 <= n <= UINT64_MAX:
        raise ValueError("n must be a composite unsigned 64-bit integer")
    if batch_size < 1 or max_restarts < 1:
        raise ValueError("batch size and restart limit must be positive")
    if n % 2 == 0:
        return 2
    if is_prime_u64(n):
        raise ValueError("n is prime")

    rng = rng if rng is not None else SystemRandom()

    for _ in range(max_restarts):
        y = rng.randrange(1, n)
        c = rng.randrange(1, n)
        block = 1
        factor = 1
        checkpoint = y
        saved_y = y

        def step(value):
            return (value * value + c) % n

        while factor == 1:
            checkpoint = y
            for _ in range(block):
                y = step(y)

            offset = 0
            while offset < block and factor == 1:
                saved_y = y
                product = 1
                width = min(batch_size, block - offset)
                for _ in range(width):
                    y = step(y)
                    product = product * abs(checkpoint - y) % n
                factor = gcd(product, n)
                offset += width
            block *= 2

        if factor == n:
            factor = 1
            y = saved_y
            for _ in range(batch_size):
                y = step(y)
                factor = gcd(abs(checkpoint - y), n)
                if factor > 1:
                    break

        if 1 < factor < n:
            return factor

    raise RuntimeError("rho search exhausted its restart budget")


def factor_u64(n, rng=None):
    if not isinstance(n, int):
        raise TypeError("n must be an integer")
    if not 1 <= n <= UINT64_MAX:
        raise ValueError("n must be a positive unsigned 64-bit integer")
    if n == 1:
        return []

    rng = rng if rng is not None else SystemRandom()
    factors = []
    pending = [n]
    while pending:
        value = pending.pop()
        if value == 1:
            continue
        if is_prime_u64(value):
            factors.append(value)
            continue
        divisor = brent_factor(value, rng=rng)
        pending.extend((divisor, value // divisor))
    factors.sort()
    return factors


def verify_factorization(n, factors):
    product = 1
    for factor in factors:
        if not is_prime_u64(factor):
            return False
        product *= factor
    return product == n


assert factor_u64(1, Random(1)) == []
assert factor_u64(2, Random(1)) == [2]
assert factor_u64(2**20, Random(2)) == [2] * 20
assert factor_u64(3**12, Random(3)) == [3] * 12

small_semiprime = 83 * 97
assert factor_u64(small_semiprime, Random(4)) == [83, 97]

hard_semiprime = 1_000_003 * 1_000_033
hard_factors = factor_u64(hard_semiprime, Random(5))
assert hard_factors == [1_000_003, 1_000_033]
assert verify_factorization(hard_semiprime, hard_factors)

mixed = 2**8 * 3**5 * 101 * 10_007
mixed_factors = factor_u64(mixed, Random(6))
assert verify_factorization(mixed, mixed_factors)
assert mixed_factors == sorted(mixed_factors)

try:
    brent_factor(101, Random(7))
    raise AssertionError("prime input was accepted")
except ValueError:
    pass

try:
    factor_u64(0, Random(8))
    raise AssertionError("zero was accepted")
except ValueError:
    pass

print("all tests passed")

Engineering choice	Effect	Expected cost	Failure mode	Mitigation
Random polynomial constant	Changes the orbit	Constant setup per restart	Fixed bad constant repeats failure	Resample constant and starting value
Brent detection	Reduces comparisons and gcd calls	Expected $O(\sqrt p)$ steps	Wrong checkpoint loses recovery data	Preserve the last sub-block start
Gcd batching	Amortizes gcd operations	One gcd per batch	Product yields $g=n$	Replay the failed batch individually
Batch size	Trades gcd count against recovery work	Workload-dependent	Huge batches overfill frequently	Benchmark values such as 64 or 128
Safe multiplication	Preserves modular residues	Backend-dependent	Native multiplication overflows	Use wider or reduced multiplication
Deterministic test RNG	Reproduces walks	Same algorithmic work	Shipped as security randomness	Inject RNG and use system randomness in production
Explicit work stack	Splits cofactors iteratively	Proportional to factor count	Deep recursion exhausts stack	Store pending cofactors explicitly
Sorted output	Canonicalizes results	$O(k\log k)$	Discovery order reaches callers	Sort and verify the product

Applications, Case Studies and Exercises

Pollard’s rho works well for machine-sized composites whose factors are too large for trial division but still fit efficient modular arithmetic. It does not replace primality testing: every returned divisor must be split until each output factor is certified prime.

Practice

Apply Pollard's Rho with Brent Cycle Detection

Case study 1 — trace a small semiprime. Factor $8051=83\cdot97$ with a fixed seed. Record $y_0$, $c$, checkpoints, batch products, and gcd results. Identify the first proper gcd and verify both cofactors are prime. Success means the factors reconstruct $8051$ exactly.

Case study 2 — factor a production-shaped semiprime. Factor $1000003\cdot1000033$ using several seeds. Record polynomial steps, gcd calls, recovery events, restarts, and elapsed time. Verify every factor with deterministic Miller–Rabin, sort the result, and multiply it back. Success means every run returns the same prime multiset despite different walks.

Hands-on restart experiment. Generate seeded semiprimes in several bit ranges. Run Brent rho with batch sizes $1,16,64,$ and $128$. Measure median steps, gcd calls, recovery frequency, and restarts. Include deliberately poor parameters in a diagnostic variant. Success means the report distinguishes stochastic runtime from correctness.

Tiered exercises. Recall: explain how a collision modulo $p$ yields a gcd factor. Derive: justify the $O(\sqrt p)$ birthday heuristic. Implement: count steps and gcd calls without changing the orbit. Debug: remove restart randomization and exhibit repeatable failure. Extend: add bounded trial division and measure its crossover point.

Use Pollard’s rho after removing tiny factors and certifying that the remaining value is composite; use Brent batching to reduce gcd overhead, randomized restarts to escape failed walks, and deterministic primality testing to certify every leaf. Never return $1$ or $n$ as progress, or trust factors without sorting, primality checks, and product reconstruction.

Cookie Consent