viii. Factorials and Binomials mod p

How to get $n!$ and $\left(\genfrac{}{}{0ex}{}{n}{k}\right)$ modulo $m$ in every regime. With a prime $p>n$ use factorial tables (Modular Arithmetic, Comb); the rest is here.

Situation	Method	Cost
$m=p$ prime, $n<p$, many queries	factorial + inverse factorial tables	$O(n)$ build, $O(1)$ query
$m=p$ prime, $n<p$ huge ($p\sim 10^9$), many queries	fast factorial (shift of sampling points)	$\approx O(\sqrt{p}\log p)$ build
$m=p$ small prime, $n$ huge	Lucas	$O(p+{\log }_{p}n)$
$n!$ with the $p$-factors removed, $n\ge p$	Wilson recursion	$O(p+{\log }_{p}n)$
$m=p^e$	generalised Lucas (Granville): $n!=p^v\cdot F(n)$	$O(p^e+\log n)$
arbitrary $m$	factor $m$, solve each $p^e$, CRT	$O(\sum p^e)$ build
exponent of $p$ in $n!$ / $\left(\genfrac{}{}{0ex}{}{n}{k}\right)$	Legendre / Kummer	$O({\log }_{p}n)$

Legendre, Kummer, trailing zeros

$v_p(n!)={\sum }_{i\ge 1}\lfloor n/p^i\rfloor =\frac{n-s_p(n)}{p-1}$ ($s_p$ = digit sum in base $p$). Kummer: $v_p\left(\genfrac{}{}{0ex}{}{n}{k}\right)$ = number of carries when adding $k$ and $n-k$ in base $p$. Trailing zeros of $n!$ in base $b=\prod q_i^{a_i}$: ${\min }_i\lfloor v_{q_i}(n!)/a_i\rfloor$.

ll legendreExp(ll n, ll p) { ll e = 0; while (n) n /= p, e += n; return e; }                // v_p(n!)
ll kummerCarries(ll n, ll k, ll p) { ll c = 0, carry = 0, a = k, b = n - k; while (a || b || carry) { ll d = a % p + b % p + carry; carry = d >= p; c += carry; a /= p, b /= p; } return c; }
ll trailingZerosBase(ll n, ll b) {                          // zeros at the end of n! written in base b
    ll res = LLONG_MAX;
    for (ll q = 2; q * q <= b; q++) if (b % q == 0) { int a = 0; while (b % q == 0) b /= q, a++; res = min(res, legendreExp(n, q) / a); }
    if (b > 1) res = min(res, legendreExp(n, b));
    return res;
}
ll smallestWithZeros(ll z) {                               // smallest n with exactly z zeros at the end of n! (base 10), -1 if none
    ll lo = 0, hi = 5 * (z + 1);
    while (lo < hi) { ll mid = (lo + hi) / 2; if (legendreExp(mid, 5) >= z) hi = mid; else lo = mid + 1; }
    return legendreExp(lo, 5) == z ? lo : -1;
}

Proof: Legendre and Kummer

$v_p(n!)={\sum }_{j\le n}{v}_p(j)={\sum }_{j\le n}{\sum }_{i\ge 1}[p^i\mid j]={\sum }_i\lfloor n/p^i\rfloor$. With $n=\sum d_i{p}^i$, $\lfloor n/p^i\rfloor ={\sum }_{t\ge i}{d}_t{p}^{t-i}$, and summing over $i$ gives ${\sum }_t{d}_t\frac{p^t-1}{p-1}=\frac{n-s_p(n)}{p-1}$. Kummer: $v_p\left(\genfrac{}{}{0ex}{}{n}{k}\right)=\frac{s_p(k)+s_p(n-k)-s_p(n)}{p-1}$, and each carry in the addition $k+(n-k)$ lowers the digit sum by exactly $p-1$. Base $b$: a trailing zero in base $b$ is a factor $b$, limited by the scarcest prime power; for base 10, $v_2\ge v_5$ so only 5 matters, and the count of zeros jumps by $v_5(n)$ at multiples of 5, which is why some $z$ are skipped (-1).

Lucas and factorials with $p$ removed

