二項係数(nCk)を素数(p)で割った余りの計算(Lucas の定理)

2020年3月22日数学剰余,二項係数,動的計画法,Lucasの定理,素数の剰余

Lucas の定理を利用すると、\(_{n}\mathrm{C}_{k}\)% \(p\) が \(O( p^2 \log_p n)\) で計算できます。素数 \(p\) が小さい場合は十分高速です。

Lucas の定理

任意の素数 \(p\) と、非負整数 \(n, k\) について
$$ _{n}\mathrm{C}_{k} \equiv \prod_{i=0}^m{}_{n_i}\mathrm{C}_{k_i}\:(\mathrm{mod}\:p)$$
が成立する。ただし、
\begin{align}
n &= n_m p^m+ n_{m-1} p^{m-1} + n_{m-2} p^{m-2} + \cdots + n_1 p^1 + n_0 \\
k&= k_m p^m + k_{m-1} p^{m-1} + k_{m-2} p^{m-2} + \cdots + k_1 p^1 + k_0
\end{align}
である。(それぞれ \(p\) 進数で表した時の数。)

※ \(n < k\) の時は \(_{n}\mathrm{C}_{r} = 0\) とします。

アルゴリズム

Lucas の定理を用いた二項係数の剰余の計算:

  • \({}_{n_i}\mathrm{C}_{k_i}\:(\mathrm{mod}\:p)\) の計算法:
    • \({}_{i}\mathrm{C}_{j} = {}_{i-1}\mathrm{C}_{j-1}+{}_{i-1}\mathrm{C}_{j} \) を利用して計算(動的計画法)
  • メイン:
    • \(_{n}\mathrm{C}_{k} \equiv \prod_{i=0}^m{}_{n_i}\mathrm{C}_{k_i}\:(\mathrm{mod}\:p)\) を計算する

※ \({}_{n_i}\mathrm{C}_{k_i}\) の計算は、パスカルの三角形を上から動的計画法で作っていくイメージです。

※ その他の nCk % p の計算法はこちらを参照してください。

計算量

  • 前処理をしない時:
    • 前処理:なし
    • クエリ:\(O( p^2 \log_p n)\)
  • 前処理をする時:
    • 前処理:\(O(p^2)\)
    • クエリ:\(O(\log_p n)\)

桁数は \(m = O(\log_p n)\) となるので、その回数だけ \({}_{n_i}\mathrm{C}_{k_i}\) の計算をすることになります。

\({}_{n_i}\mathrm{C}_{k_i}\) の計算には \(O(p^2)\) だけかかります(\(n_i, k_i < p\) のため)。計算結果は使い回せるので、前処理して保存しておけば高速に計算ができます。

C++ での実装例

前処理をする場合の実装例です。構造体を用いています。

/* Com:nCk % p の計算のための構造体
    前処理・初期化: O(p^2)
    nCk(n,k): nCk % p の計算。O(log n)
*/
struct Comb {
    vector<vector<long long>> com;  // 前計算の結果を保存
    long long p;                    // p は素数である必要がある
    Comb(long long _p) : p(_p) {
        init(p);
    }
    void init(long long p) {  // 動的計画法で前処理
        com.assign(p, vector<long long>(p));
        com[0][0] = 1;
        for (int i = 1; i < p; i++) {
            com[i][0] = 1;
            for (int j = i; j > 0; j--) {
                com[i][j] = (com[i - 1][j - 1] + com[i - 1][j]) % p;
            }
        }
    }
    long long nCk(long long n, long long k) {
        long long ret = 1;
        while (n > 0) {  // 下から一桁ずつ計算する
            int ni = n % p;
            int ki = k % p;
            ret *= com[ni][ki];
            ret %= p;
            n /= p;
            k /= p;
        }
        return ret;
    }
};