> ## Documentation Index
> Fetch the complete documentation index at: https://resources.devweekends.com/llms.txt
> Use this file to discover all available pages before exploring further.

# Number Theory for CP

> Master modular arithmetic, primes, combinatorics, and mathematical techniques for CP

# 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+.

<Note>
  **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
</Note>

***

## 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 $2^{1000000}$, we keep everything "wrapped around" within $[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) \mod m = ((a \mod m) + (b \mod m)) \mod m$
* $(a \times b) \mod m = ((a \mod m) \times (b \mod m)) \mod m$
* $(a - b) \mod m = ((a \mod m) - (b \mod m) + m) \mod m$ (add $m$ to handle negatives)

```cpp theme={null}
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; }
```

<Warning>
  **Division is Different!**
  `(a / b) % m ≠ ((a % m) / (b % m)) % m`

  For division, multiply by **modular inverse**: `a / b ≡ a × b^(-1) (mod m)`
</Warning>

### Binary Exponentiation

**The Problem**: Computing $a^n$ directly takes O(n) multiplications. For $n = 10^{18}$, that's impossible.

**The Insight**: Use the property that:

* If n is even: $a^n = (a^{n/2})^2$
* If n is odd: $a^n = a \times a^{n-1}$

**Example**: Computing $3^{13}$:

* $3^{13} = 3 \times 3^{12}$
* $3^{12} = (3^6)^2$
* $3^6 = (3^3)^2$
* $3^3 = 3 \times 3^2$
* $3^2 = (3^1)^2$
* $3^1 = 3 \times 3^0$

This reduces O(n) to **O(log n)** multiplications!

**Bit Interpretation**: 13 in binary is 1101. We compute $3^1, 3^2, 3^4, 3^8$ and multiply those where the bit is 1: $3^{13} = 3^8 \times 3^4 \times 3^1$.

```cpp theme={null}
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 $a$ modulo $m$ is a number $a^{-1}$ such that:
$a \times a^{-1} \equiv 1 \pmod{m}$

Think of it as "division in modular world." To compute $\frac{x}{a} \mod m$, we compute $x \times a^{-1} \mod m$.

**Fermat's Little Theorem**: For prime $p$ and $a$ not divisible by $p$:
$a^{p-1} \equiv 1 \pmod{p}$

Rearranging: $a \times a^{p-2} \equiv 1 \pmod{p}$, so $a^{-1} \equiv a^{p-2} \pmod{p}$

**When does inverse exist?**

* For prime modulus: Always exists if $\gcd(a, p) = 1$ (i.e., $a \not\equiv 0$)
* For general modulus: Only if $\gcd(a, m) = 1$

```cpp theme={null}
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, y$ such that $ax + by = \gcd(a, b)$

**Why it matters**: This is **Bézout's Identity**—such $x, y$ always exist!

**How it works**:

1. Base case: When $b = 0$, we have $\gcd(a, 0) = a$, so $x = 1, y = 0$ works: $a(1) + 0(0) = a$
2. Recursive case: If we know $x', y'$ for $(b, a \mod b)$ such that $bx' + (a \mod b)y' = g$, then:
   * $a \mod b = a - \lfloor a/b \rfloor \cdot b$
   * Substituting: $bx' + (a - \lfloor a/b \rfloor \cdot b)y' = g$
   * Rearranging: $ay' + b(x' - \lfloor a/b \rfloor \cdot y') = g$
   * So $x = y'$ and $y = x' - \lfloor a/b \rfloor \cdot y'$

**Application**: Finding modular inverse when modulus is not prime. If $\gcd(a, m) = 1$, then from $ax + my = 1$, we get $ax \equiv 1 \pmod{m}$, so $x$ is the inverse.

```cpp theme={null}
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 $i$, start crossing from $i^2$ (smaller multiples already crossed by smaller primes)
5. Repeat until $i^2 > N$

**Why O(N log log N)?** Each prime $p$ crosses out $N/p$ numbers. The sum $\sum_{p \leq N} N/p = N \sum_{p \leq N} 1/p$ is $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
```

```cpp theme={null}
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 $i$:

* We multiply $i$ by primes $p \leq \text{spf}[i]$
* This ensures $i \cdot p$ is marked by its SPF, which is $p$
* We stop when $p > \text{spf}[i]$ because then $i \cdot p$ would have SPF = $\text{spf}[i]$, not $p$

**Why SPF is useful**:

* O(log n) factorization: repeatedly divide by SPF
* Euler's totient: compute multiplicatively using prime factorization

```cpp theme={null}
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))