$\left(\genfrac{}{}{0ex}{}{n}{k}\right)\equiv {\prod }_i\left(\genfrac{}{}{0ex}{}{n_i}{k_i}\right)(\mathrm{mod}p)$ over base-$p$ digits. For $n\ge p$: $n!=p^{v_p(n!)}\cdot F_p(n)$ with $F_p(n)\equiv (-1{)}^{\lfloor n/p\rfloor }(n\mathrm{mod}p)!F_p(\lfloor n/p\rfloor )(\mathrm{mod}p)$ (Wilson).

struct LucasP {                                            // small prime p, n and k up to 1e18
    ll p; vector<ll> f, fi;
    LucasP(ll p) : p(p), f(p), fi(p) {
        f[0] = 1; for (ll i = 1; i < p; i++) f[i] = f[i - 1] * i % p;
        fi[p - 1] = power(f[p - 1], p - 2, p); for (ll i = p - 1; i > 0; i--) fi[i - 1] = fi[i] * i % p;
    }
    ll smallC(ll n, ll k) const { return k < 0 || k > n ? 0 : f[n] * fi[k] % p * fi[n - k] % p; }
    ll C(ll n, ll k) const { ll r = 1; for (; (n || k) && r; n /= p, k /= p) r = r * smallC(n % p, k % p) % p; return k < 0 ? 0 : r; }
    ll factNoP(ll n) const { ll r = 1; for (; n > 0; n /= p) { r = r * f[n % p] % p; if ((n / p) & 1) r = (p - r) % p; } return r; }   // n! / p^v mod p
};

Proof: Lucas and the Wilson recursion

Lucas: in ${\mathbb{F}}_p[x]$, $(1+x{)}^p=1+x^p$, so $(1+x{)}^n={\prod }_i(1+x^{p^i}{)}^{n_i}$. The coefficient of $x^k=\prod x^{k_i{p}^i}$ on the right is ${\prod }_i\left(\genfrac{}{}{0ex}{}{n_i}{k_i}\right)$, because base-$p$ representations with digits $<p$ are unique. $F_p$: the numbers $\le n$ not divisible by $p$ form $\lfloor n/p\rfloor$ full blocks of residues $\mathrm{1..}p-1$ (each contributing $(p-1)!\equiv -1$ by Wilson) and a partial block $(n\mathrm{mod}p)!$; the multiples $p,2p,\dots$ contribute $p^{\lfloor n/p\rfloor }\lfloor n/p\rfloor !$, whose $p$-free part is $F_p(\lfloor n/p\rfloor )$.

Fast factorial for a large prime

$n!\mathrm{mod}p$ for $p$ up to $2^{31}$ and many queries: with block size $B$, $h(i)={\prod }_{j=1}^B(Bi+j)$ is a polynomial of degree $B$ in $i$; its values at $i=\mathrm{0..}B$ are computed directly and extended to all $i<p/B$ by shifting sampling points (Lagrange interpolation written as a convolution, done with a 3-prime NTT so any prime modulus works). Prefix products give $(kB)!$ and a query finishes the partial block in $O(B)$. Build $\approx O(B^2+\frac{p}{B}\log \frac{p}{B})$, about 0.15 s for $p\approx 10^9$.

// usage: FastFactorial<mint> F; F(n) = n! mod P (0 for n >= P); F.binom(n, k)
// ---------- modular integer (number-theory/modint.cpp) ----------
template <unsigned P>
struct MInt {
    unsigned x = 0;
    static inline unsigned dynamicMod = 998244353;
 
    static constexpr unsigned mod() {
        if constexpr (P > 0) return P;
        else return dynamicMod;
    }
    static void setMod(unsigned m) {
        static_assert(P == 0, "setMod is only for MInt<0>");
        dynamicMod = m;
    }
 
    constexpr MInt() {}
    constexpr MInt(i64 v) {
        v %= (i64)mod();
        if (v < 0) v += mod();
        x = (unsigned)v;
    }
 
    constexpr unsigned val() const { return x; }
    constexpr explicit operator unsigned() const { return x; }
 
    constexpr MInt operator-() const {
        MInt r;
        r.x = x ? mod() - x : 0;
        return r;
    }
    constexpr MInt pow(i64 e) const {
        MInt base = *this, res = 1;
        if (e < 0) base = base.inv(), e = -e;
        for (; e; e >>= 1, base *= base) {
            if (e & 1) res *= base;
        }
        return res;
    }
    constexpr MInt inv() const {
        i64 a = x, b = mod(), u = 1, v = 0;
        while (b) {
            i64 t = a / b;
            a -= t * b; swap(a, b);
            u -= t * v; swap(u, v);
        }
        assert(a == 1 && "value is not invertible modulo mod()");
        return MInt(u);
    }
 
