Prime Factorisation — 10^9 Calls Killed API

RSA's security depends on one fact: factoring the product of two large primes is computationally hard. The best known classical algorithm (General Number Field Sieve) runs in sub-exponential but super-polynomial time. A 2048-bit RSA modulus would take longer than the age of the universe to factor on current hardware.

For small numbers — up to ~10^12 — trial division works fine. For competitive programming — up to ~10^6 per query — the SPF sieve turns factorization from O(√n) per query to O(log n). Knowing which regime you're in and choosing the right tool is the practical skill.

The real trap? Assuming one algorithm works for all sizes. Engineers waste hours because they don't detect when trial division becomes impractical or when Pollard's rho loops indefinitely on a prime input.

Trial Division

Check divisibility by all primes up to √n. O(√n) per number. This is the simplest and most intuitive algorithm. For a single n ≤ 10^12, trial division completes in microseconds. But scale it to 10^6 numbers and you're looking at days.

The key optimization: after checking 2, only check odd numbers. That halves the work. For even better performance, precompute a list of primes up to √max_n using a simple sieve, then only divide by those primes.

def prime_factors(n: int) -> dict[int, int]:
    """Returns {prime: exponent} factorization."""
    factors = {}
    d = 2
    while d * d <= n:
        while n % d == 0:
            factors[d] = factors.get(d, 0) + 1
            n //= d
        d += 1
    if n > 1:
        factors[n] = factors.get(n, 0) + 1
    return factors

print(prime_factors(360))   # {2:3, 3:2, 5:1}
print(prime_factors(1001))  # {7:1, 11:1, 13:1}
print(prime_factors(97))    # {97:1} — prime

Number of Divisors and Sum of Divisors

These formulas are direct consequences of the multiplicative property. For example, τ(360) = τ(2³)τ(3²)τ(5¹) = (3+1)(2+1)(1+1) = 4×3×2 = 24 divisors. The sum σ(360) = (2^(3+1)-1)/(2-1) × (3^(2+1)-1)/(3-1) × (5^(1+1)-1)/(5-1) = (16-1)/1 × (27-1)/2 × (25-1)/4 = 15 × 13 × 6 = 1170.

def num_divisors(n: int) -> int:
    factors = prime_factors(n)
    result = 1
    for exp in factors.values():
        result *= (exp + 1)
    return result

def sum_divisors(n: int) -> int:
    factors = prime_factors(n)
    result = 1
    for p, exp in factors.items():
        result *= (p**(exp+1) - 1) // (p - 1)
    return result

def euler_totient(n: int) -> int:
    factors = prime_factors(n)
    result = n
    for p in factors:
        result -= result // p
    return result

print(num_divisors(360))   # 24
print(sum_divisors(360))   # 1170
print(euler_totient(36))   # 12

Listing All Divisors

To list all divisors, iterate up to √n and collect pairs. Complexity O(√n) – but you can't beat that because output size is O(τ(n)) which is at most O(√n). For n=10^12, τ(n) is at most ~6720, so enumerating divisors is feasible.

The algorithm: for each divisor d ≤ √n, if n % d == 0, add d and n//d. Sort the result at the end.

def all_divisors(n: int) -> list[int]:
    divisors = []
    for i in range(1, int(n**0.5) + 1):
        if n % i == 0:
            divisors.append(i)
            if i != n // i:
                divisors.append(n // i)
    return sorted(divisors)

print(all_divisors(36))  # [1,2,3,4,6,9,12,18,36]
print(len(all_divisors(360)))  # 24

SPF Sieve for Batch Factorization

When factorizing many numbers up to n, precompute the Smallest Prime Factor (SPF) sieve. Each factorization then takes O(log n) instead of O(√n).

Construction: Create an array spf[0..n] filled with 0. For i from 2 to n, if spf[i]==0 (i is prime), set spf[i]=i and for j in range(i*i, n+1, i), if spf[j]==0, set spf[j]=i. Then to factor a number x, repeatedly divide by spf[x] and record the count.

This is the standard for competitive programming and any scenario where you need to factorise a batch of numbers up to ~10^7.

def spf_sieve(n: int) -> list[int]:
    spf = [0] * (n + 1)
    for i in range(2, n + 1):
        if spf[i] == 0:
            spf[i] = i
            if i * i <= n:
                for j in range(i * i, n + 1, i):
                    if spf[j] == 0:
                        spf[j] = i
    return spf

def factorize_with_spf(x: int, spf: list[int]) -> dict[int, int]:
    factors = {}
    while x > 1:
        p = spf[x]
        factors[p] = factors.get(p, 0) + 1
        x //= p
    return factors

spf = spf_sieve(100)
print(factorize_with_spf(72, spf))  # {2:3, 3:2}

Pollard's Rho Algorithm for Large Composite Numbers

When n exceeds ~10^12 and you need to factorise a single number, trial division becomes impractical. Pollard's rho is a probabilistic algorithm that finds a non-trivial factor in O(n^(1/4)) expected time.

It works by iterating a pseudo-random function f(x) = (x^2 + c) mod n and detecting cycles using Floyd's cycle detection. When a cycle is found, gcd(|x-y|, n) gives a factor.

It's not guaranteed – for prime n it will loop forever. Always pair it with a primality test (Miller-Rabin) before running rho.

import math

def pollard_rho(n: int) -> int:
    if n % 2 == 0:
        return 2
    c = 1
    x = 2
    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:
            c += 1
            x = 2
            y = 2
            d = 1
    return d

n = 1000000000000000003  # a large semiprime
print(pollard_rho(n))

Frequently Asked Questions