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

数学剰余, 二項係数

競技プログラミングの問題などでは、二項係数を非常に大きい素数 P で割った余りを出力させる問題が出題されることがあります。 \(P = 1000000007 = 10^9 + 7\) の素数を使用することなどが多いです。

単純に計算すると、\(\frac{n!}{k!(n-k)!}\)を計算してから割ることになりますが、値が非常に大きくなってしまいます。32bitや64bitの整数には収まり切らないかもしれません。

また、上手く行ったとしても、1つの二項係数を計算するのに\(O(n)\)だけかかってしまいます。

しかし、始めに\(O(n)\)の前計算で「mod P」での階乗とその逆元を求めておくと、1つの二項係数を素数 P で割った余りを定数時間で求められるようになります。

※ 以下では「mod P」での \(a\) の逆元を \(a^{-1}\) などと書くことにします。

※ 逆元や数学的証明について詳しくは後述

アルゴリズム

\(_nC_k\) におけるn や k のサイズ・二項係数を何回計算するのか、などによって最適なアルゴリズムが変わってきますが、ここではよく使われる手法について説明します。

\(_nC_k\) を素数 P で割った余りを求めるためには、i の階乗を P で割った余りを fact[i]、その逆元を fact_inv[i]、i の逆元を inv[i] とした時に、以下のように計算します。

  1. i が 0 から n のときまで、fact[i], fact_inv[i], inv[i] を計算して結果を保存しておく(メモ化)
  2. \(_nC_k\) を素数 P で割った余りは、「 fact[n]・ fact_inv[k]・fact_inv[n-k] 」を P で割った余りとなる

※「fact[n]・ fact_inv[k]・fact_inv[n-k]」を P で割った余りを計算する際は、掛け算を1回行うごとに P で割った余りを求めないと、32bit や 64 bit の整数をオーバーしてしまうかもしれないので注意が必要です。

階乗と逆元の求め方

fact[i] , fact_inv[i], inv[i] はメモ化や動的計画法と呼ばれる方法で、高速に計算することができます。i 番目の計算結果を、(i+1)番目の計算のときに利用することで、早い計算が可能になるのです。

数式では以下のようになります(説明は後述)。

  • fact[i+1]:
    • \((i+1)! \equiv i! \times (i+1) ~(mod ~P)\)
  • inv[i+1]:
    • \((i+1)^{-1} \equiv ~ -(P\%(i+1))^{-1} \times [\frac{P}{i+1}] ~(mod ~P)\)
  • fact_inv[i+1]:
    • \(((i+1)!)^{-1} \equiv (i!)^{-1} \times (i+1)^{-1} ~(mod ~P)\)

※ \([\frac{P}{i+1}] \)は P を i+1 で割った数の小数点以下を切り捨てた数です。

よって以下のようにプログラムを書けば良いです。

  • fact[i+1] = fact[i] * (i+1) % P
  • inv[i+1] = P – inv[P%(i+1)] * (P / (i+1)) % P
  • fact_inv[i+1] = fact_inv[i] * inv[i+1] % P

計算量

階乗と逆元を1つずつ求めるのに \(O(1)\) かかるので、n まで求めると\(O(n)\)かかります。

\(_nC_k\) は、階乗と逆元の値さえわかっていれば\(O(1)\)で計算できます。仮に n 個の二項係数を求めたとしても \(O(n)\) です。

合わせて計算量は \(O(n)\) となります。ただし、n が非常に大きいと階乗と逆元を覚えておくためのメモリが足りなくなってしまうので注意が必要です。

数学的な説明と原理

そもそも逆元とは

逆元とは、ある要素の対になるもので、足し算や掛け算をすると「お互いに打ち消し合う」ようなものです。

例えば、通常の掛け算においては、\(3\) の逆元は \(\frac{1}{3}\)となります。掛け算をすると打ち消し合って \(1\) になります。

$$3 \times \frac{1}{3} = 1 $$

同様に、素数Pで割った余りを考える「mod P の世界」でも逆元を考えることができます。 \(a\) の逆元を \(a^{-1}\) のように書くことにすると以下のようになります。

$$ a \times a^{-1} \equiv 1 ~(mod ~ P) $$

例えば、

$$2\times 3 \equiv 1~(mod ~ 5)$$

という式が成り立つので、「mod 5 の世界」においては、 2 の逆元が 3 ということになります。

二項係数がなぜ逆元を使って求まるのか

なぜ上述のアルゴリズムで、二項係数を素数で割った余りが求められるのかを説明します。

\(_nC_k\)を素数 P で割った余りが、fact[n]・fact_inv[k]・fact_inv[n-k]であること、つまり、