    constexpr MInt& operator+=(const MInt& o) {
        x += o.x;
        if (x >= mod()) x -= mod();
        return *this;
    }
    constexpr MInt& operator-=(const MInt& o) {
        x += mod() - o.x;
        if (x >= mod()) x -= mod();
        return *this;
    }
    constexpr MInt& operator*=(const MInt& o) {
        x = (unsigned)((u64)x * o.x % mod());
        return *this;
    }
    constexpr MInt& operator/=(const MInt& o) { return *this *= o.inv(); }
 
    friend constexpr MInt operator+(MInt a, const MInt& b) { return a += b; }
    friend constexpr MInt operator-(MInt a, const MInt& b) { return a -= b; }
    friend constexpr MInt operator*(MInt a, const MInt& b) { return a *= b; }
    friend constexpr MInt operator/(MInt a, const MInt& b) { return a /= b; }
    friend constexpr bool operator==(const MInt& a, const MInt& b) { return a.x == b.x; }
    friend constexpr bool operator!=(const MInt& a, const MInt& b) { return a.x != b.x; }
 
    friend istream& operator>>(istream& in, MInt& a) {
        i64 v;
        in >> v;
        a = MInt(v);
        return in;
    }
    friend ostream& operator<<(ostream& out, const MInt& a) { return out << a.x; }
};
 
using mint = MInt<998244353>;   // change to MInt<1000000007> if the problem asks for 1e9+7
using dmint = MInt<0>;
 
// ---------- number theoretic transform (polynomial/ntt.cpp) ----------
constexpr unsigned nttPowMod(u64 b, u64 e, unsigned m) {
    u64 r = 1;
    for (b %= m; e; e >>= 1, b = b * b % m) {
        if (e & 1) r = r * b % m;
    }
    return (unsigned)r;
}
 
// smallest quadratic non-residue g mod P; g^((P - 1) / 2^k) is then a primitive 2^k-th root of unity.
template <unsigned P> constexpr unsigned nttNonResidue() {
    unsigned g = 2;
    while (nttPowMod(g, (P - 1) / 2, P) != P - 1) g++;
    return g;
}
 
template <unsigned P> void ntt(vector<MInt<P>>& a, bool invert = false) {
    using Z = MInt<P>;
    int n = a.size();
    if (n <= 1) return;
    assert((n & (n - 1)) == 0 && (P - 1) % n == 0);
    // roots[k + j] = w^j for a primitive 2k-th root of unity w (k a power of two, 0 <= j < k).
    static vector<Z> roots{0, 1};
    if ((int)roots.size() < n) {
        constexpr unsigned g = nttNonResidue<P>();
        int k = __builtin_ctz(roots.size());
        roots.resize(n);
        for (; (1 << k) < n; k++) {
            Z e = Z(g).pow((P - 1) >> (k + 1));
            for (int i = 1 << (k - 1); i < (1 << k); i++) {
                roots[2 * i] = roots[i];
                roots[2 * i + 1] = roots[i] * e;
            }
        }
    }
    for (int i = 1, j = 0; i < n; i++) {
        int bit = n >> 1;
        for (; j & bit; bit >>= 1) j ^= bit;
        j ^= bit;
        if (i < j) swap(a[i], a[j]);
    }
    for (int k = 1; k < n; k *= 2) {
        for (int i = 0; i < n; i += 2 * k) {
            for (int j = 0; j < k; j++) {
                Z u = a[i + j], v = a[i + j + k] * roots[k + j];
                a[i + j] = u + v;
                a[i + j + k] = u - v;
            }
        }
    }
    if (invert) {
        reverse(a.begin() + 1, a.end());
        Z inv = Z(n).inv();
        for (Z& x : a) x *= inv;
    }
}
 
