二項係数(mod p^q)
(math/binom_mod_prime_power.cpp)

• View this file on GitHub
• Last update: 2026-03-08 22:25:54+09:00

説明

素数冪 p^q を法とする二項係数 C(n, k) を扱う。 前計算 $O(p^q)$、各クエリ $O(log_p n)$。

できること

• BinomModPrimePower(ll p, int q) 法 p^q 用のテーブルを作る
• ll C(ll n, ll k) C(n, k) mod p^q を返す。k < 0 または k > n なら 0
• ll modulus() 法 p^q を返す

使い方

BinomModPrimePower binom(p, q); を作り、binom.C(n, k) を呼ぶ。 一般の合成数 mod m で使いたいときは、m を素因数分解して各 p^q ごとに値を出し、CRT で併合する。

実装上の補足

n! の p を除いた部分と p の指数を再帰的に分けて計算する。 メモリは $O(p^q)$ 使うので、法が大きすぎる場合には向かない。

Depends on

• math/modinv.cpp

Verified with

• test/yosupo_binomial_coefficient.test.cpp

Code

#include "modinv.cpp"

struct BinomModPrimePower {
    ll p, mod;
    int q;
    ll block_prod;
    vector<ll> ppow;
    vector<int> prod;

    explicit BinomModPrimePower(ll prime, int exponent) : p(prime), mod(1), q(exponent), ppow(exponent + 1, 1) {
        for (int i = 0; i < q; ++i) {
            mod *= p;
            ppow[i + 1] = mod;
        }
        block_prod = (p == 2 && q >= 3 ? 1 : mod - 1);
        prod.assign(mod + 1, 1);
        for (int i = 1; i <= mod; ++i) {
            prod[i] = prod[i - 1];
            if (i % p != 0) prod[i] = (ull)prod[i] * i % mod;
        }
    }

    pair<ll, ll> factorial(ll n) const {
        ll x = 1, e = 0;
        while (n) {
            if (block_prod != 1 && (n / mod) & 1) x = mod - x;
            x = (ull)x * prod[n % mod] % mod;
            n /= p;
            e += n;
        }
        return {x, e};
    }

    ll C(ll n, ll k) const {
        if (k < 0 || k > n) return 0;
        auto [a, ea] = factorial(n);
        auto [b, eb] = factorial(k);
        auto [c, ec] = factorial(n - k);
        ll e = ea - eb - ec;
        if (e >= q) return 0;
        ll x = (ull)b * c % mod;
        return (ull)a * mod_inv(x, mod) % mod * ppow[e] % mod;
    }

    ll modulus() const {
        return mod;
    }
};

/**
 * @brief 二項係数(mod p^q)
 */

#line 1 "math/modinv.cpp"
template<typename T>  
T mod_inv(T x, T M){  
   T u = 1, t = 1, v = 0, s = 0, m = M;  
   while (x) { T q = m/x; swap(s -= q*u, u); swap(t -= q*v, v); swap(m -= q*x, x); }  
   if(s < 0) s += M;  
   return s;  
}
#line 2 "math/binom_mod_prime_power.cpp"

struct BinomModPrimePower {
    ll p, mod;
    int q;
    ll block_prod;
    vector<ll> ppow;
    vector<int> prod;

    explicit BinomModPrimePower(ll prime, int exponent) : p(prime), mod(1), q(exponent), ppow(exponent + 1, 1) {
        for (int i = 0; i < q; ++i) {
            mod *= p;
            ppow[i + 1] = mod;
        }
        block_prod = (p == 2 && q >= 3 ? 1 : mod - 1);
        prod.assign(mod + 1, 1);
        for (int i = 1; i <= mod; ++i) {
            prod[i] = prod[i - 1];
            if (i % p != 0) prod[i] = (ull)prod[i] * i % mod;
        }
    }

    pair<ll, ll> factorial(ll n) const {
        ll x = 1, e = 0;
        while (n) {
            if (block_prod != 1 && (n / mod) & 1) x = mod - x;
            x = (ull)x * prod[n % mod] % mod;
            n /= p;
            e += n;
        }
        return {x, e};
    }

    ll C(ll n, ll k) const {
        if (k < 0 || k > n) return 0;
        auto [a, ea] = factorial(n);
        auto [b, eb] = factorial(k);
        auto [c, ec] = factorial(n - k);
        ll e = ea - eb - ec;
        if (e >= q) return 0;
        ll x = (ull)b * c % mod;
        return (ull)a * mod_inv(x, mod) % mod * ppow[e] % mod;
    }

    ll modulus() const {
        return mod;
    }
};

/**
 * @brief 二項係数(mod p^q)
 */

Back to top page