$$_nC_k \equiv n! \cdot (k!)^{-1}\cdot ((n-k)!)^{-1} ~ (mod ~ P)$$

となる理由を考えてみましょう。

まず二項係数の定義から以下のような式が成り立ちます。

$$_nC_k = \frac{n!}{k!(n-k)!} = n!\cdot \frac{1}{k!} \cdot \frac{1}{(n-k)!}$$

ここで、「mod P の世界」を考えてみます。すると、

  • mod P における \(n!\)
  • mod P における \(\frac{1}{k!}\) つまり \((k!)^{-1} \)
  • mod P における \(\frac{1}{(n-k)!}\) つまり \(((n-k)!)^{-1} \)

の3つが分かれば、あとは掛け算を使用して計算をすることができます。1つ目の計算は簡単ですが、2つ目と3つ目はmodでの割り算になってしまいます。なぜ逆元使って表すことができるのかを説明しましょう。

modでの割り算と逆元の関係

  • mod P における \(\frac{1}{k!}\)
  • mod P における \(\frac{1}{(n-k)!}\)

を計算するために「mod P の世界」における割り算とは何かを考えてみましょう。

例えば、$$2\times 3 \equiv 1~(mod ~ 5)$$という式を考えてみます。両辺を\(3\)で割ることを考えると、$$2\equiv \frac{1}{3}~(mod ~ 5)$$となるはずです。

同様にmod P における \(\frac{1}{i}\) の値を \(a\) とすると、

$$ \frac{1}{i} \equiv a ~(mod ~P)$$

$$ 1 \equiv a\cdot i ~(mod ~P)$$

となるので、実は「mod P の世界」における \(\frac{1}{i}\) の値 \(a\) は、 「mod P の世界」における \(i\) の逆元 \(i^{-1}\) となります。

つまり

$$ \frac{1}{i} \equiv (i)^{-1} $$

以上より、\(\frac{1}{k!}\) については以下の式が成立します。

$$ \frac{1}{k!} \equiv (k!)^{-1} $$

逆元の高速な計算方法

  • 拡張ユークリッドの互除法
  • フェルマーの小定理と繰り返し2乗法

などの方法を使うと、i の逆元 を \(O(log P)\) で求めることができますが、0からnまでの逆元を求める場合では \(O(n log P)\) だけかかってしまいます。

ですが、以下のように「i の逆元を、P % i の逆元を利用して求める」ことで、1つあたりO(1)で求めることが可能になります。

$$i^{-1} \equiv ~ -(P\%i)^{-1} \times [\frac{P}{i}] ~(mod ~P)$$

これを示してみましょう。

まずP を i で割った商は\([\frac{P}{i}]\) なので、

$$ P = [\frac{P}{i}] \times i + (P\%i) $$

などとなります。mod P で考えると、

$$ 0 \equiv [\frac{P}{i}] \times i + (P\%i) ~(mod ~P) $$

全体に \(i^{-1}\) を掛けて変形すると

$$ 0 \equiv [\frac{P}{i}] + (P\%i) \times i^{-1} ~(mod ~P)$$

$$ (P\%i) \times i^{-1} \equiv ~ -[\frac{P}{i}] ~(mod ~P)$$

$$ i^{-1} \equiv ~ -(P\%i)^{-1} \times [\frac{P}{i}] ~(mod ~P) $$

となって示されました。あとはこれをプログラムで書き直すと、

  • inv[i] = P – inv[P%i] * (P / i) % P

のようになります。

プログラム例

C++ の例

#include <bits/stdc++.h>
using namespace std;
const int MAX = 1000000;
const int MOD = 1e9 + 7;
long long fact[MAX], fact_inv[MAX], inv[MAX];
void pre_calculation_combi() { // fact[], fact_inv[], inv[]を前計算
    fact[0] = fact[1] = 1;
    fact_inv[0] = fact_inv[1] = 1;
    inv[1] = 1;
    for (int i = 1; i < MAX - 1; i++) {
        fact[i + 1] = fact[i] * (i + 1) % MOD;
        inv[i + 1] = MOD - inv[MOD % (i + 1)] * (MOD / (i + 1)) % MOD;
        fact_inv[i + 1] = fact_inv[i] * inv[i + 1] % MOD;
    }
}
long long combination_mod(int n, int k) { // fact[] and fact_inv[] で nCk % MOD を求める
    if (n < k) {
        return 0;
    } else if (n < 0 || k < 0) {
        return 0;
    } else {
        return ((fact[n] * fact_inv[k]) % MOD) * fact_inv[n - k] % MOD;
    }
}
int main() {
    pre_calculation_combi();
    cout << combination_mod(10, 1) << endl;
    return 0;
}

関連問題

参考になる記事

数学剰余, 二項係数

Posted by cs-learner