template <unsigned P> vector<MInt<P>> convolution(vector<MInt<P>> a, vector<MInt<P>> b) {
    if (a.empty() || b.empty()) return {};
    int need = a.size() + b.size() - 1;
    if (min(a.size(), b.size()) <= 32) {
        vector<MInt<P>> c(need);
        for (int i = 0; i < (int)a.size(); i++) {
            for (int j = 0; j < (int)b.size(); j++) c[i + j] += a[i] * b[j];
        }
        return c;
    }
    int n = 1;
    while (n < need) n *= 2;
    a.resize(n);
    b.resize(n);
    ntt(a);
    ntt(b);
    for (int i = 0; i < n; i++) a[i] *= b[i];
    ntt(a, true);
    a.resize(need);
    return a;
}
 
// exact convolution through three NTT primes + CRT (Garner): no modulus restriction.
constexpr unsigned NTT_P1 = 754974721, NTT_P2 = 167772161, NTT_P3 = 469762049;   // product ~ 2^85.6
 
template <unsigned P> vector<MInt<P>> toMInt(const vector<i64>& a) {
    vector<MInt<P>> r(a.size());
    for (int i = 0; i < (int)a.size(); i++) r[i] = MInt<P>(a[i]);
    return r;
}
 
vector<u128> convolutionCRT(const vector<i64>& a, const vector<i64>& b) {   // c_k mod P1 P2 P3
    if (a.empty() || b.empty()) return {};
    auto c1 = convolution(toMInt<NTT_P1>(a), toMInt<NTT_P1>(b));
    auto c2 = convolution(toMInt<NTT_P2>(a), toMInt<NTT_P2>(b));
    auto c3 = convolution(toMInt<NTT_P3>(a), toMInt<NTT_P3>(b));
    constexpr u64 M12 = (u64)NTT_P1 * NTT_P2;
    const MInt<NTT_P2> inv1 = MInt<NTT_P2>(NTT_P1).inv();
    const MInt<NTT_P3> inv12 = MInt<NTT_P3>((i64)(M12 % NTT_P3)).inv();
    vector<u128> res(c1.size());
    for (int i = 0; i < (int)c1.size(); i++) {
        u64 x1 = c1[i].val();
        u64 t2 = ((c2[i] - MInt<NTT_P2>((i64)x1)) * inv1).val();
        u64 x12 = x1 + t2 * NTT_P1;   // c mod P1 P2
        u64 t3 = ((c3[i] - MInt<NTT_P3>((i64)(x12 % NTT_P3))) * inv12).val();
        res[i] = x12 + (u128)t3 * M12;
    }
    return res;
}
 
vector<i64> convolutionLL(const vector<i64>& a, const vector<i64>& b) {   // exact if every |c_k| < 2^63
    const u128 M = (u128)NTT_P1 * NTT_P2 * NTT_P3;
    vector<i64> res;
    for (u128 x : convolutionCRT(a, b)) res.push_back(x > M / 2 ? (i64)((i128)x - (i128)M) : (i64)x);
    return res;
}
 
vector<i64> convolutionAnyMod(vector<i64> a, vector<i64> b, i64 m) {   // any m < 2^31 (prime or not)
    for (auto& x : a) x = (x % m + m) % m;
    for (auto& x : b) x = (x % m + m) % m;
    vector<i64> res;
    for (u128 x : convolutionCRT(a, b)) res.push_back((i64)(x % m));
    return res;
}
 
// ---------- fast factorial mod a prime (n up to 1e9, many queries) ----------
template <class Z>
struct FastFactorial {   // n! mod P for n >= 0, P = Z::mod() a prime < 2^31 (3-prime NTT inside, any prime)
    i64 P;
    int B;
    vector<Z> block;   // block[k] = (k B)!
 
    static vector<Z> conv(const vector<Z>& a, const vector<Z>& b) {
        vector<i64> x(a.size()), y(b.size());
        for (int i = 0; i < (int)a.size(); i++) x[i] = a[i].val();
        for (int i = 0; i < (int)b.size(); i++) y[i] = b[i].val();
        vector<i64> c = convolutionAnyMod(x, y, Z::mod());
        return vector<Z>(c.begin(), c.end());
    }
 
