Binomial Coefficients Modulo Composite Numbers

Intuition, Vocabulary and Foundations

For a prime modulus $p$, factorial-based binomial formulas often use $\binom nk\equiv n!\,(k!)^{-1}\,((n-k)!)^{-1}\pmod p$. This fails for a composite modulus when a denominator shares a factor with the modulus. Modulo $8$, for example, $2!$ has no inverse even though $\binom42=6$ is well defined.

The remedy separates every factorial into powers of a prime $p$ and a remaining unit invertible modulo $p^q$. Legendre’s formula gives $v_p(n!)=\sum_{j\ge1}\lfloor n/p^j\rfloor$, so the exponent of $p$ in the binomial coefficient is $e=v_p(n!)-v_p(k!)-v_p((n-k)!)$.

After removing those prime factors, the factorial units may be divided modulo $p^q$. If $e\ge q$, the binomial coefficient is divisible by $p^q$ and its local residue is zero.

A composite modulus has a unique prime-power decomposition $m=\prod_i p_i^{q_i}$. Compute the binomial residue modulo every prime power and reconstruct the unique result modulo $m$ with CRT. Each prime power supplies one independent local constraint.

Never cancel factorials modulo a composite number until all noninvertible prime factors have been separated. For each $p^q$, track the valuation, invert only unit parts, return zero when the valuation reaches $q$, and combine local answers with CRT.

Concept	Definition or invariant	Concrete example	Non-example or boundary	Question answered
Prime-power factorization	$m=\prod_i p_i^{q_i}$	$72=2^3\cdot3^2$	Treating $72$ as prime	Which local moduli must be solved?
$p$-adic valuation	$v_p(n!)=\sum_{j\ge1}\lfloor n/p^j\rfloor$	$v_2(10!)=5+2+1=8$	Counting only multiples of $p$ misses higher powers	How many noninvertible factors occur?
Unit factorial	Factorial product with factors of $p$ removed	Remove factors of $2$ modulo $8$	Omitting multiples of $p$ only once is insufficient	Which portion may be inverted?
Binomial valuation	$v_p(n!)-v_p(k!)-v_p((n-k)!)$	$v_2(\binom42)=1$	Valuations are not modular inverses	What power of $p$ remains?
Zero residue	Valuation at least $q$	A multiple of $8$ is zero modulo $2^3$	Unit factorials alone cannot detect this	When may computation stop early?
Local residue	$\binom nk\bmod p^q$	Solve modulo $8$ and $9$ independently	Local answers are not yet modulo $72$	What does one prime power require?
CRT reconstruction	Unique class modulo the product	Combine residues modulo $8$ and $9$	Components must be coprime	How are local answers recombined?
Query preprocessing	Prefix products excluding multiples of $p$	Reuse one table for many queries	A table of size $p^q$ may be too large	When is preprocessing worthwhile?

Step-by-Step Derivations and Algorithms

Factor the modulus into prime powers. For one component $P=p^q$, compute $e=v_p(n!)-v_p(k!)-v_p((n-k)!)$. If $e\ge q$, return zero modulo $P$.

Otherwise compute factorials with all factors of $p$ removed. Let $U_P(t)=\prod_{1\le i\le t,\ p\nmid i}i\pmod P$ for $0\le t\le P$. Products over complete blocks of length $P$ repeat modulo $P$. Factors divisible by $p$ are handled recursively after division by $p$, giving a logarithmic-depth unit-factorial routine.

If $F_p(n)$ denotes the factorial unit, then

$$ \binom nk \equiv F_p(n)F_p(k)^{-1}F_p(n-k)^{-1}p^e \pmod P. $$

Both denominator units are coprime to $P$, so their inverses exist. This is the exact point where modular division becomes valid.

After computing residues $a_i$ modulo pairwise coprime $P_i=p_i^{q_i}$, merge them incrementally. If $x\equiv a\pmod M$ and $x\equiv b\pmod P$, solve $Mt\equiv b-a\pmod P$ and replace $a$ by $a+Mt$ modulo $MP$. The merged class satisfies every processed component.

For one modulus and many queries, factor it and build each unit-prefix table once. Query cost becomes logarithmic in $n$ for each component plus a small CRT reconstruction. Table size depends on the modulus, not on $n$.

Binomial Coefficients Modulo Composite Numbers algorithm map

