Skip to main content

Number Theory for CP

Why Number Theory Matters

Many competitive programming problems require mathematical insights. Modular arithmetic, prime factorization, and combinatorics appear in 30%+ of CF problems rated 1400+.
Pattern Recognition Signals:
  • “Output answer modulo 10^9+7” → Modular arithmetic
  • “Count ways” → Combinatorics (nCr, nPr)
  • “GCD/LCM” → Euclidean algorithm
  • “Divisors” → Prime factorization
  • Large exponents → Binary exponentiation

Modular Arithmetic

The Basics

When numbers get too large, we work modulo some prime (usually 10^9 + 7). Why modulo 10^9 + 7?
  • It’s a prime number (required for modular inverse via Fermat’s theorem)
  • It fits in a 32-bit integer
  • Two numbers modulo this can be multiplied without overflow in 64-bit
The Core Idea: Instead of computing huge numbers like 210000002^{1000000}, we keep everything “wrapped around” within [0,MOD1][0, MOD-1]. Think of it like a clock—after 12, we go back to 1. Key Properties (these let us modulo at every step):
  • (a+b)modm=((amodm)+(bmodm))modm(a + b) \mod m = ((a \mod m) + (b \mod m)) \mod m
  • (a×b)modm=((amodm)×(bmodm))modm(a \times b) \mod m = ((a \mod m) \times (b \mod m)) \mod m
  • (ab)modm=((amodm)(bmodm)+m)modm(a - b) \mod m = ((a \mod m) - (b \mod m) + m) \mod m (add mm to handle negatives)
const ll MOD = 1e9 + 7;

// Properties:
// (a + b) % m = ((a % m) + (b % m)) % m
// (a * b) % m = ((a % m) * (b % m)) % m
// (a - b) % m = ((a % m) - (b % m) + m) % m  // Add m to handle negative

ll add(ll a, ll b) { return ((a % MOD) + (b % MOD)) % MOD; }
ll sub(ll a, ll b) { return ((a % MOD) - (b % MOD) + MOD) % MOD; }
ll mul(ll a, ll b) { return ((a % MOD) * (b % MOD)) % MOD; }
Division is Different! (a / b) % m ≠ ((a % m) / (b % m)) % mFor division, multiply by modular inverse: a / b ≡ a × b^(-1) (mod m)

Binary Exponentiation

The Problem: Computing ana^n directly takes O(n) multiplications. For n=1018n = 10^{18}, that’s impossible. The Insight: Use the property that:
  • If n is even: an=(an/2)2a^n = (a^{n/2})^2
  • If n is odd: an=a×an1a^n = a \times a^{n-1}
Example: Computing 3133^{13}:
  • 313=3×3123^{13} = 3 \times 3^{12}
  • 312=(36)23^{12} = (3^6)^2
  • 36=(33)23^6 = (3^3)^2
  • 33=3×323^3 = 3 \times 3^2
  • 32=(31)23^2 = (3^1)^2
  • 31=3×303^1 = 3 \times 3^0
This reduces O(n) to O(log n) multiplications! Bit Interpretation: 13 in binary is 1101. We compute 31,32,34,383^1, 3^2, 3^4, 3^8 and multiply those where the bit is 1: 313=38×34×313^{13} = 3^8 \times 3^4 \times 3^1.
ll power(ll base, ll exp, ll mod = MOD) {
    ll result = 1;
    base %= mod;
    while (exp > 0) {
        if (exp & 1) result = result * base % mod;
        base = base * base % mod;
        exp >>= 1;
    }
    return result;
}

Modular Inverse

What is it? The modular inverse of aa modulo mm is a number a1a^{-1} such that: a×a11(modm)a \times a^{-1} \equiv 1 \pmod{m} Think of it as “division in modular world.” To compute xamodm\frac{x}{a} \mod m, we compute x×a1modmx \times a^{-1} \mod m. Fermat’s Little Theorem: For prime pp and aa not divisible by pp: ap11(modp)a^{p-1} \equiv 1 \pmod{p} Rearranging: a×ap21(modp)a \times a^{p-2} \equiv 1 \pmod{p}, so a1ap2(modp)a^{-1} \equiv a^{p-2} \pmod{p} When does inverse exist?
  • For prime modulus: Always exists if gcd(a,p)=1\gcd(a, p) = 1 (i.e., a≢0a \not\equiv 0)
  • For general modulus: Only if gcd(a,m)=1\gcd(a, m) = 1