    // h(n), ..., h(n + cnt - 1) from h(0), ..., h(n - 1) where deg h < n (needs n + cnt < P)
    static vector<Z> extend(const vector<Z>& ys, int cnt) {
        int n = ys.size(), len = n + cnt;
        vector<Z> inv(len + 1, Z(1)), ifact(n, Z(1));
        for (int i = 2; i <= len; i++) inv[i] = -Z(Z::mod() / i) * inv[Z::mod() % i];
        for (int i = 1; i < n; i++) ifact[i] = ifact[i - 1] * inv[i];
        vector<Z> c(n), h(len - 1);
        for (int i = 0; i < n; i++) {
            c[i] = ys[i] * ifact[i] * ifact[n - 1 - i];
            if ((n - 1 - i) % 2) c[i] = -c[i];
        }
        for (int s = 0; s < len - 1; s++) h[s] = inv[s + 1];   // 1 / (x - n + 1 + s) at x = n
        vector<Z> cv = conv(c, h), res(cnt);
        Z prod = 1;   // prod_{j<n} (n + t - j)
        for (int j = 1; j <= n; j++) prod *= Z(j);
        for (int t = 0; t < cnt; t++) {
            res[t] = prod * cv[t + n - 1];
            prod *= Z(n + t + 1) * inv[t + 1];
        }
        return res;
    }
 
    FastFactorial(int blockSize = 1 << 11) : P(Z::mod()), B(blockSize) {
        int cnt = (P - 1) / B + 1;   // blocks 0..cnt-1 cover [0, P)
        vector<Z> h(B + 1);          // h(i) = prod_{j=1..B} (B i + j), a polynomial of degree B in i
        for (int i = 0; i <= B; i++) {
            Z p = 1;
            for (int j = 1; j <= B; j++) p *= Z((i64)B * i + j);
            h[i] = p;
        }
        if (cnt > B + 1) {
            vector<Z> more = extend(h, cnt - (B + 1));
            h.insert(h.end(), more.begin(), more.end());
        }
        block.assign(cnt + 1, Z(1));
        for (int k = 0; k < cnt; k++) block[k + 1] = block[k] * h[k];
    }
 
    Z operator()(i64 n) const {
        if (n >= P) return Z(0);
        i64 k = n / B;
        if (n - k * B <= B / 2 || (k + 1) * B >= P) {   // walk up from the block start
            Z r = block[k];
            for (i64 i = k * B + 1; i <= n; i++) r *= Z(i);
            return r;
        }
        Z d = 1;   // walk down from the next block start
        for (i64 i = n + 1; i <= (k + 1) * B; i++) d *= Z(i);
        return block[k + 1] / d;
    }
 
    Z binom(i64 n, i64 k) const {   // n < P
        if (k < 0 || k > n) return Z(0);
        return (*this)(n) / ((*this)(k) * (*this)(n - k));
    }
};

Proof sketch: sampling-point shift

A polynomial $h$ of degree $d$ is determined by $h(0),\dots ,h(d)$; Lagrange’s formula $h(m+t)={\sum }_{i=0}^{d}h(i){\prod }_{j\ne i}\frac{m+t-j}{i-j}$ factors as $({\prod }_j(m+t-j)){\sum }_i\frac{h(i)}{i!(d-i)!(-1{)}^{d-i}}\cdot \frac{1}{m+t-i}$, a convolution of the sequence $\frac{h(i)(-1{)}^{d-i}}{i!(d-i)!}$ with $\frac{1}{m-d+s}$. So one convolution produces $h(m),\dots ,h(m+d)$ from $h(0),\dots ,h(d)$; doubling the degree and the sample range repeatedly (or one large shift) gives the block products for all $i<p/B$. Multiplying the block products in order yields $(kB)!$, and the remainder of a query is at most $B$ direct multiplications. All arithmetic is in ${\mathbb{F}}_p$; the 3-prime NTT with CRT computes the integer convolution exactly before reduction.

Binomials modulo any $m$

For each $p^e\Vert m$: $n!=p^{v_p(n!)}F(n)$ where $F(n)$ is the product of the numbers $\le n$ coprime to $p$, computed as $F(n)=U^{\lfloor n/p^e\rfloor }\cdot \text{units}(n\mathrm{mod}{p}^e)\cdot F(\lfloor n/p\rfloor )$ with $U=\text{units}(p^e-1)=\pm 1$. Then $\left(\genfrac{}{}{0ex}{}{n}{k}\right)\equiv p^{v(n)-v(k)-v(n-k)}F(n)F(k{)}^{-1}F(n-k{)}^{-1}(\mathrm{mod}{p}^e)$, and CRT combines the prime powers.

