競プロでよく使う二項係数(nCk)を素数(p)で割った余りの計算と逆元のまとめ

2019年12月4日数学剰余,二項係数,テーマ記事,競プロ,素数,逆元

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

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

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

この記事ではよく使われる二項係数の計算方法をいくつかご紹介しつつ、数学的な理解をすることを目指します。

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

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

様々な二項係数のアルゴリズム

\(_nC_k\) におけるn や k のサイズ、二項係数を何回計算するのか、などによって最適なアルゴリズムが変わってきます。

\(k \leq n \leq 10^7\) でPが素数のとき

まずは一番よく使われる計算方法について説明します。（以下では P は n より大きいものとしています）

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

\(k \leq n \leq 10^7\) でPが素数のとき：

  • 前処理
    1. i が 0 から n のときまで、fact[i], fact_inv[i], inv[i] を計算して結果を保存しておく(動的計画法)
  • クエリ
    1. \(_nC_k\) を素数 P で割った余りは、「 fact[n]･ fact_inv[k]･fact_inv[n-k] 」を P で割った余りとなる

※ 始めに\(O(n)\)の前計算で「mod P」での階乗とその逆元を求めておくと、1つの二項係数を素数 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 で割った数の小数点以下を切り捨てた数です。

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

for (int i = 2; i < 10000000; i++) {
    fact[i] = fact[i - 1] * i % P;
    inv[i] = P - inv[P % i] * (P / i) % P;
    fact_inv[i] = fact_inv[i - 1] * inv[i] % P;
}

計算量

  • 前処理：\(O(n)\)
  • クエリ：\(O(1)\)

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

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

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

C++ での実装例

#include <bits/stdc++.h>
using namespace std;
const int MOD = 1e9 + 7;

vector<long long> fact, fact_inv, inv;

/*  init_nCk :二項係数のための前処理
    計算量:O(n)
*/
void init_nCk(int SIZE) {
    fact.resize(SIZE + 5);
    fact_inv.resize(SIZE + 5);
    inv.resize(SIZE + 5);
    fact[0] = fact[1] = 1;
    fact_inv[0] = fact_inv[1] = 1;
    inv[1] = 1;
    for (int i = 2; i < SIZE + 5; i++) {
        fact[i] = fact[i - 1] * i % MOD;
        inv[i] = MOD - inv[MOD % i] * (MOD / i) % MOD;
        fact_inv[i] = fact_inv[i - 1] * inv[i] % MOD;
    }
}

/*  nCk :MODでの二項係数を求める(前処理 int_nCk が必要)
    計算量:O(1)
*/
long long nCk(int n, int k) {
    assert(!(n < k));
    assert(!(n < 0 || k < 0));
    return fact[n] * (fact_inv[k] * fact_inv[n - k] % MOD) % MOD;
}

n が巨大 (\(k \leq 10^7, n \leq 10^9\)) でPが素数のとき

n が非常に大きい場合は、先程のような \(O(n)\) の前処理計算が時間内にできません。そこで以下のように\(O(k)\) の計算をします。

n が巨大 (\(k \leq 10^7, n \leq 10^9\)) でPが素数のとき：

  • 前処理
    1. i が 0 から k のときまで、fact_inv[i], inv[i] を計算して結果を保存しておく(動的計画法)
  • クエリ
    1. \(_nC_k = \frac{n}{1} × \frac{n-1}{2} × … × \frac{n-k+1}{k} \) と fact_inv[k], inv[i] を利用して計算

fact_inv[i], inv[i] を求めるアルゴリズムは前述の通りです。この記事の後半で数学的に詳しい解説をします。

\(_nC_k = \frac{n}{1} × \frac{n-1}{2} × … × \frac{n-k+1}{k} \) の計算

階上の逆元である

  • fact_inv[ k ] := \((k!)^{-1}\)

を利用することで、\(_nC_k\) は以下のように \(O(k)\) で計算することができます。

  • \(\frac{n}{1} × \frac{n-1}{2} × … × \frac{n-k+1}{k} = n(n-1)(n-2)…(n-k+1)\) × fact_inv[k]

