[競プロ] nCrの(mod 素数)を計算したい

競プロをしていたらよく、nCrを1000000007で割った余りを求めたいことがあると思う。
これの求め方を自分の復習がてら説明する。
以下、一般に M で割った余り、すなわちmod M について書く。

前提知識

この記事でも軽く説明するが、高校数学で習う nCrの意味とmodについての知識があると良い。

まず、nCrは「n個の中からr個選ぶ時の選び方の総数」である。

数式にすると,

{}_nC_r=\frac{n!}{r!(n−r)!}

である。
これをコードで表すために軽く式変形すると

{}_nC_r=\frac{n!}{r!(n−r)!}=\frac{n(n-1)(n-2) …(n-r-1)}{1・2・3・…(r-1)r}

のようになる。

ここまでをC++で書くと、分母をy,分子をxとして

for (int i = 0; i < r; i++) {
	x = x * (n - i)  ;
	y = y * (i + 1);
}

{}_nC_r=\frac{x}{y}

の用になる。
これをMODで取ってやりたいのだが、MODの性質として足し算、引き算、掛け算は同じMOD同士のもので演算して良いのだが、割り算はとある条件の時以外してはいけないという物がある。

つまり,

for (int i = 0; i < r; i++) {
	x = x * (n - i)%M;
	y = y * (i + 1)%M;
}

のように分母分子をmod Mで取った後に x/yを計算してはだめなのだ。

ここで

一般に、yで割るというのは 1/yをかける　と言い換えられるのはいいだろうか。
ここでy * (1/y)=1である。
これをmodの世界でも適用してあげると yで割るという演算の代わりに、

y * y'=1 (mod M)

となるような y'を見つけてあげてかければよいのである。

ここで、†フェルマーの小定理†という小難しい名前の定理がある。
定理の主張自体は全く難しいものではなく、

pを素数、aを pの倍数でない任意の整数としたとき、

a^{p−1} ≡ 1 \  (mod\ p)

が成立するというものだ。

ここでpをMに、aをyに置き換えて、少しだけ式変形をすると

y*y^{M−2} ≡ 1 \  (mod\ M)

となり、上記のy'が y^(M-2)に相当することがわかる。

ゆえに求めたいものは

for (int i = 0; i < r; i++) {
	x = x * (n - i)%M;
	y = y * (i + 1)%M;
}

// ans=( x * pow(y,M-2) )% M

となる。なお、 pow(a,b)はa^bを表す
以上が本質部分であり、今回の主題であったが、実装に置いては細かいことが残っている。
それはC++の標準のpow()だと遅い上に、計算過程でオーバーフローしてしまうことである。
ゆえに逐次modを取りながらべき乗の計算をしなければならない。
ここでは詳しく語らないが、pow()高速化をするのに繰り返し二乗法というアルゴリズムがある。
n乗の計算を、標準のアルゴリズムでは o(n) かかるのに対し、 o(log n)　で計算するものである。

以上のものを使って書いたコードが以下だ。

//繰り返し二乗法
long long pow_mod(long long x, long long n) {
	long long ret = 1;
	while (n > 0) {
		if (n & 1) ret = ret * x % MOD;  // n の最下位bitが1ならば(nが奇数ならば) x^(2^i) をかける
		x = x * x % MOD;
		n >>= 1;  // n を1bit 左にずらす
	}
	return ret;
}


long long nCr_mod(long long n, long long r) {
	long long x = 1, y = 1;

	for (int i = 0; i < r; i++) {
		x = x * (n - i) % MOD ;
		y = y * (i + 1) % MOD;
	}

	return x * pow_mod(y, MOD - 2) % MOD;
}