Number Theory & Combinatorics

Number Theory & Combinatorics | Algoryva

Intuition

How many ways? (and the answer will not fit in any variable)

How many ways can you pick 3 friends from 10 for a trip? 120. How many ways can you pick 50,000 from 100,000? A number with about 30,000 digits. Contest problems ask for such counts modulo 10⁹ + 7: you are only asked for the remainder, which fits comfortably in a variable.

The catch: counting formulas contain division (n! divided by r! and (n − r)!), and division does not survive taking remainders the way addition and multiplication do. This lesson fixes that with modular inverses, then builds the standard toolkit of counting.

First, the naive way

Take remainders, then divide

The formula is C(n, r) = n! / (r! · (n − r)!). Since the numbers are huge, a natural idea is to reduce each factorial modulo p first and then divide.

That is wrong. Remainders respect +, − and ×, but not /. For n = 2000, r = 1000 and the prime p = 1,000,003, the true answer is 311,237, while “reduce, then divide” gives 0. Computing exact big integers with math.comb works for one query, but a factorial of 10⁵ has about 456,000 digits, so thousands of queries are far too slow.

python

from math import comb, factorial

p = 1000003                      # a prime
n, r = 2000, 1000
exact = comb(n, r) % p
tempting = (factorial(n) % p) // ((factorial(r) % p) * (factorial(n - r) % p) % p)
print(exact, tempting)   # 311237 0
print(exact == tempting)   # False

Cost: Wrong as written, and exact big integers cost far too much for many queries.

The observation

Dividing by b means multiplying by b's inverse

Modulo a prime p, every non-zero number b has a modular inverse b⁻¹ with b · b⁻¹ ≡ 1. Then “divide by b” becomes “multiply by b⁻¹”, which does respect remainders. By Fermat's little theorem, b⁻¹ ≡ b^(p−2) (mod p), computable in O(log p) with fast exponentiation.

For many queries compute the factorials once, take one inverse (of n!), and walk down to get the inverse of every smaller factorial for free using inv[i−1] = inv[i] · i. Then each C(n, r) is two multiplications.

The concept

The toolkit

Number theory gives you the arithmetic; combinatorics gives you the counting rules.

Modular inverse: a⁻¹ mod m exists exactly when gcd(a, m) = 1. For a prime modulus use Fermat (pow(a, p − 2, p)); for any modulus use extended Euclid or pow(a, -1, m) (Python 3.8+).
Extended Euclid: finds x, y with a·x + b·y = gcd(a, b); when the gcd is 1, x mod m is the inverse of a modulo m.
Euler's totient φ(n): how many numbers from 1 to n are coprime to n. Euler's theorem: a^φ(m) ≡ 1 (mod m) when gcd(a, m) = 1.
Product and sum rules: independent choices multiply; alternatives add.
Permutations and combinations: arrangements of r out of n are n!/(n−r)!; choices (order ignored) are C(n, r) = n!/(r!(n−r)!). Arrangements with repeated items: n! / (k₁! k₂! …).
Pascal's rule: C(n, r) = C(n−1, r−1) + C(n−1, r). It needs no division, so it works for any modulus, in O(n²).
Stars and bars: distributing n identical items into k boxes (boxes may be empty) can be done in C(n + k − 1, k − 1) ways.
Inclusion–exclusion: to count elements with *at least one* of several properties, add singles, subtract pairs, add triples, and so on.
Catalan numbers: C(2n, n)/(n + 1) counts balanced bracket strings of length 2n (and many similar structures): 1, 1, 2, 5, 14, 42, …

See it move

Pascal's triangle

Cell (n, r) is C(n, r). Every inner cell is the sum of the two cells above it (the two yellow cells): either the first item is chosen, or it is not. Notice the symmetry C(n, r) = C(n, n − r) and that each row adds up to 2ⁿ, the number of subsets of n items.

Worked by hand

C(5, 2) three ways

Choosing 2 items out of 5, by the formula, by Pascal's rule, and with inverses modulo 7.

method	working	result
formula	5! / (2! · 3!) = 120 / (2 · 6)	10
Pascal's rule	C(4, 1) + C(4, 2) = 4 + 6	10
modulo 7, with inverses	5! = 120 ≡ 1; 2! = 2 has inverse 4; 3! = 6 has inverse 6; 1 · 4 · 6 = 24 ≡ 3	3 (and 10 mod 7 = 3)