n が固定値の場合は

  • Com[ k ] := \(_nC_k\) の値

などとすることで、\(O(k)\) で以下のように前処理計算をすることもできます。

  • Com[ i ] = Com[ i-1 ] × \(\frac{n-i+1}{i}\)

計算量

nが固定値でない場合は

  • 前処理：\(O(k)\)
  • クエリ：\(O(k)\)

nが固定値の場合は

  • 前処理：\(O(k)\)
  • クエリ：\(O(1)\)

で計算することができます。

\(k \leq 10^5\) 程度なら簡単に実装できる

k の値がもう少し小さければ、前処理計算無しでさらに簡単に計算することができます。こちらの方が理解はしやすいでしょう。

$$_nC_k = \frac{n}{1} × \frac{n-1}{2} × … × \frac{n-k+1}{k}$$

を直接計算することを考えれば良いです。

「i の割り算」は P で割った余りを考えると「i の逆元の掛け算」のことであり、これは 1つあたり \(O(log( i ))\) で計算することができます。

$$\frac{n}{1} × \frac{n-1}{2} × … × \frac{n-k+1}{k} = n(1)^{-1} × (n-1)(2)^{-1} × … × (n-k+1)(k)^{-1}$$

を直接計算すれば良いことになります。

C++ での実装例

まずは n が固定値ではない場合です。

#include <bits/stdc++.h>
using namespace std;
const int MOD = 1e9 + 7;

vector<long long> fact_inv, inv;

/*  init_nCk :二項係数のための前処理
    計算量:O(k)
*/
void init_nCk(int SIZE) {
    fact_inv.resize(SIZE + 5);
    inv.resize(SIZE + 5);
    fact_inv[0] = fact_inv[1] = 1;
    inv[1] = 1;
    for (int i = 2; i < SIZE + 5; i++) {
        inv[i] = MOD - inv[MOD % i] * (MOD / i) % MOD;
        fact_inv[i] = fact_inv[i - 1] * inv[i] % MOD;
    }
}

/*  nCk :MODでの二項係数を求める(前処理 int_nCk が必要)
    計算量:O(k)
*/
long long nCk(int n, int k) {
    assert(!(n < k));
    assert(!(n < 0 || k < 0));
    long long ans = 1;
    for (int i = n; i >= n - k + 1; i--) {
        ans *= i;
        ans %= MOD;
    }
    return ans * fact_inv[k] % MOD;
}

n が固定値の場合は以下のようになります。

#include <bits/stdc++.h>
using namespace std;
const int MOD = 1e9 + 7;

vector<long long> fact_inv, inv, Com;

/*  init_nCk :二項係数のための前処理
    計算量:O(n)
*/
void init_nCk(int n, int SIZE) {
    fact_inv.resize(SIZE + 5);
    inv.resize(SIZE + 5);
    fact_inv[0] = fact_inv[1] = 1;
    inv[1] = 1;
    for (int i = 2; i < SIZE + 5; i++) {
        inv[i] = MOD - inv[MOD % i] * (MOD / i) % MOD;
        fact_inv[i] = fact_inv[i - 1] * inv[i] % MOD;
    }
    Com.resize(SIZE + 5);
    Com[0] = 1;
    for (int i = 1; i < SIZE + 5; i++) {
        Com[i] = Com[i - 1] * ((n - i + 1) * inv[i] % MOD) % MOD;
    }
}

/*  nCk :MODでの二項係数を求める(前処理 int_nCk が必要)
    計算量:O(1)
*/
long long nCk(int k) {
    assert(!(k < 0));
    return Com[k];
}

\(k \leq n \leq 5000\) でPが素数でないとき

P が素数ではない時も考えられます。アルゴリズムは以下のとおりです。

\(k \leq n \leq 5000\) でPが素数でないとき：

  • 前処理
    1. \(_nC_k = {}_{n-1}C_{k-1} +{} _{n-1}C_{k}\) （動的計画法）
  • クエリ
    1. 先程の計算結果を利用