```cpp theme={null}
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 x p2^a2 x ... x pk^ak, then:

* Number of divisors = (a1+1) x (a2+1) x ... x (ak+1)
* Sum of divisors = product of ((pi^(ai+1) - 1) / (pi - 1)) for each prime

**Why?** Each divisor is formed by choosing an exponent for each prime factor independently. For prime p with exponent a, we can choose 0, 1, 2, ..., or a copies of p. That is (a+1) choices. Since the choices for different primes are independent, we multiply them.

**Example**: 12 = 2^2 x 3^1. Divisor count = (2+1) x (1+1) = 6. Divisors: 1, 2, 3, 4, 6, 12.

<Tip>
  **Contest tip**: The maximum number of divisors of any n up to 10^9 is 1344 (the number 735134400). Up to 10^18, it is 103680. This means iterating over all divisors is always feasible if you can find them in O(sqrt(n)).
</Tip>

```cpp theme={null}
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 $\binom{n}{r} = \frac{n!}{r!(n-r)!}$ for many queries.

**Naive approach**: Compute each query in O(n)—too slow for $10^5$ queries.

**Smart approach**: Precompute $n!$ and $(n!)^{-1}$ for all $n$.

**The Trick for Inverse Factorials**: Instead of computing $n$ modular inverses (each costs O(log MOD)):

1. Compute $(MAXN-1)!^{-1}$ using Fermat's theorem
2. Use the recurrence: $(n!)^{-1} = (n+1)! \cdot (n+1)^{-1}$... but even simpler:
   * $(n!)^{-1} = ((n+1)!)^{-1} \cdot (n+1)$
   * This is because $(n+1)! = n! \cdot (n+1)$, so $(n!)^{-1} = ((n+1)!)^{-1} \cdot (n+1)$

This gives **O(N) precomputation** and **O(1) per query**.

```cpp theme={null}
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

| Formula        | Expression          | Use Case                                   |
| -------------- | ------------------- | ------------------------------------------ |
| nCr            | n! / (r! × (n-r)!)  | Choose r from n                            |
| nPr            | n! / (n-r)!         | Arrange r from n                           |
| Stars and Bars | C(n+r-1, r-1)       | Distribute n identical items into r groups |
| Catalan(n)     | C(2n, n) / (n+1)    | Valid parentheses, BST count               |
| Derangements   | n! × Σ((-1)^i / i!) | Permutations with no fixed points          |

### Pascal's Triangle (Small n, r)

```cpp theme={null}
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 $\binom{n}{r} \mod p$ when $n$ is huge (like $10^{18}$) but $p$ is a small prime.

**Why standard approach fails**: We can't precompute factorials up to $10^{18}$.

**Lucas Theorem**:
$\binom{n}{r} \equiv \prod_{i=0}^{k} \binom{n_i}{r_i} \pmod{p}$

where $n = n_k p^k + ... + n_1 p + n_0$ and $r = r_k p^k + ... + r_1 p + r_0$ are base-$p$ 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**: $\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)
* $\binom{12}{5} \equiv \binom{1}{0} \cdot \binom{1}{1} \cdot \binom{0}{2} = 1 \cdot 1 \cdot 0 = 0 \pmod{3}$

```cpp theme={null}
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

```cpp theme={null}
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

```cpp theme={null}
ll gcdArray(vector<ll>& arr) {
    ll g = 0;  // gcd(0, x) = x, so starting with 0 is correct
    for (ll x : arr) g = gcd(g, x);
    return g;
}
```

### Properties

* gcd(a, b) x lcm(a, b) = a x b (useful for computing LCM without overflow: `a / gcd(a,b) * b`)
* gcd(a, b, c) = gcd(gcd(a, b), c) (associative -- works for arrays)
* gcd(a, 0) = a (identity element -- this is why the array GCD starts at 0)
* If gcd(a, b) = g, then gcd(a/g, b/g) = 1 (coprime after dividing by GCD)
* gcd(a, b) = gcd(a - b, b) for a > b (subtractive form, basis of Euclidean algorithm)
* gcd does not change under modular reduction: gcd(a, b) = gcd(a mod b, b)