ll modinv(ll a, ll mod = MOD) {
    return power(a, mod - 2, mod);
}

ll divide(ll a, ll b, ll mod = MOD) {
    return mul(a, modinv(b, mod));
}

Extended Euclidean Algorithm

The Problem: Find integers x,yx, y such that ax+by=gcd(a,b)ax + by = \gcd(a, b) Why it matters: This is Bézout’s Identity—such x,yx, y always exist! How it works:
  1. Base case: When b=0b = 0, we have gcd(a,0)=a\gcd(a, 0) = a, so x=1,y=0x = 1, y = 0 works: a(1)+0(0)=aa(1) + 0(0) = a
  2. Recursive case: If we know x,yx', y' for (b,amodb)(b, a \mod b) such that bx+(amodb)y=gbx' + (a \mod b)y' = g, then:
    • amodb=aa/bba \mod b = a - \lfloor a/b \rfloor \cdot b
    • Substituting: bx+(aa/bb)y=gbx' + (a - \lfloor a/b \rfloor \cdot b)y' = g
    • Rearranging: ay+b(xa/by)=gay' + b(x' - \lfloor a/b \rfloor \cdot y') = g
    • So x=yx = y' and y=xa/byy = x' - \lfloor a/b \rfloor \cdot y'
Application: Finding modular inverse when modulus is not prime. If gcd(a,m)=1\gcd(a, m) = 1, then from ax+my=1ax + my = 1, we get ax1(modm)ax \equiv 1 \pmod{m}, so xx is the inverse.
ll extgcd(ll a, ll b, ll &x, ll &y) {
    if (b == 0) {
        x = 1; y = 0;
        return a;
    }
    ll x1, y1;
    ll g = extgcd(b, a % b, x1, y1);
    x = y1;
    y = x1 - (a / b) * y1;
    return g;
}

// Modular inverse for any coprime a, m (not just prime m)
ll modinv_general(ll a, ll m) {
    ll x, y;
    ll g = extgcd(a, m, x, y);
    if (g != 1) return -1;  // No inverse exists
    return (x % m + m) % m;
}

Prime Numbers

Sieve of Eratosthenes

The Problem: Find all prime numbers up to N. The Insight: Instead of testing each number for primality (expensive), we “cross out” multiples of each prime. Algorithm:
  1. Start with all numbers marked as potentially prime
  2. The first unmarked number (2) is prime—cross out all its multiples (4, 6, 8, …)
  3. Next unmarked (3) is prime—cross out 6, 9, 12, …
  4. For number ii, start crossing from i2i^2 (smaller multiples already crossed by smaller primes)
  5. Repeat until i2>Ni^2 > N
Why O(N log log N)? Each prime pp crosses out N/pN/p numbers. The sum pNN/p=NpN1/p\sum_{p \leq N} N/p = N \sum_{p \leq N} 1/p is O(NloglogN)O(N \log \log N) by the prime harmonic series. Visual Example for N = 20:
Start:  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18 19 20
After 2: ✓  3  X  5  X  7  X  9  X 11  X 13  X 15  X 17  X 19  X
After 3: ✓  ✓  X  5  X  7  X  X  X 11  X 13  X  X  X 17  X 19  X
After 5: (5² = 25 > 20, so we stop)
Primes: 2, 3, 5, 7, 11, 13, 17, 19
const int MAXN = 1e7 + 5;
bool isPrime[MAXN];
vector<int> primes;

void sieve() {
    fill(isPrime, isPrime + MAXN, true);
    isPrime[0] = isPrime[1] = false;
    
    for (int i = 2; i < MAXN; i++) {
        if (isPrime[i]) {
            primes.push_back(i);
            for (ll j = (ll)i * i; j < MAXN; j += i) {
                isPrime[j] = false;
            }
        }
    }
}

Linear Sieve with SPF