flowchart TD
    A["Input n, k, and composite m"] --> B{"k outside valid range?"}
    B -->|Yes| C["Return zero"]
    B -->|No| D["Factor m into prime powers"]
    D --> E["Select p^q"]
    E --> F["Compute binomial valuation"]
    F --> G{"Valuation at least q?"}
    G -->|Yes| H["Local residue is zero"]
    G -->|No| I["Compute unit factorials"]
    I --> J["Invert unit denominators"]
    J --> K["Restore p to valuation power"]
    H --> L{"More prime powers?"}
    K --> L
    L -->|Yes| E
    L -->|No| M["Reconstruct with CRT"]

For each prime power, the algorithm preserves the exact exponent of $p$ separately from a unit invertible modulo $p^q$; multiplying them reconstructs the local binomial residue. With unit prefixes totaling $\sum p_i^{q_i}$ entries, each query uses $O(\sum_i\log_{p_i} n)$ recursive steps plus $O(r)$ CRT merges for $r$ components.

Implementation, Complexity and Engineering

The reference implementation supports large $n$ when the prime-power components are small enough to precompute. It caches one prefix table per prime power, evaluates unit factorials recursively, and reconstructs with pairwise-coprime CRT.

The main engineering boundary is explicit: $O(p^q)$ preprocessing is unsuitable when a prime-power component is enormous. In that regime, use a segmented or product-tree method, exploit a small query count, or consider direct exact computation.

from math import comb, gcd


def factor_prime_powers(modulus):
    if modulus <= 0:
        raise ValueError("modulus must be positive")
    factors = []
    remaining = modulus
    prime = 2
    while prime * prime <= remaining:
        if remaining % prime == 0:
            power = 1
            exponent = 0
            while remaining % prime == 0:
                remaining //= prime
                power *= prime
                exponent += 1
            factors.append((prime, exponent, power))
        prime = 3 if prime == 2 else prime + 2
    if remaining > 1:
        factors.append((remaining, 1, remaining))
    return factors


def valuation_factorial(n, prime):
    total = 0
    while n:
        n //= prime
        total += n
    return total