※ もちろんPが素数の時もこのアルゴリズムは有効ですが、その場合は前述のアルゴリズムのほうが高速です。

動的計画法での前処理

二項係数の有名公式

$$_nC_k = {}_{n-1}C_{k-1} +{} _{n-1}C_{k}$$

を利用して前処理をしておきます。

以下のような配列を定義しましょう。

  • Com[ n ][ k ] := \(_nC_k\) の値

動的計画法での計算は以下の通りです。

for (int j = 1; j < MAX_N; j++) {
    Com[i][j] = (Com[i - 1][j - 1] + Com[i - 1][j]) % P;
}

計算量

  • 前処理：\(O(n^2)\)
  • クエリ：\(O(1)\)

二次元配列上での動的計画法を考えるので、前処理計算に \(O(n^2)\) だけ必要になります。

一度前処理計算してしまえば、あとは \(O(1)\) のクエリで二項係数を利用することができます。

C++ での実装例

#include <bits/stdc++.h>
using namespace std;

const int MAX_N = 5000;
const long long MOD = 1000000007;

long long Com[MAX_N][MAX_N];

/* init_nCk:二項係数のための前処理
    計算量:O(MAX_N*MAX_N)
*/
void init_nCk() {
    memset(Com, 0, sizeof(Com));
    Com[0][0] = 1;
    for (int i = 1; i < MAX_N; ++i) {
        Com[i][0] = 1;
        for (int j = 1; j < MAX_N; j++) {
            Com[i][j] = (Com[i - 1][j - 1] + Com[i - 1][j]) % MOD;
        }
    }
}

/*  nCk :MODでの二項係数を求める(前処理 int_nCk が必要)
    計算量:O(1)
*/
long long nCk(int n, int k) {
    assert(!(n < k));
    assert(!(n < 0 || k < 0));
    return Com[n][k];
}

数学的な説明と原理

そもそも逆元とは

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

例えば、通常の掛け算においては、\(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} $$

逆元の計算方法 (拡張ユークリッドの互除法)

拡張ユークリッドの互除法を利用すると、 i の逆元 を \(O(log P)\) で求めることができます。

$$ax \equiv 1 ~(mod ~P)$$

となる \(x\) は a の逆元そのものなので、これを求める方法を考えましょう。上の式は

$$ax = 1 + Py$$

となる \(y\) が存在することと同値なので、これを変形して

\begin{align}
ax = 1 + Py \\
\Leftrightarrow ax – Py = 1
\end{align}

となる、 \(x\) を求めれば良いです。これは a と P が互いに素の時(特に P が素数でaより大きい時)には、\(gcd(a,P) = 1 \) となるので拡張ユークリッドの互除法で計算が可能です。

long long gcd_ext(long long a, long long b, long long &x, long long &y) {
    if (b == 0) {
        x = 1;
        y = 0;
        return a;
    } else {
        long long d = gcd_ext(b, a % b, y, x);
        y -= a / b * x;
        return d;
    }
}

long long inverse(long long a, long long P) { // aの逆元を求める
    long long x, y;
    gcd_ext(a, P, x, y);
    return (P + x % P) % P;
}

逆元の計算方法(フェルマーの小定理と繰り返し2乗法)

フェルマーの小定理と繰り返し二乗法を利用することでも、 i の逆元 を \(O(log P)\) で求めることもできます。

フェルマーの小定理

P が素数、 x が任意の自然数の時
$$ x^P \equiv x ~(mod ~P) $$
特に、P が素数、 Pとx が互いに素の時に
$$ x^{P-1} \equiv 1 ~(mod ~P) $$

フェルマーの小定理を利用すると、

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

であることが分かるので、繰り返し二乗法によって \(a^{P-2}\) の計算が\(O(logn)\)で求まります。

前処理付きでの逆元の高速な計算方法

さきほど説明した

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

では、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

のようになります。

関連問題

参考になる記事

数学剰余,二項係数,テーマ記事,競プロ,素数,逆元