The Improvement: Standard sieve can mark a composite multiple times (12 is marked by both 2 and 3). The linear sieve ensures each composite is marked exactly once—by its smallest prime factor (SPF). Key Invariant: When processing number ii:
  • We multiply ii by primes pspf[i]p \leq \text{spf}[i]
  • This ensures ipi \cdot p is marked by its SPF, which is pp
  • We stop when p>spf[i]p > \text{spf}[i] because then ipi \cdot p would have SPF = spf[i]\text{spf}[i], not pp
Why SPF is useful:
  • O(log n) factorization: repeatedly divide by SPF
  • Euler’s totient: compute multiplicatively using prime factorization
int spf[MAXN];  // Smallest Prime Factor

void linearSieve() {
    for (int i = 2; i < MAXN; i++) {
        if (spf[i] == 0) {
            spf[i] = i;
            primes.push_back(i);
        }
        for (int p : primes) {
            if (p > spf[i] || (ll)i * p >= MAXN) break;
            spf[i * p] = p;
        }
    }
}

// O(log n) factorization using SPF
vector<pair<int, int>> factorize(int n) {
    vector<pair<int, int>> factors;
    while (n > 1) {
        int p = spf[n], cnt = 0;
        while (n % p == 0) {
            n /= p;
            cnt++;
        }
        factors.push_back({p, cnt});
    }
    return factors;
}

Prime Factorization (O(√n))

vector<pair<ll, int>> factorize(ll n) {
    vector<pair<ll, int>> factors;
    for (ll p = 2; p * p <= n; p++) {
        if (n % p == 0) {
            int cnt = 0;
            while (n % p == 0) {
                n /= p;
                cnt++;
            }
            factors.push_back({p, cnt});
        }
    }
    if (n > 1) factors.push_back({n, 1});
    return factors;
}

Counting Divisors

If n = p1^a1 × p2^a2 × … × pk^ak, then:
  • Number of divisors = (a1+1) × (a2+1) × … × (ak+1)
  • Sum of divisors = ∏((pi^(ai+1) - 1) / (pi - 1))
ll countDivisors(ll n) {
    ll cnt = 1;
    for (ll p = 2; p * p <= n; p++) {
        if (n % p == 0) {
            int exp = 0;
            while (n % p == 0) {
                n /= p;
                exp++;
            }
            cnt *= (exp + 1);
        }
    }
    if (n > 1) cnt *= 2;
    return cnt;
}

// Get all divisors
vector<ll> getDivisors(ll n) {
    vector<ll> divs;
    for (ll i = 1; i * i <= n; i++) {
        if (n % i == 0) {
            divs.push_back(i);
            if (i != n / i) divs.push_back(n / i);
        }
    }
    sort(divs.begin(), divs.end());
    return divs;
}

Combinatorics

Factorial and Inverse Factorial

The Problem: Computing (nr)=n!r!(nr)!\binom{n}{r} = \frac{n!}{r!(n-r)!} for many queries. Naive approach: Compute each query in O(n)—too slow for 10510^5 queries. Smart approach: Precompute n!n! and (n!)1(n!)^{-1} for all nn. The Trick for Inverse Factorials: Instead of computing nn modular inverses (each costs O(log MOD)):
  1. Compute (MAXN1)!1(MAXN-1)!^{-1} using Fermat’s theorem
  2. Use the recurrence: (n!)1=(n+1)!(n+1)1(n!)^{-1} = (n+1)! \cdot (n+1)^{-1}… but even simpler:
    • (n!)1=((n+1)!)1(n+1)(n!)^{-1} = ((n+1)!)^{-1} \cdot (n+1)
    • This is because (n+1)!=n!(n+1)(n+1)! = n! \cdot (n+1), so (n!)1=((n+1)!)1(n+1)(n!)^{-1} = ((n+1)!)^{-1} \cdot (n+1)
This gives O(N) precomputation and O(1) per query.
const int MAXN = 2e5 + 5;
ll fact[MAXN], inv_fact[MAXN];

void precompute() {
    fact[0] = 1;
    for (int i = 1; i < MAXN; i++) {
        fact[i] = fact[i-1] * i % MOD;
    }
    
    // Only ONE modular inverse computation!
    inv_fact[MAXN-1] = power(fact[MAXN-1], MOD - 2);
    // The rest computed in O(1) each
    for (int i = MAXN - 2; i >= 0; i--) {
        inv_fact[i] = inv_fact[i+1] * (i+1) % MOD;
    }
}

