Primes, fast power & combinatorics

Sieve of Eratosthenes, binary exponentiation, modular inverse, and nCr mod p. · 12 min read

Sieve of Eratosthenes

To find all primes up to n, start with every number marked prime, then for each prime p cross out its multiples p², p²+p, …. What remains unmarked is prime. It runs in O(n log log n).

vector<bool> sieve(int n) {
    vector<bool> is_prime(n + 1, true);
    is_prime[0] = is_prime[1] = false;
    for (int p = 2; (long long)p * p <= n; p++)
        if (is_prime[p])
            for (int q = p * p; q <= n; q += p)
                is_prime[q] = false;
    return is_prime;
}

Interactive visualization: Sieve of Eratosthenes — cross out each prime's multiples — step through it frame by frame on this page.

Practice: Count Primes (Medium) — How many primes are strictly less than n? — the sieve.

Binary (fast) exponentiation

Computing base^exp by multiplying exp times is O(exp). Square the base and halve the exponent instead: if the current exponent bit is 1, fold the base into the result. That's O(log exp) multiplications — essential for big exponents under a modulus.

Exponentiation by squaring

base^exp: square base, halve exp — exp in binary drives the multiplications.
multiply into result when a bit is 1 — At most log₂ exp folds.
take mod at every multiply — Keeps values bounded.

Result: O(log exp)

long long mod_pow(long long b, long long e, long long m) {
    long long r = 1 % m; b %= m;
    while (e) {
        if (e & 1) r = r * b % m;
        b = b * b % m;
        e >>= 1;
    }
    return r;
}

Modular inverse

You can't divide under a modulus directly, but you can multiply by the inverse. When m is prime, Fermat's little theorem gives a^(m-1) ≡ 1 (mod m), so a⁻¹ ≡ a^(m-2) (mod m) — compute it with mod_pow(a, m-2, m).

Combinatorics mod a prime

C(n, r) = n! / (r!·(n-r)!). Under a prime modulus, compute the numerator and denominator mod p, then multiply by the modular inverse of the denominator. Precomputing factorials and inverse factorials answers many queries in O(1) each.

Extended Euclid & the general inverse

The extended Euclidean algorithm finds integers x, y with a·x + b·y = gcd(a, b). When gcd(a, m) = 1, that x (mod m) is the modular inverse of a — and unlike Fermat's trick, it works for any modulus, not just primes.

// returns gcd; sets x, y so that a*x + b*y = gcd(a, b)
long long ext_gcd(long long a, long long b, long long& x, long long& y) {
    if (b == 0) { x = 1; y = 0; return a; }
    long long x1, y1;
    long long g = ext_gcd(b, a % b, x1, y1);
    x = y1; y = x1 - (a / b) * y1;
    return g;
}

More number theory to know

Euler's totient φ(n) counts integers in 1..n coprime to n; a^φ(n) ≡ 1 (mod n) generalises Fermat.
Chinese Remainder Theorem (CRT) reconstructs a value from its residues modulo pairwise-coprime moduli.
Factorisation & divisors — trial division to √n, or factor fast with a sieve of smallest prime factors; the divisor-count/sum are multiplicative over the prime factorisation.

Number-theory resources

cp-algorithms — Number Theory: dense, reliable write-ups of the sieve, fast power, extended Euclid, modular inverse, CRT, Euler's totient, and more, with code.
Forthright48 — CPPS 101: a number-theory-heavy CP roadmap with worked blog posts (also linked from the CP intro).

Practice: Modular Exponentiation (Medium) — Compute base^exp mod m fast with binary exponentiation.

Practice: Binomial Coefficient mod p (Hard) — C(n, r) mod a prime, using factorials and a modular inverse.