1
Why the inverses are 4 and 6
Modulo 7: 2 · 4 = 8 ≡ 1, and 6 · 6 = 36 ≡ 1. So dividing by 2 means multiplying by 4, and dividing by 6 means multiplying by 6.
2
A word of caution about small primes
This is only valid because 5 < 7. If n ≥ p, the factorial n! is divisible by p and has no inverse: the method breaks (see edge cases).

Implementation

Modular inverses

Extended Euclid returns the gcd together with the coefficients of a combination that produces it.

Plan in plain steps

pseudocode

ext_gcd(a, b):
    if b == 0: return (a, 1, 0)
    (g, x, y) = ext_gcd(b, a mod b)
    return (g, y, x − (a // b) · y)

python

def ext_gcd(a, b):
    if b == 0:
        return a, 1, 0
    g, x, y = ext_gcd(b, a % b)
    return g, y, x - (a // b) * y          # a*x + b*y = g

def mod_inverse(a, m):
    g, x, _ = ext_gcd(a % m, m)
    return x % m if g ==          

((, ))   
((, ), (, ), (, ))   
((, -, ), (,  - , ))   
((a * (a, ) %  ==   a  (, )))

Implementation

C(n, r) mod p in O(1) per query

Precompute factorials and inverse factorials once, then every query is two multiplications.

python

from math import comb

MOD = 10**9 + 7

def build(n):
    fact = [1] * (n + 1)
    for i in range(1, n + 1):
        fact[i] = fact[i - 1] * i % MOD
    inv = [1] * (n + 1)
    inv[n] = pow(fact[n], MOD - 2, MOD)          # one modular inverse (Fermat) ...
    for i in range(n, 0, -1):
        inv[i - 1] = inv[i] * i % MOD            # ... then every other one for free
    return fact, inv

fact, inv = (**)

 (n, r):
     r <   r > n:
         
     fact[n] * inv[r] % MOD * inv[n - r] % MOD

((, ), (, ), (, ), (, ))   
((, ))   
(((n, r) == (n, r) % MOD  n  (, )  r  (, n + )))   
((, ) == (, ) % MOD)

Implementation

Pascal's triangle for any modulus

If the modulus is not prime (or n is small and you want no inverses), build the triangle directly with Pascal's rule: only additions.

python

def pascal(n, m):
    t = [[0] * (n + 1) for _ in range(n + 1)]
    for i in range(n + 1):
        t[i][0] = 1
        for j in range(1, i + 1):
            t[i][j] = (t[i - 1][j - 1] + t[i - 1][j]) % m
    return t

t = pascal(6, 10**9 + 7)
print(t[6][:7])   # [1, 6, 15, 20, 15, 6, 1]
print(sum(t[6]) == 2**6)   # True
((, )[][])   
((, ** + )[][])

Implementation

Euler's totient with a sieve

φ(n) counts the numbers from 1 to n that are coprime to n. A sieve computes φ for every number up to n: for each prime p, multiply every multiple of p by (1 − 1/p), which is phi[j] -= phi[j] // p.

python

from math import gcd

def totients(n):
    phi = list(range(n + 1))
    for i in range(2, n + 1):
        if phi[i] == i:                        # i has not been touched, so it is prime
            for j in range(i, n + 1, i):
                phi[j] -= phi[j] // i
    return phi

phi = totients(30)
print(phi[1:13])   # [1, 1, 2, 2, 4, 2, 6, 4, 6, 4, 10, 4]
print(phi[30], phi[29])   # 8 28
print(all(phi[n] == (  k  (, n + )  (k, n) == )  n  (, )))   
((, phi[], ))

Implementation

Counting rules, checked by brute force

Stars and bars, Catalan numbers and arrangements with repeats, each checked against direct enumeration on small inputs.

python

from itertools import product
from math import comb, factorial

def stars_and_bars(n, k):                      # n identical items into k boxes
    return comb(n + k - 1, k - 1)

def stars_brute(n, k):
    return sum(1 for xs in product(range(n + 1), repeat=k) if sum(xs) == n)

def catalan(n):
    return ( * n, n) // (n + )

 (n):                               
    count = 
     s  (, repeat= * n):
        depth = 
         ch  s:
            depth +=   ch ==   -
             depth < :
                
        :
            count += depth == 
     count

((, ), (, ))   
([(i)  i  ()])   
(((n) == (n)  n  (, )))   
(() // (() * () * ()))

Complexity

Costs

The theme: pay a one-time cost (tables, inverses, factorisation) so that each query is nearly free.

Fast modular exponentiation (one inverse by Fermat)

time O(log p)space O(1)

Square-and-multiply on the exponent p − 2.

Extended Euclid / pow(a, -1, m)

time O(log m)space O(log m) recursion

Same number of steps as ordinary Euclid.

Factorial and inverse-factorial tables up to N

time O(N + log p)space O(N)

Two linear passes and a single modular exponentiation.

One query C(n, r) mod p after the tables

time O(1)space O(1)

Two multiplications and two reductions.

Pascal's triangle up to n

time O(n²)space O(n²)

One addition per cell; works for any modulus.

Totients for 1..n with a sieve

time O(n log log n)space O(n)

Each prime updates its multiples, like the sieve of Eratosthenes.

Inclusion–exclusion over k distinct primes

time O(2^k)space O(k)

Edge cases

r < 0 or r > nC(n, r) = 0. Return 0 explicitly rather than indexing a table out of range.
r = 0 or r = nC(n, r) = 1. The table formula handles it because fact[0] = inv[0] = 1.
n ≥ p (n at least the prime modulus)n! is divisible by p, so it has no inverse and the factorial method fails. Use Lucas's theorem, or keep n below p (10⁹ + 7 is bigger than any table you can build).
The modulus is not primeFermat inverses do not apply. Use extended Euclid (only works if gcd = 1), Pascal's triangle, or factor the modulus and combine results.
Table too smallBuild tables up to the largest n in any query, not an arbitrary 10⁵.
Negative valuesPython's % gives a non-negative result; in C++ use ((x % m) + m) % m after a subtraction.
Overflow in other languagesProducts of two residues below 10⁹ + 7 need 64-bit integers; reduce after every multiplication.

Common mistakes

Mistakes to avoid

Reducing factorials mod p and then using / or //.
Multiply by the modular inverse instead.
Computing one pow(..., p - 2, p) per query for every factorial.
Take one inverse for the largest factorial and derive the rest in one backward pass.
Forgetting % MOD after each multiplication.
Reduce after every product; in other languages an unreduced product overflows.
Using Fermat's inverse with a composite modulus.
It only works for a prime modulus. Use extended Euclid or Pascal's triangle otherwise.
Treating the 0 case as having an inverse.
No inverse exists for a multiple of the modulus; check before inverting.

Recognize it in a contest

How to spot a counting problem

Reach for it when the problem says or implies

“Count the ways … modulo 10⁹ + 7” (the answer is too large to print exactly).
“Choose k of n”, “the number of subsets of size k”, “committees”, “teams”: C(n, k) with tables.
“Number of paths on a grid moving only right and down”: C(R + C, R), possibly split at a required cell.
“Distribute n identical items among k people”: stars and bars.
“Arrange letters, some of them equal”: multinomial n!/(k₁!k₂!…).
“Count numbers up to n that are (not) divisible by any of these primes”: inclusion–exclusion, or totients for coprimality.

Worked problem

Paths through a required cell

The problem

A robot starts at the top-left cell (1, 1) of an R × C grid and may move only right or down. How many different routes reach the bottom-right cell (R, C) while passing through the cell (a, b)? Print the count modulo 10⁹ + 7.

A route from (1, 1) to (a, b) has a − 1 down moves and b − 1 right moves in some order: C(a + b − 2, a − 1) ways.
From (a, b) to (R, C) it is C((R − a) + (C − b), R − a).
Passing through (a, b) means one route of each kind, so multiply the two counts (the product rule).
Check the formula with an independent table-based count on small grids.

python

MOD = 10** + 

 (n):
    fact = [] * (n + )
     i  (, n + ):
        fact[i] = fact[i - ] * i % MOD
    inv = [] * (n + )
    inv[n] = (fact[n], MOD - , MOD)          
     i  (n, , -):
        inv[i - ] = inv[i] * i % MOD            
     fact, inv

fact, inv = (**)

 (n, r):
     r <   r > n:
         
     fact[n] * inv[r] % MOD * inv[n - r] % MOD

 (R, C, a, b):                 
     (a + b - , a - ) * (R - a + C - b, R - a) % MOD

 (R, C, a, b):                      
     (r1, c1, r2, c2):
         r2 < r1  c2 < c1:
             
        dp = [[] * (c2 - c1 + )  _  (r2 - r1 + )]
         i  ((dp)):
             j  ((dp[])):
                dp[i][j] =   i ==   j ==   (dp[i - ][j]  i  ) + (dp[i][j - ]  j  )
         dp[-][-]
     (, , a, b) * (a, b, R, C) % MOD

((, , , ))   
((, , , ))   
(((R, C, a, b) == (R, C, a, b)  R  (, )  C  (, )  a  (, R + )  b  (, C + )))

Worked problem

Numbers coprime to m, by inclusion–exclusion

The problem

Count the integers k with 1 ≤ k ≤ n that are coprime to m, where n and m can be as large as 10¹².

Looping over all k is impossible for n = 10¹². Instead count the *bad* numbers, those sharing a prime factor with m.
Factorise m into its distinct primes (at most 11 of them for m ≤ 10¹²).
A number is bad if it is divisible by at least one of those primes. Inclusion–exclusion: add the multiples of each prime, subtract the multiples of each pair (the product), add triples, and so on. Each term is n // product.
The answer is n minus the bad count; the signed sum over all subsets of primes (including the empty subset, which contributes +n) gives it directly.

python

from math import gcd

def distinct_primes(m):
    ps = []
    d = 2
    while d * d <= m:
        if m % d == 0:
            ps.(d)
             m % d == :
                m //= d
        d += 
     m > :
        ps.(m)
     ps

 (n, m):                       
    primes = (m)
    total = 
     mask  ( << (primes)):
        prod, bits = , 
         i, p  (primes):
             mask >> i & :
                prod *= p
                bits += 
        total += (-) ** bits * (n // prod)    
     total

((, ))   
((, ))   
((**, ))   
(((n, m) == ((k, m) ==   k  (, n + ))  n  (, )  m  (, )))

Wrap-up

Remember this

Remainders respect +, − and ×, but not division: divide by multiplying with a modular inverse.
Fermat gives a⁻¹ = a^(p−2) for a prime modulus; extended Euclid or pow(a, -1, m) works whenever gcd(a, m) = 1.
Precompute factorials and inverse factorials with a single modular exponentiation to answer C(n, r) mod p in O(1).
Pascal's rule needs no division and works for any modulus; totients and inclusion–exclusion count by divisibility.
Map counting problems onto formulas: paths and subsets (binomials), distributions (stars and bars), brackets (Catalan).

Cheat sheet

Inverse: a⁻¹ mod p = pow(a, p − 2, p) · any m: pow(a, −1, m) if gcd = 1
C(n, r) mod p: fact[n] · inv[r] · inv[n − r] with one pow() for the tables
Pascal's rule: C(n, r) = C(n−1, r−1) + C(n−1, r): any modulus, O(n²)
Stars and bars: n items, k boxes: C(n + k − 1, k − 1)
Catalan: C(2n, n) / (n + 1)
Inclusion–exclusion: sum over subsets of primes of (−1)^|S| · n // product(S)

Top traps to avoid

Reducing factorials mod p and then using / or //.
Computing one pow(..., p - 2, p) per query for every factorial.
Forgetting % MOD after each multiplication.

Mastery check

0 of 15 answered

Answer from understanding, not memory of the page. You need 12 of 15. Each answer is explained, and you can retry as often as you like.

1.Why can you not compute C(n, r) mod p by reducing each factorial mod p and then dividing normally?
2.What is the modular inverse of 3 modulo 11?
3.How many ways are there to split 5 identical sweets among 3 children (a child may get none)?
4.After precomputing factorial tables, what does one C(n, r) mod p query cost?
5.“Count the monotone grid paths from the top-left to the bottom-right modulo 10⁹ + 7” suggests:
6.When does the factorial-based method for C(n, r) mod p fail?
7.What does print(pow(3, 200, 13)) print?
8.What is the modular inverse of 3 modulo 7, as printed by print(pow(3, -1, 7))?
9.How many positive divisors does 36 have?
10.What does math.comb(10, 3) return? (Import math first.)
11.What does Fermat's little theorem say for a prime p and a not divisible by p?
12.What is the running time of fast modular exponentiation with exponent e?
13.What is the sum of the divisors of 12 (sigma of 12)?
14.Why are binomial coefficients modulo a prime computed with inverse factorials?
15.What is gcd(1071, 462)?

Number Theory & Combinatorics

After this lesson you will be able to

Builds on

Unlocks

How many ways? (and the answer will not fit in any variable)

Take remainders, then divide

Dividing by b means multiplying by b's inverse

The toolkit

Pascal's triangle

C(5, 2) three ways

Modular inverses

C(n, r) mod p in O(1) per query

Pascal's triangle for any modulus

Euler's totient with a sieve

Counting rules, checked by brute force

Costs

Edge cases

Mistakes to avoid

How to spot a counting problem

Paths through a required cell

Numbers coprime to m, by inclusion–exclusion

Wrap-up

Remember this

Cheat sheet

Top traps to avoid

Mastery check

Practice this topic

Core59 problems

Application1 problem

Reading the code

Reading the code

Reading the code

Reading the code

Reading the code