ll C(int n, int r) {
    if (r < 0 || r > n) return 0;
    return fact[n] * inv_fact[r] % MOD * inv_fact[n-r] % MOD;
}

ll P(int n, int r) {
    if (r < 0 || r > n) return 0;
    return fact[n] * inv_fact[n-r] % MOD;
}

Common Formulas

FormulaExpressionUse Case
nCrn! / (r! × (n-r)!)Choose r from n
nPrn! / (n-r)!Arrange r from n
Stars and BarsC(n+r-1, r-1)Distribute n identical items into r groups
Catalan(n)C(2n, n) / (n+1)Valid parentheses, BST count
Derangementsn! × Σ((-1)^i / i!)Permutations with no fixed points

Pascal’s Triangle (Small n, r)

const int MAXN = 5005;
ll C[MAXN][MAXN];

void buildPascal() {
    for (int i = 0; i < MAXN; i++) {
        C[i][0] = 1;
        for (int j = 1; j <= i; j++) {
            C[i][j] = (C[i-1][j-1] + C[i-1][j]) % MOD;
        }
    }
}

Lucas Theorem (Large n, small p)

The Problem: Computing (nr)modp\binom{n}{r} \mod p when nn is huge (like 101810^{18}) but pp is a small prime. Why standard approach fails: We can’t precompute factorials up to 101810^{18}. Lucas Theorem: (nr)i=0k(niri)(modp)\binom{n}{r} \equiv \prod_{i=0}^{k} \binom{n_i}{r_i} \pmod{p} where n=nkpk+...+n1p+n0n = n_k p^k + ... + n_1 p + n_0 and r=rkpk+...+r1p+r0r = r_k p^k + ... + r_1 p + r_0 are base-pp representations. Intuition: Break down n and r into their “digits” in base p. Each digit contributes a small binomial coefficient (computable since values are < p). Example: (125)mod3\binom{12}{5} \mod 3
  • 12 in base 3: 110 → digits are (1, 1, 0)
  • 5 in base 3: 012 → digits are (0, 1, 2)
  • (125)(10)(11)(02)=110=0(mod3)\binom{12}{5} \equiv \binom{1}{0} \cdot \binom{1}{1} \cdot \binom{0}{2} = 1 \cdot 1 \cdot 0 = 0 \pmod{3}
ll lucas(ll n, ll r, ll p) {
    if (r == 0) return 1;
    return C(n % p, r % p) * lucas(n / p, r / p, p) % p;
}

GCD and LCM

Basic Functions

ll gcd(ll a, ll b) { return b ? gcd(b, a % b) : a; }
ll lcm(ll a, ll b) { return a / gcd(a, b) * b; }  // Prevent overflow

GCD of Array

ll gcdArray(vector<ll>& arr) {
    ll g = 0;
    for (ll x : arr) g = gcd(g, x);
    return g;
}

Properties

  • gcd(a, b) × lcm(a, b) = a × b
  • gcd(a, b, c) = gcd(gcd(a, b), c)
  • gcd(a, 0) = a
  • If gcd(a, b) = g, then gcd(a/g, b/g) = 1

Euler’s Totient Function

Definition: ϕ(n)\phi(n) = count of integers in [1,n][1, n] that are coprime to nn (i.e., gcd(k,n)=1\gcd(k, n) = 1). Examples:
  • ϕ(1)=1\phi(1) = 1 (just 1 itself)
  • ϕ(6)=2\phi(6) = 2 (1 and 5 are coprime to 6)
  • ϕ(7)=6\phi(7) = 6 (for prime p: all of 1, 2, …, p-1 are coprime)
  • ϕ(p)=p1\phi(p) = p - 1 for any prime pp
