Lucas Theorem

Mathematical Statement

Given prime $p$, the following congruence holds:

$$\binom{n}{m} \equiv \binom{\lfloor n/p \rfloor}{\lfloor m/p \rfloor} \cdot \binom{n \bmod p}{m \bmod p} \pmod{p}$$

Proof Foundation

Lemma 1

For prime $p$:

$$\binom{p}{k} \equiv 0 \pmod{p} \text{ when } 0 < k < p$$

The numerator contains factor $p$, making the expression divisible by $p$ unless $k = 0$ or $k = p$.

Lemma 2

Using binomial theorem:

$$(x + y)^p \equiv x^p + y^p \pmod{p}$$

This follows from Fermat's Little Theorem where $a^p \equiv a \pmod{p}$ for prime $p$.

Complete Proof

Consider $(1 + x)^n$. The coefficient of $x^m$ is $\binom{n}{m}$.

By Lemma 2:

$$(1 + x)^n \equiv (1 + x^p)^{\lfloor n/p \rfloor} \cdot (1 + x)^{n \bmod p} \pmod{p}$$

To obtain term $x^m$ where $m = qp + r$:

• Select $q$ terms of $x^p$ from $(1 + x^p)^{\lfloor n/p \rfloor}$
• Select $r$ terms of $x$ from $(1 + x)^{n \bmod p}$

Thus coefficient becomes:

$$\binom{\lfloor n/p \rfloor}{q} \cdot \binom{n \bmod p}{r} \equiv \binom{\lfloor n/p \rfloor}{\lfloor m/p \rfloor} \cdot \binom{n \bmod p}{m \bmod p} \pmod{p}$$

Implementation Example

#include <iostream>
using namespace std;

const int MAX_PRIME = 100005;
typedef long long int64;

int64 factorial[MAX_PRIME], inverse[MAX_PRIME];

int64 combination(int64 n, int64 k, int64 mod) {
    if (k > n) return 0;
    return factorial[n] * inverse[k] % mod * inverse[n - k] % mod;
}

int64 lucas(int64 n, int64 k, int64 mod) {
    if (k == 0) return 1;
    return combination(n % mod, k % mod, mod) * lucas(n / mod, k / mod, mod) % mod;
}

void preprocess(int64 mod) {
    factorial[0] = factorial[1] = 1;
    inverse[0] = inverse[1] = 1;
    
    for (int i = 2; i <= mod; i++)
        factorial[i] = factorial[i-1] * i % mod;
        
    for (int i = 2; i <= mod; i++)
        inverse[i] = (mod - mod/i) * inverse[mod % i] % mod;
        
    for (int i = 2; i <= mod; i++)
        inverse[i] = inverse[i-1] * inverse[i] % mod;
}

int main() {
    int test_cases, n, m;
    int64 prime_mod;
    
    cin >> test_cases;
    while (test_cases--) {
        cin >> n >> m >> prime_mod;
        preprocess(prime_mod);
        cout << lucas(n + m, n, prime_mod) << endl;
    }
    return 0;
}

Extended Lucas Theorem

Problem Definition

Compute $\binom{n}{m} \bmod p$ where $p$ is not necessarily prime.

Solution Approach

Decompose $p$ into prime power factors: $p = \prod p_i^{a_i}$.

Since these factors are pairwise coprime, compute $\binom{n}{m} \bmod p_i^{a_i}$ separately and combine using Chinese Remainder Theorem.

Express binomial coefficient as:

$$\binom{n}{m} = \frac{n!}{m!(n-m)!}$$

Factor out powers of prime $q$ from factorial:

$$n! = q^{\lfloor n/q \rfloor} \cdot \lfloor n/q \rfloor! \cdot \prod_{q \nmid i, i \leq n} i$$

Components:

• $q^{\lfloor n/q \rfloor}$ computed via fast exponentiation
• $\lfloor n/q \rfloor!$ computed recursively
• Remaining product has periodicity modulo $q^{a_i}$

Periodic Product Computation

For modulus $q^{a}$, the product $\prod_{q \nmid i}^n i$ exhibits cycles of length less then $q^{a}$.

Example with $n = 19, q = 3, a = 2$:

$$19! \bmod 9 = (1 \times 2 \times 4 \times 5 \times 7 \times 8)^2 \times 19 \times 3^6 \times 6! \pmmod{9}$$

First cycle computed explicitly, remaining cycles handled via exponentiation.

Implementation

#include <bits/stdc++.h>
using namespace std;

typedef long long int64;

int64 extended_gcd(int64 a, int64 b, int64 &x, int64 &y) {
    if (b == 0) {
        x = 1; y = 0;
        return a;
    }
    int64 result = extended_gcd(b, a % b, x, y);
    tie(x, y) = make_pair(y, x - (a/b) * y);
    return result;
}

int64 modular_inverse(int64 value, int64 mod) {
    int64 x, y;
    extended_gcd(value, mod, x, y);
    return (x % mod + mod) % mod;
}

int64 fast_power(int64 base, int64 exp, int64 mod) {
    int64 result = 1;
    while (exp > 0) {
        if (exp & 1) result = result * base % mod;
        base = base * base % mod;
        exp >>= 1;
    }
    return result;
}

int64 factorial_mod(int64 n, int64 prime, int64 power) {
    if (n == 0) return 1;
    
    int64 cycle_result = 1;
    for (int i = 2; i <= power; i++)
        if (i % prime != 0)
            cycle_result = cycle_result * i % power;
            
    cycle_result = fast_power(cycle_result, n / power, power);
    
    for (int i = 2; i <= n % power; i++)
        if (i % prime != 0)
            cycle_result = cycle_result * i % power;
            
    return cycle_result * factorial_mod(n / prime, prime, power) % power;
}

int64 combination_prime_power(int64 n, int64 m, int64 prime, int64 power) {
    if (n < m) return 0;
    
    int64 numerator = factorial_mod(n, prime, power);
    int64 denominator_m = factorial_mod(m, prime, power);
    int64 denominator_nm = factorial_mod(n - m, prime, power);
    
    int64 prime_count = 0;
    for (int64 i = prime; i <= n; i *= prime) prime_count += n / i;
    for (int64 i = prime; i <= m; i *= prime) prime_count -= m / i;
    for (int64 i = prime; i <= n - m; i *= prime) prime_count -= (n - m) / i;
    
    return numerator * modular_inverse(denominator_m, power) % power * 
           modular_inverse(denominator_nm, power) % power * 
           fast_power(prime, prime_count, power) % power;
}

int64 chinese_remainder_combine(int64 residue, int64 mod_component) {
    static int64 total_mod;
    return residue * (total_mod / mod_component) % total_mod * 
           modular_inverse(total_mod / mod_component, mod_component) % total_mod;
}

int64 extended_lucas(int64 n, int64 m, int64 composite_mod) {
    int64 temp = composite_mod, result = 0;
    chinese_remainder_combine(0, 0); // Initialize static variable
    
    for (int64 factor = 2; factor * factor <= temp; factor++) {
        if (temp % factor == 0) {
            int64 prime_power = 1;
            while (temp % factor == 0) {
                prime_power *= factor;
                temp /= factor;
            }
            result = (result + chinese_remainder_combine(
                combination_prime_power(n, m, factor, prime_power), 
                prime_power)) % composite_mod;
        }
    }
    
    if (temp > 1)
        result = (result + chinese_remainder_combine(
            combination_prime_power(n, m, temp, temp), 
            temp)) % composite_mod;
            
    return result;
}

int main() {
    int64 n, m, mod;
    cin >> n >> m >> mod;
    cout << extended_lucas(n, m, mod) << endl;
    return 0;
}