// usage: BinomialAnyMod binom(m); binom.C(n, k)   (n, k < 2^63)
// ---------- extended Euclid and linear Diophantine equations (number-theory/extended-gcd.cpp) ----------
i64 extGcd(i64 a, i64 b, i64& x, i64& y) {
    i64 x1 = 0, y1 = 1;
    x = 1, y = 0;
    while (b != 0) {
        i64 q = a / b, t;
        t = a - q * b, a = b, b = t;
        t = x - q * x1, x = x1, x1 = t;
        t = y - q * y1, y = y1, y1 = t;
    }
    if (a < 0) a = -a, x = -x, y = -y;
    return a;
}
 
i64 invMod(i64 a, i64 m) {
    i64 x, y;
    a %= m;
    if (a < 0) a += m;
    if (extGcd(a, m, x, y) != 1) return -1;
    x %= m;
    return x < 0 ? x + m : x;
}
 
bool diophantine(i64 a, i64 b, i64 c, i64& x, i64& y, i64& g) {
    if (a == 0 && b == 0) {
        x = y = g = 0;
        return c == 0;
    }
    i64 x1, y1;
    g = extGcd(a, b, x1, y1);
    if (c % g != 0) return false;
    if (b == 0) {
        x = c / a, y = 0;
        return true;
    }
    i64 step = abs(b / g);
    i128 xx = (i128)x1 * (c / g) % step;
    if (xx < 0) xx += step;
    x = (i64)xx;
    y = (i64)(((i128)c - (i128)a * x) / b);
    return true;
}
 
i64 countDiophantine(i64 a, i64 b, i64 c, i64 xlo, i64 xhi, i64 ylo, i64 yhi) {
    i64 x0, y0, g;
    if (xlo > xhi || ylo > yhi || !diophantine(a, b, c, x0, y0, g)) return 0;
    if (a == 0 && b == 0) return (xhi - xlo + 1) * (yhi - ylo + 1);
    auto floorDivide = [](i128 p, i128 q) { return p >= 0 ? p / q : -((-p + q - 1) / q); };
    i128 lo = -((i128)1 << 100), hi = (i128)1 << 100;   // range of t
    // restrict t so that vlo <= v0 + t * d <= vhi.
    auto clip = [&](i128 v0, i128 d, i128 vlo, i128 vhi) {
        if (d == 0) {
            if (v0 < vlo || v0 > vhi) lo = 1, hi = 0;
            return;
        }
        if (d < 0) d = -d, v0 = -v0, swap(vlo, vhi), vlo = -vlo, vhi = -vhi;
        lo = max(lo, -floorDivide(v0 - vlo, d));
        hi = min(hi, floorDivide(vhi - v0, d));
    };
    clip(x0, b / g, xlo, xhi);
    clip(y0, -(a / g), ylo, yhi);
    return hi < lo ? 0 : (i64)(hi - lo + 1);
}
 
// ---------- Chinese remainder theorem (number-theory/crt.cpp) ----------
pair<i64, i64> crt(i64 r1, i64 m1, i64 r2, i64 m2) {
    r1 %= m1;
    if (r1 < 0) r1 += m1;
    r2 %= m2;
    if (r2 < 0) r2 += m2;
    i64 p, q;
    i64 g = extGcd(m1, m2, p, q);
    if ((r2 - r1) % g != 0) return {-1, -1};
    i64 step = m2 / g;
    // x = r1 + m1 * t with (m1 / g) * t = (r2 - r1) / g (mod m2 / g)
    i128 t = (i128)((r2 - r1) / g) * p % step;
    if (t < 0) t += step;
    return {(i64)(r1 + (i128)m1 * t), m1 / g * m2};
}
 
pair<i64, i64> crt(const vector<i64>& r, const vector<i64>& m) {
    pair<i64, i64> cur = {0, 1};
    for (int i = 0; i < (int)r.size(); i++) {
        cur = crt(cur.first, cur.second, r[i], m[i]);
        if (cur.second == -1) break;
    }
    return cur;
}
 