<Tip>
  **Contest tip**: When computing LCM of an array, always divide before multiplying to avoid overflow: `lcm = lcm / gcd(lcm, arr[i]) * arr[i]`. Also, if any element is 0, the LCM is 0 (or undefined, depending on convention) -- check the problem statement.
</Tip>

***

## Euler's Totient Function

**Definition**: $\phi(n)$ = count of integers in $[1, n]$ that are coprime to $n$ (i.e., $\gcd(k, n) = 1$).

**Examples**:

* $\phi(1) = 1$ (just 1 itself)
* $\phi(6) = 2$ (1 and 5 are coprime to 6)
* $\phi(7) = 6$ (for prime p: all of 1, 2, ..., p-1 are coprime)
* $\phi(p) = p - 1$ for any prime $p$

**Formula**: If $n = p_1^{a_1} \times p_2^{a_2} \times ... \times p_k^{a_k}$, then:
$\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 $n$ numbers. For each prime factor $p$, exactly $1/p$ of numbers are divisible by $p$, so we keep $(1 - 1/p)$ of them.

**Example**: $\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 ✓

```cpp theme={null}
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^phi(n) = 1 (mod n)

**Application**: a^b mod m when b is huge

* If gcd(a, m) = 1: a^b = a^(b mod phi(m)) (mod m)

**Why this matters**: When b is astronomically large (like 10^(10^18)), you cannot even store b. But Euler's theorem lets you reduce the exponent: compute b mod phi(m) first, then use binary exponentiation on the reduced exponent.

**Extended version** (when gcd(a, m) != 1): If b >= log2(m), then a^b = a^(phi(m) + b mod phi(m)) (mod m). This handles cases where a and m share common factors.

<Warning>
  **Edge case**: When b = 0, a^b = 1 regardless of a (by convention, even 0^0 = 1 in most CP contexts). However, some problems define 0^0 differently--always check the problem statement.
</Warning>

***

## Pattern: Counting with Inclusion-Exclusion

**The Problem**: Count elements satisfying **at least one** of k conditions.

**The Principle**:
$|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):

* $|A \cup B|$ counts elements in either A or B
* $|A| + |B|$ overcounts elements in both by 1
* So $|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]$ divisible by at least one of the given primes.

* Divisible by $p$: there are $\lfloor n/p \rfloor$ such numbers
* Divisible by both $p$ and $q$: there are $\lfloor n/(p \cdot q) \rfloor$ such numbers

```cpp theme={null}
// 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)

| Problem           | Topic         | Link                                                         |
| ----------------- | ------------- | ------------------------------------------------------------ |
| Exponentiation    | Binary exp    | [CSES](https://cses.fi/problemset/task/1095)                 |
| Counting Divisors | Factorization | [CSES](https://cses.fi/problemset/task/1713)                 |
| GCD Query         | GCD           | [CF 1350C](https://codeforces.com/problemset/problem/1350/C) |

### Intermediate (1300-1600)

| Problem               | Topic         | Link                                         |
| --------------------- | ------------- | -------------------------------------------- |
| Binomial Coefficients | nCr           | [CSES](https://cses.fi/problemset/task/1079) |
| Creating Strings II   | Combinatorics | [CSES](https://cses.fi/problemset/task/1715) |
| Sum of Divisors       | Number theory | [CSES](https://cses.fi/problemset/task/1082) |

### Advanced (1600-1900)

| Problem                | Topic        | Link                                         |
| ---------------------- | ------------ | -------------------------------------------- |
| Counting Coprime Pairs | Euler phi    | [CSES](https://cses.fi/problemset/task/2417) |
| Christmas Party        | Derangements | [CSES](https://cses.fi/problemset/task/1717) |

***

## Key Takeaways

<CardGroup cols={2}>
  <Card title="Modular Inverse" icon="divide">
    Use a^(p-2) for prime p, or extended GCD for general.
  </Card>

  <Card title="Precompute Factorials" icon="table">
    O(n) precomputation gives O(1) nCr queries.
  </Card>

  <Card title="SPF Sieve" icon="filter">
    Linear sieve enables O(log n) factorization.
  </Card>

  <Card title="Inclusion-Exclusion" icon="circle-nodes">
    Count by adding/subtracting overlaps.
  </Card>
</CardGroup>

***

## Next Up

<Card title="Chapter 24: Advanced Graph Algorithms" icon="arrow-right" href="/courses/competitive-programming/24-advanced-graphs">
  Master SCC, bridges, articulation points, and advanced graph techniques.
</Card>