class PrimePowerBinomial:
    def __init__(self, prime, exponent):
        if prime < 2 or exponent < 1:
            raise ValueError("invalid prime power")
        self.prime = prime
        self.exponent = exponent
        self.modulus = prime**exponent
        self.unit_prefix = [1] * (self.modulus + 1)
        for value in range(1, self.modulus + 1):
            self.unit_prefix[value] = self.unit_prefix[value - 1]
            if value % prime != 0:
                self.unit_prefix[value] = (
                    self.unit_prefix[value] * value
                ) % self.modulus
        self.block_product = self.unit_prefix[self.modulus]

    def unit_factorial(self, n):
        if n == 0:
            return 1
        complete_blocks, remainder = divmod(n, self.modulus)
        result = pow(self.block_product, complete_blocks, self.modulus)
        result = result * self.unit_prefix[remainder] % self.modulus
        return result * self.unit_factorial(n // self.prime) % self.modulus

    def choose(self, n, k):
        if n < 0:
            raise ValueError("n must be nonnegative")
        if k < 0 or k > n:
            return 0
        valuation = (
            valuation_factorial(n, self.prime)
            - valuation_factorial(k, self.prime)
            - valuation_factorial(n - k, self.prime)
        )
        if valuation >= self.exponent:
            return 0
        numerator = self.unit_factorial(n)
        denominator = (
            self.unit_factorial(k) * self.unit_factorial(n - k)
        ) % self.modulus
        unit_part = numerator * pow(denominator, -1, self.modulus) % self.modulus
        return unit_part * pow(self.prime, valuation, self.modulus) % self.modulus


def crt_pairwise(residues, moduli):
    value = 0
    combined_modulus = 1
    for residue, modulus in zip(residues, moduli):
        if gcd(combined_modulus, modulus) != 1:
            raise ValueError("CRT moduli must be pairwise coprime")
        step = (
            (residue - value) * pow(combined_modulus, -1, modulus)
        ) % modulus
        value += combined_modulus * step
        combined_modulus *= modulus
        value %= combined_modulus
    return value, combined_modulus


class CompositeBinomial:
    def __init__(self, modulus, table_limit=2_000_000):
        if modulus <= 0:
            raise ValueError("modulus must be positive")
        self.modulus = modulus
        if modulus == 1:
            self.components = []
            return
        factors = factor_prime_powers(modulus)
        if any(power > table_limit for _, _, power in factors):
            raise ValueError("prime-power table would exceed limit")
        self.components = [
            PrimePowerBinomial(prime, exponent)
            for prime, exponent, _ in factors
        ]

    def choose(self, n, k):
        if n < 0:
            raise ValueError("n must be nonnegative")
        if k < 0 or k > n or self.modulus == 1:
            return 0
        residues = [component.choose(n, k) for component in self.components]
        moduli = [component.modulus for component in self.components]
        value, rebuilt_modulus = crt_pairwise(residues, moduli)
        assert rebuilt_modulus == self.modulus
        return value


for modulus in range(1, 101):
    calculator = CompositeBinomial(modulus)
    for n in range(31):
        for k in range(-1, n + 2):
            expected = 0 if k < 0 or k > n else comb(n, k) % modulus
            assert calculator.choose(n, k) == expected

calculator_72 = CompositeBinomial(72)
assert calculator_72.choose(4, 2) == 6
assert calculator_72.choose(100, 50) == comb(100, 50) % 72

large_n = 10**18 + 123
calculator_1000 = CompositeBinomial(1000)
assert calculator_1000.choose(large_n, 0) == 1
assert calculator_1000.choose(large_n, 1) == large_n % 1000
assert calculator_1000.choose(large_n, large_n) == 1

mod_8 = PrimePowerBinomial(2, 3)
assert mod_8.choose(8, 1) == 0
assert mod_8.choose(4, 2) == 6

try:
    CompositeBinomial(0)
    raise AssertionError("zero modulus was accepted")
except ValueError:
    pass

print("all tests passed")

Implementation choice	Best workload	Cost	Failure mode	Mitigation
Exact `comb` then reduce	Small $n$, few queries	Depends on exact integer size	Huge intermediate integer	Keep as an oracle or baseline
Prime-power unit tables	Large $n$, repeated queries, moderate $p^q$	$O(p^q)$ preprocessing	Table exceeds memory	Enforce a table-size budget
Recursive unit factorial	Large $n$ after preprocessing	$O(\log_p n)$ per unit	Omitting recursive multiples of $p$	Test dense small prefixes
Valuation tracking	Every local query	$O(\log_p n)$	Inverting a denominator divisible by $p$	Remove powers first
Zero-residue shortcut	High valuations	Constant after valuation	Returning zero when $e\lt q$	Compare with the exponent $q$
CRT reconstruction	Composite modulus factors	Linear in component count	Combining non-coprime components	Use distinct prime powers
Fixed-width implementation	Bounded modulus	Fast arithmetic	Products overflow before reduction	Use wide or safe multiplication
Cached calculator	Many queries for one modulus	Amortized setup	Rebuilding tables per query	Cache factorization and tables

Applications, Case Studies and Exercises

Composite-modulus binomial coefficients arise in combinatorial dynamic programming, periodic counting, polynomial algorithms, and systems using nonprime machine-friendly moduli. A huge $n$ can remain manageable when every prime-power component is modest.

Practice

Apply Binomial Coefficients Modulo Composite Numbers

Case study 1 — compute $\binom{10}{4}\bmod72$. Factor $72=2^3\cdot3^2$. Compute valuation and unit parts modulo $8$ and $9$, then combine them with CRT. Checkpoint: compare the result with $210\bmod72=66$. Success means no factor divisible by $2$ or $3$ is inverted.

Case study 2 — repeated large-$n$ queries. Build one calculator for $1000=2^3\cdot5^3$. Evaluate thousands of queries with $n$ up to $10^{18}$. Record preprocessing memory, latency, zero local residues, and CRT time. Verify $\binom n0=\binom nn=1$ and $\binom n1\equiv n\pmod{1000}$, then compare smaller queries with exact arithmetic.

Hands-on differential test. For every modulus from $1$ through $200$, every $0\le n\le60$, and every relevant $k$, compare with math.comb(n, k) % m. Include modulus one, prime powers, squarefree values, highly composite values, and invalid $k$. Success means zero mismatches and no accidental inversion of nonunits.

Tiered exercises. Recall: explain why factorial division fails modulo $8$. Derive: prove the valuation formula. Implement: cache calculators by modulus. Debug: remove the valuation-zero shortcut and locate a false nonzero result. Extend: replace the full $p^q$ table with segmented products for a large component.

Factor the modulus first, solve each prime power by separating valuations from invertible unit factorials, and reconstruct with CRT. Use exact arithmetic for small one-off inputs, cached prime-power tables for repeated large-$n$ queries with moderate moduli, and another product strategy when a prime-power table is too large.

Cookie Consent