// Garner's algorithm: the solution in [0, prod m) reduced modulo mod. moduli must be pairwise coprime.
i64 crtMod(const vector<i64>& r, const vector<i64>& m, i64 mod) {
    int k = m.size();
    vector<i64> ms = m;
    ms.push_back(mod);
    vector<i64> coef(k + 1, 1), val(k + 1, 0);   // coef[j] = m[0] * ... * m[i-1] (mod ms[j])
    for (int i = 0; i < k; i++) {
        i64 t = (r[i] % m[i] - val[i]) % m[i];
        if (t < 0) t += m[i];
        t = (i64)((i128)t * invMod(coef[i], m[i]) % m[i]);
        for (int j = i + 1; j <= k; j++) {
            val[j] = (i64)((val[j] + (i128)t * coef[j]) % ms[j]);
            coef[j] = (i64)((i128)coef[j] * m[i] % ms[j]);
        }
    }
    return val[k];
}
 
// ---------- binomial coefficients modulo any integer ----------
struct BinomialAnyMod {
    struct PrimePower {
        i64 p, q;             // q = p^e
        int e;
        vector<i64> units;    // units[i] = product of 1 <= j <= i with p not dividing j, modulo q
    };
    vector<PrimePower> parts;
 
    BinomialAnyMod(i64 m) {
        for (i64 p = 2; m > 1; p++) {
            if (p * p > m) p = m;   // what remains of m is a prime
            if (m % p != 0) continue;
            PrimePower t{p, 1, 0, {}};
            while (m % p == 0) m /= p, t.q *= p, t.e++;
            t.units.assign(t.q, 1);
            for (i64 i = 1; i < t.q; i++) t.units[i] = i % p == 0 ? t.units[i - 1] : t.units[i - 1] * i % t.q;
            parts.push_back(move(t));
        }
    }
 
    // n! with every factor p removed, modulo q.
    static i64 factorialUnits(const PrimePower& t, i64 n) {
        i64 res = 1;
        for (; n > 0; n /= t.p) {
            if ((n / t.q) % 2 == 1) res = res * t.units[t.q - 1] % t.q;
            res = res * t.units[n % t.q] % t.q;
        }
        return res;
    }
 
    // exponent of p in n! (Legendre's formula).
    static i64 factorialExponent(i64 n, i64 p) {
        i64 e = 0;
        while (n > 0) n /= p, e += n;
        return e;
    }
 
    static i64 modPrimePower(const PrimePower& t, i64 n, i64 k) {
        i64 e = factorialExponent(n, t.p) - factorialExponent(k, t.p) - factorialExponent(n - k, t.p);
        if (e >= t.e) return 0;
        i64 res = factorialUnits(t, n) * invMod(factorialUnits(t, k), t.q) % t.q;
        res = res * invMod(factorialUnits(t, n - k), t.q) % t.q;
        while (e-- > 0) res = res * t.p % t.q;
        return res;
    }
 
    i64 C(i64 n, i64 k) const {
        if (k < 0 || k > n) return 0;
        vector<i64> r, mods;
        for (const PrimePower& t : parts) {
            r.push_back(modPrimePower(t, n, k));
            mods.push_back(t.q);
        }
        return crt(r, mods).first;
    }
};

Proof: the prime-power factorial recursion (Granville)

Among $\mathrm{1..}n$, the numbers coprime to $p$ come in $\lfloor n/p^e\rfloor$ full periods of the residues coprime to $p$ modulo $p^e$ (each period’s product is $U$) plus a partial period, giving $U^{\lfloor n/p^e\rfloor }\text{units}(n\mathrm{mod}{p}^e)$. The multiples of $p$ are $p\cdot 1,\dots ,p\lfloor n/p\rfloor$, contributing $p^{\lfloor n/p\rfloor }\lfloor n/p\rfloor !$, whose $p$-free part is $F(\lfloor n/p\rfloor )$ (recursively). $U\equiv -1$ except for $p=2,e\ge 3$ where it is $+1$ (the group of units modulo $p^e$ is cyclic for odd $p$ and for $4$, and ${\mathbb{Z}}_2\times {\mathbb{Z}}_{2^{e-2}}$ otherwise; the product of all elements of a finite abelian group is the product of its elements of order 2). The $p$-parts give the exponent $v(n)-v(k)-v(n-k)\ge 0$ (Kummer), and the $p$-free parts are units, so the quotient is well defined modulo $p^e$; CRT reassembles $m$.