Formula: If n=p1a1×p2a2×...×pkakn = p_1^{a_1} \times p_2^{a_2} \times ... \times p_k^{a_k}, then: ϕ(n)=n×i=1k(11pi)=n×p11p1×p21p2×...×pk1pk\phi(n) = n \times \prod_{i=1}^{k} \left(1 - \frac{1}{p_i}\right) = n \times \frac{p_1 - 1}{p_1} \times \frac{p_2 - 1}{p_2} \times ... \times \frac{p_k - 1}{p_k} Intuition: Start with nn numbers. For each prime factor pp, exactly 1/p1/p of numbers are divisible by pp, so we keep (11/p)(1 - 1/p) of them. Example: ϕ(12)=ϕ(22×3)=12×(11/2)×(11/3)=12×1/2×2/3=4\phi(12) = \phi(2^2 \times 3) = 12 \times (1 - 1/2) \times (1 - 1/3) = 12 \times 1/2 \times 2/3 = 4 Indeed, the numbers coprime to 12 in [1,12] are: 1, 5, 7, 11 → 4 numbers ✓
ll phi(ll n) {
    ll result = n;
    for (ll p = 2; p * p <= n; p++) {
        if (n % p == 0) {
            while (n % p == 0) n /= p;
            result -= result / p;
        }
    }
    if (n > 1) result -= result / n;
    return result;
}

// Sieve for phi
int phi_sieve[MAXN];

void computePhi() {
    iota(phi_sieve, phi_sieve + MAXN, 0);
    for (int i = 2; i < MAXN; i++) {
        if (phi_sieve[i] == i) {  // i is prime
            for (int j = i; j < MAXN; j += i) {
                phi_sieve[j] -= phi_sieve[j] / i;
            }
        }
    }
}

Euler’s Theorem

If gcd(a, n) = 1, then a^φ(n) ≡ 1 (mod n) Application: a^b mod m when b is huge
  • If gcd(a, m) = 1: a^b ≡ a^(b mod φ(m)) (mod m)

Pattern: Counting with Inclusion-Exclusion

The Problem: Count elements satisfying at least one of k conditions. The Principle: A1A2...Ak=AiAiAj+AiAjAk...|A_1 \cup A_2 \cup ... \cup A_k| = \sum |A_i| - \sum |A_i \cap A_j| + \sum |A_i \cap A_j \cap A_k| - ... Alternately: Add single sets, subtract pairs, add triples, subtract quadruples, … Why it works (for 2 sets):
  • AB|A \cup B| counts elements in either A or B
  • A+B|A| + |B| overcounts elements in both by 1
  • So AB=A+BAB|A \cup B| = |A| + |B| - |A \cap B|
Bitmask Implementation: With k conditions, use a k-bit mask where bit i = 1 means “include condition i.” For each non-empty subset:
  • Count elements satisfying ALL conditions in the subset
  • Add if odd number of conditions, subtract if even
Classic Application: Count numbers in [1,n][1, n] divisible by at least one of the given primes.
  • Divisible by pp: there are n/p\lfloor n/p \rfloor such numbers
  • Divisible by both pp and qq: there are n/(pq)\lfloor n/(p \cdot q) \rfloor such numbers
// Count numbers in [1, n] divisible by at least one prime in list
ll countDivisible(ll n, vector<ll>& primes) {
    int k = primes.size();
    ll result = 0;
    
    for (int mask = 1; mask < (1 << k); mask++) {
        ll product = 1;
        int bits = 0;
        
        for (int i = 0; i < k; i++) {
            if (mask & (1 << i)) {
                product *= primes[i];
                bits++;
                if (product > n) break;
            }
        }
        
        if (product <= n) {
            if (bits % 2 == 1) result += n / product;
            else result -= n / product;
        }
    }
    
    return result;
}

Practice Problems

Beginner (1000-1300)

ProblemTopicLink
ExponentiationBinary expCSES
Counting DivisorsFactorizationCSES
GCD QueryGCDCF 1350C

Intermediate (1300-1600)

ProblemTopicLink
Binomial CoefficientsnCrCSES
Creating Strings IICombinatoricsCSES
Sum of DivisorsNumber theoryCSES

Advanced (1600-1900)

ProblemTopicLink
Counting Coprime PairsEuler phiCSES
Christmas PartyDerangementsCSES

Key Takeaways

Modular Inverse

Use a^(p-2) for prime p, or extended GCD for general.

Precompute Factorials

O(n) precomputation gives O(1) nCr queries.

SPF Sieve

Linear sieve enables O(log n) factorization.

Inclusion-Exclusion

Count by adding/subtracting overlaps.

Next Up

Chapter 24: Advanced Graph Algorithms

Master SCC, bridges, articulation points, and advanced graph techniques.