AtCoderのABC042Dで使う、組み合わせの数\({}_n C_r\)とその余りの効率的な求め方についてまとめました。合同式やフェルマーの小定理について書いています。
問題設定
\({}_n C_r = \frac{n!}{(n-r)! r!}\) を求める時には、その数がとても大きくなることが多いため、\( 10^9 + 7\)で割ったあまりの数を使います。
\({}_n C_r\)の計算ではそのまま階乗を使った計算をすると、実行時間が間に合わなくなるため、あらかじめ\(0 \leq x \leq n\)に対して、\(x! (mod 10^9 + 7)\)と\((x!)^{-1} (mod 10^9 + 7)\)を計算してテーブルにしておくことで、その後の\({}_n C_r\) の計算に時間をかけないようにします。
では、その\(x! (mod 10^9 + 7)\)と\((x!)^{-1} (mod 10^9 + 7)\)をどのように計算していくか解説していきます。
\(x! (mod 10^9 + 7)\)のテーブル作成
\(0 \leq x \leq n\)に対して\(x! (mod 10^9 + 7)\)のテーブル作成には、合同式の性質を利用します。
合同式では
$$ a \equiv c, b \equiv d (mod p)$$
ならば、
$$ ab \equiv cd (mod p)$$
が成立します。
よって、
\( 0! \equiv 1 (mod 10^9 +7)\)の両辺に1をかけて、 \( 1! \equiv 1 \times 1 \equiv 1 (mod 10^9 +7)\)を求め、それの両辺に2をかけて\( 2! \equiv 1 \times 2 \equiv 2 (mod 10^9 +7)\)を求め、というのを繰り返して\(x! (mod 10^9 + 7)\)のテーブル作成を行います。
Python3でのコードの例は下の通りです。
MOD = 10**9 + 7
table = [1]
for i in range(1, N+1):
table.append((table[-1]*i) % MOD)\(x\)の逆元\(x^{-1}\)とは？
\( x ^{-1} (mod 10^9 + 7)\)とは何か。
それは、\( x x^{-1} = 1\)より、\(x (mod10^9 + 7)\)にかけて\(mod 10^9 + 7\)で1となる数です。\(x^{-1}\)は\( x \equiv 0 (mod 10^9 +7)\)でない任意の\(x\)に対して存在します。
その\(x^{-1}\)を求めるために、フェルマーの小定理
” \(p\)が素数で、整数\(a\)が\(p\)と互いに素であるとき、\(a^{p-1} \equiv 1 (mod p)\) “
を利用します。
\(10^9 + 7\)は素数なので、
$$a^{10^9 + 7 – 1} = a^{10^9 + 6} = a a^{10^9 + 5}\equiv 1 (mod 10^9+7) $$
より、\(a^{-1} = a^{10^9 + 5}\)が導けます。
以上より、\(0 \leq x \leq n\)に対して\((x!)^{-1}\)を求める時には、\((n!)^{-1} (mod 10^9 + 7)\)を求めて、それに\(n\)をかけると\(((n-1)!)^{-1} (mod 10^9 + 7)\)が求められ、それに\(n-1\)をかけると\(((n-2)!)^{-1} (mod 10^9 + 7)\)が求まるというように、高速に逆元のテーブル作成ができます。
Python3でのコードの例は下の通りです。途中までは上のソースコードと同じです。
MOD = 10**9 + 7
table = [1]
for i in range(1, N+1):
table.append((table[-1]*i) % MOD)
inverse_table = []
inverse_table.append(pow(table[-1], MOD-2, MOD)) % N!(mod MOD) のMOD-2乗のMODでの余り
for i in range(N):
inverse_table.insert(0, (inverse_table[0]*(N-i)) % MOD)
inverse_table.insert(0, 1)また、簡単な書き方として\(a^{-1} = a^{10^9 + 5}\)をそれぞれの\(a\)で適応する以下のようなコードも考えられます。
MOD = 10**9 + 7
table = [1]
inverse_table = [1]
for i in range(1, N+1):
table.append((table[-1]*i) % MOD)
inverse_table.append(pow(table[-1], MOD-2, MOD))