※上記の広告は60日以上更新のないWIKIに表示されています。更新することで広告が下部へ移動します。

これはCompetitive Programming Advent Calendar Div2012に投稿された記事のようです。
今年ははてなではなくatwiki(こちら)のほうでやります。はてなはTeXが汚かったりソースコードの折りたたみができない等色々アレなので・・。

さて、今回は数え上げ・数学ゲーの常連、Combinationの剰余の計算について書きます。正確にはBinomial Coefficient(2項係数)ですね。ここでは1秒以下で求められることを目安に制約をつけてそれぞれについて解説します。
{}_nC_r=\binom{n}{r}=\frac{n!}{r!(n-r)!}\ (\text{mod}\ m). TeX的に楽な\binom{n}{r}式で書きたいとおもいます。
0o_
※一応Javaコードを添えてはいますが、完全なものとは限らないので、全部欲しい場合は関数名から察するか@uwitenpenにmentionください。


mがsqaurefree(4以上の平方数で割り切れない)のとき

mを素因数分解して素数p_iの積にばらしたあと、各p_iについて\binom{n}{r}を求めることができれば、拡張互除法によりmでの剰余も求められます。拡張互除法がわからなくなったら僕はここに行きます。

コード
+...
public static long[] exGCD(long a, long b)
{
	if(a == 0 || b == 0)return null;
	int as = Long.signum(a);
	int bs = Long.signum(b);
	a = Math.abs(a); b = Math.abs(b);
	long p = 1, q = 0, r = 0, s = 1;
	while(b > 0){
		long c = a / b;
		long d;
		d = a; a = b; b = d % b;
		d = p; p = q; q = d - c * q;
		d = r; r = s; s = d - c * s;
	}
	return new long[]{a, p * as, r * bs};
}
 
public static long crt(long[] divs, long[] mods)
{
	long div = divs[0];
	long mod = mods[0];
	for(int i = 1;i < divs.length;i++){
		long[] apr = exGCD(div, divs[i]);
		if((mods[i] - mod) % apr[0] != 0)return -1;
		long ndiv = div * divs[i] / apr[0];
		long nmod = (apr[1] * (mods[i] - mod) / apr[0] * div + mod) % ndiv;
		if(nmod < 0)nmod += ndiv;
		div = ndiv;
		mod = nmod;
	}
	return mod;
}
 

値(階乗+逆元) (n<10^7, r<10^7, n,r<m, m:素数)

前計算O(\max(n,r)\log \max(n,r)), 本計算O(1).
基本に忠実なコードです。まずmax(n,r)までの階乗x!, 階乗の逆元(x!)^-1を列挙しておいて、
(n!)*((n-r)!)^{-1}*(r!)^{-1}\ \text{mod}\ m
を計算するだけです。この計算はmの剰余類が群にならなければいけない(どの要素にも逆元が一意に決まる)ことが前提なので、汎用的に使うのであれば、mは素数でなければいけません。

コード
+...
public static long C(int n, int r)
{
	if(n < 0 || r < 0 || r > n)return 0;
	if(r > n / 2)r = n - r;
	return FACT[n]*IFACT[n-r]%mod*IFACT[r]%mod;
}
 
static int mod;
static long[] FACT, IFACT;
 
	N = 5000;
	FACT = new long[N+1];
	IFACT = new long[N+1];
	FACT[0] = 1;
	IFACT[0] = 1;
	for(int i = 1;i <= N;i++){
		FACT[i] = FACT[i-1] * i % mod;
		IFACT[i] = invl(FACT[i], mod);
	}
 

値(逆元) (r<10^7, n,r<m, m:素数)

O(r).
これも基本に忠実。n!と(n-r)!の部分を相殺して
\frac{n\cdot (n-1)\cdot\cdots \cdot(n-r+1)}{1\cdot 2\cdot\cdots \cdot r}
を計算します。これならnが大きくても大丈夫。

コード
+...
public static long C(int n, int r, int p)
{
	long ret = 1;
	while(true){
		if(r == 0)break;
		int N = n % p;
		int R = r % p;
		if(N < R)return 0;
 
		for(int i = 0;i < R;i++){
			ret = ret * (N-i) % p;
		}
		long imul = 1;
		for(int i = 0;i < R;i++){
			imul = imul * (i+1) % p;
		}
		ret = ret * Modulo.inv(imul, p) % p;
		n /= p; r /= p;
	}
	return ret;
}
 

列挙(テーブル型) (n*r<10^7, m:任意)

O(nr).
2次元配列で列挙するあれです。次の性質を利用しています。
\binom{n}{r}=\binom{n-1}{r}+\binom{n-1}{r-1}
和は2項までしかとらないので%mをとらずにif文で-mするのもあり。

コード
+...
long[][] C = new long[n+1][n+1];
for(int i = 0;i <= n;i++){
	C[i][0] = 1;
	for(int j = 1;j <= i;j++){
		C[i][j] = (C[i-1][j-1] + C[i-1][j]) % m;
	}
}
return ret;
Javaでは0で初期化されるのでC[2][5]みたいなのにアクセスされても0を返すのですがC++まわりは初期化が必要です。

値(行列冪乗) (r<100, n,m:任意)

rが小さい場合行列の冪乗で書けます。上のテーブル型をDPとした場合の行列DPですね。この行列は下三角Toeplitz行列なのでO(r^2\log n)で計算できます。

\begin{pmatrix} \binom{n}{r} \\ \binom{n}{r-1} \\ \binom{n}{r-2} \\ \vdots \\ \binom{n}{0} \end{pmatrix}=\begin{pmatrix} 1&amp;0&amp;0&amp;\hdots &amp;0&amp;0\\ 1&amp;1&amp;0&amp;\hdots &amp;0&amp;0 \\ 0&amp;1&amp;1&amp;\hdots &amp;0&amp;0 \\ \vdots &amp;\vdots &amp;\ddots &amp;\vdots &amp;\vdots \\ 0&amp;0&amp;0&amp;\hdots &amp;1&amp;1\end{pmatrix}^n \begin{pmatrix} \binom{0}{r} \\ \binom{0}{r-1} \\ \binom{0}{r-2} \\ \vdots \\ \binom{0}{0} \end{pmatrix}, \begin{pmatrix} \binom{0}{r} \\ \binom{0}{r-1} \\ \binom{0}{r-2} \\ \vdots \\ \binom{0}{0} \end{pmatrix}=\begin{pmatrix} 0 \\ 0 \\ 0 \\ \vdots \\ 1 \end{pmatrix}

と大仰なことをかきましたが、早い話がVandermonde's identityを利用して
\binom{2n}{r}=\sum_{k=0}^r \binom{n}{m}\binom{n}{r-m}
\binom{n}{r}=\binom{n-1}{r}+\binom{n-1}{r-1}で繰り返し二乗法しようという話です。

どうみても畳み込みなので、高速剰余変換とか使って高速化できるんですかね・・

コード
+...
public static int[] vander(long n, int r, int mod)
{
	int[] v = new int[r+1];
	v[0] = 1;
	for(int d = 63;d >= 0;d--){
		// Vandermonde
		int[] w = new int[r+1];
		for(int i = 0;i < r+1;i++){
			long sum = 0;
			for(int j = 0;j <= i;j++){
				sum = (sum + (long)v[j] * v[i-j]) % mod;
			}
			w[i] = (int)sum;
		}
		v = w;
 
		if(n<<63-d<0){
			// +1
			for(int i = r;i >= 1;i--){
				v[i] += v[i-1];
			}
		}
	}
	return v;
}
 

値(素因数分解) (n<10^7, r,m:任意)

前計算O(n\log\log n), 本計算O(n/log n)?
n!,(n-r)!,r!を素因数分解して各素数の乗数を求めたあと掛けあわせます。n!を素因数分解したときの素数pの乗数は[n/p]+[n/p^2]+[n/p^3]+\cdots で書けるので・・

コード
+...
public static long C(int n, int r, int mod, int[] primes)
{
	if(n < 0 || r < 0 || r > n)return 0;
	if(r > n / 2)r = n - r;
	int[] a = new int[n];
	for(int i = 0;i < r;i++)a[i] = n-i;
 
	for(int p : primes){
		if(p > r)break;
		for(long q = p;q <= r;q *= p){
			int m = n % q;
			for(int i = m, j = 0;j < r / q;i+=q,j++){
				a[i] /= p;
			}
		}
	}
 
	long mul = 1;
	for(int i = 0;i < r;i++){
		mul = mul * a[i] % mod;
	}
	return mul;
}
 

値(Lucas' Theorem) (m:素数,m<=1000, n,r:任意)

前計算O(m^2), 本計算O(\log n)
上述した乗数の式をうまく料理してやるとLucasの定理になります。n,rをm進数でn_kn_{k-1}\cdots n_1n_0, r_kr_{k-1}\cdots r_1r_0と表す時、
\binom{n}{r}=\prod_i \binom{n_i}{r_i}\ (\text{mod}\ m)と計算できます。

コード
+...
public static long C(long n, long r, int p, long[][] Cs)
{
	long mul = 1;
	while(n>0){
		mul *= Cs[(int)(n%p)][(int)(r%p)];
		mul %= p;
		n /= p; r /= p;
	}
	return mul;
}
 

列挙(乗算型SegmentTree) (n<=10^6, r:すべて, m:任意)

O((n/\log n)*\log(n/\log n)*(\log n)+n\sqrt{n})
Codeforces Round 100 Eを解いた時にTLEがとれなくて苦し紛れに書いたやつです。\binom{n}{*}を列挙します。
乗算を伝播させるSegmentTreeをつくり、各葉はひとつの素数^(乗数) mod mを表します。まっさらな状態からはじめて
  1. nを素因数分解したものを足し, 1を素因数分解したものを引き、ルートの値を取得
  2. n-1を素因数分解したものを足し, 2を素因数分解したものを引き、ルートの値を取得
ということを繰り返していくと列挙できます。結局TLEしたんですがー

+...
public static int[] enumC(int n, int[] primes, int mod)
{
	int m = primes.length;
	int[] ip = new int[n+1];
	for(int i = 0;i < m;i++){
		ip[primes[i]] = i;
	}
	SegmentTreeRMULQSimple st = new SegmentTreeRMULQSimple(m, primes, mod);
	int[] ret = new int[n+1];
	ret[0] = ret[n] = 1;
	int[][][] fs = Num.enumFactor(n, primes);
	for(int i = 1;i <= n/2;i++){
		for(int[] ps : fs[n+1-i]){
			st.add(ip[ps[0]], ps[1]);
		}
		for(int[] ps : fs[i]){
			st.add(ip[ps[0]], -ps[1]);
		}
		ret[i] = ret[n-i] = st.st[1];
	}
	return ret;
}
 
public static class SegmentTreeRMULQSimple {
	public int M, H, N, B;
	public int[] st;
	public int[] ps;
	public int[] es;
	public int mod;
 
	public SegmentTreeRMULQSimple(int n, int[] ps, int mod)
	{
		N = n;
		H = Integer.highestOneBit(Math.max(n,1))<<1;
		M = H<<1;
		B = Integer.numberOfTrailingZeros(H);
		st = new int[M];
		Arrays.fill(st, 1);
		es = new int[H];
 
		this.ps = new int[H];
		Arrays.fill(this.ps, 1);
		for(int i = 0;i < ps.length;i++){
			this.ps[i] = ps[i];
		}
		this.mod = mod;
	}
 
	public void add(int pos, int v)
	{
		es[pos] += v;
		st[H+pos] = (int)Modulo.pow(ps[pos], es[pos], mod);
		for(int i = (H+pos)>>1;i >= 1;i >>>= 1){
			st[i] = (int)((long)st[2*i]*st[2*i+1]%mod);
		}
	}
}
 

値(Lucas' Theoremの拡張) (m=p^q:素数の累乗, m<10^7, n,r:任意)

CodeSprint 3 nCrを解くときに調べたら出てきたpdfのやつです。割と有名らしい?
これまで出てきた\binom{n}{r}の値を求めるものの剰余はだいたい素数になっていました。これを素数の累乗まで拡張できると、小さい方のmではほとんど攻略できるといえるでしょう。

いくつか変数を定義します。例によってn,r,n-r=oをp進数で表します。
n=n_kn_{k-1}\cdots n_1n_0
r=r_kr_{k-1}\cdots r_1r_0
o=o_ko_{k-1}\cdots o_1o_0
また、N_iを、nのp進数表現でのi桁目からi+q-1桁目をつなげた値とします。R_i,O_iも同様に定義します。
N_i=n_i+n_{i+1}p+\cdots +n_{i+q-1}p^{q-1}=[n/p^i]\ \text{mod}\ p^q
そしてe_i
e_i=([n/p^{i+1}]+[n/p^{i+2}]+\cdots)-([o/p^{i+1}]+[o/p^{i+2}]+\cdots)-([r/p^{i+1}]+[r/p^{i+2}]+\cdots)
で定義しておきます。
最後に、(x!)_pを、x!からpの倍数の項を消したもの、つまり(x!)_p=x!/[x/p]!/p^{[x/p]}とします。
このとき、
\frac{1}{p^{e_0}}\binom{n}{r}=(\pm 1)^{e_{q-1}}\left(\frac{(N_0!)_p}{(O_0!)_p(R_0!)_p}\right)\left(\frac{(N_1!)_p}{(O_1!)_p(R_1!)_p}\right)\cdots \left(\frac{(N_k!)_p}{(O_k!)_p(R_k!)_p}\right)\ (\text{mod}\ p^q)
が成り立ちます。\pm 1の部分は、p=2かつq>3の時に限り1, それ以外は-1. \binom{n}{r}の計算において邪魔なpを左辺に移行したものだと考えれば良いです。p^{e_0}を右辺に移項すると\binom{n}{r}が求められます。

前計算はp^qまでの(x!)_pとその逆元を列挙しておけばOK.

コード
+...
public static long C(long n, long r, int p, int q, int[] fact, int[] ifact)
{
	if(n < 0 || r < 0 || r > n)return 0;
	int P = 1;
	for(int i = 0;i < q;i++)P*=p;
 
	long z = n - r;
	int e0 = 0;
	for(long u = n/p;u > 0;u /= p)e0 += u;
	for(long u = r/p;u > 0;u /= p)e0 -= u;
	for(long u = z/p;u > 0;u /= p)e0 -= u;
 
	int em = 0;
	for(long u = n/P;u > 0;u /= p)em += u;
	for(long u = r/P;u > 0;u /= p)em -= u;
	for(long u = z/P;u > 0;u /= p)em -= u;
 
	long ret = 1;
	while(n > 0){
		ret = ret * fact[(int)(n%P)] % P * ifact[(int)(r%P)] % P * ifact[(int)(z%P)] % P;
		n /= p; r /= p; z /= p;
	}
	for(int i = 0;i < e0;i++)ret = ret * p % P;
	if(!(p == 2 && q >= 3) && (em&1)==1)ret = (P-ret) % P;
	return ret;
}
 
public static int[][] makeFF(int p, int q)
{
	int P = 1;
	for(int i = 0;i < q;i++)P*=p;
	int[] fact = new int[P+1];
	int[] ifact = new int[P+1];
	fact[0] = 1;
	for(int i = 1;i <= P;i++){
		if(i % p == 0){
		fact[i] = fact[i-1];
		}else{
			fact[i] = fact[i-1] * i % P;
		}
	}
	for(int i = 0;i <= P;i++){
		long ret = 1;
		long mul = fact[i];
		for(long n = P/p*(p-1)-1;n > 0;n >>>= 1){
			if((n&1)==1){
				ret = (ret * mul) % P;
			}
			mul = (mul * mul) % P;
		}
		ifact[i] = (int)ret;
	}
	return new int[][]{fact, ifact};
}
 

値(Lucas' Theoremの拡張2) (pq<10^7)

上の方式でmが小さい時はいけるようになりましたが、qが大きいと(x!)_pが列挙できずに破綻してしまいます。次の2個の公式を使って(x!)_pの列挙範囲をx\le pqでよくします。

\left(\frac{( (up+v)!)_p}{(up!)_p(v!)_p}\right)=\prod_{j=1}^{q-1}\left(\frac{( (jp+v)!)_p}{(jp!)_p(v!)_p}\right)^{\alpha_j}\ (\text{mod}\ p^q)).
  • vは0以上p-1以下
  • \alpha_j=\frac{u}{j}\prod_{1\le i\le q-1, i\neq j}\frac{u-i}{j-i}.
\alpha_jは負までいった二項係数の計算に見えますね。整数にはなるみたいなので、桁あふれをどう抑えるかはまだ課題です。mod p^qだったので\alpha_j\varphi(p^q)=(p-1)p^{q-1}で剰余とれそうですね。

これと、もうひとつ、
(up!)_p=\pm \prod_{j=1}^r (jp!)_p^{\beta_j}\ (\text{mod}\ p^{2r+1})
が、p^r=2または2r+1=p\ or\ p^2以外のときに成り立ちます。成り立たない場合は大きめにrをとればいいだけの話です。
  • \pm はp=2のときのみマイナスになる。
  • \beta_j=\frac{u}{j}\prod_{1\le i\le r, i\neq j}\frac{u^2-i^2}{j^2-i^2}}.
この\beta_jも桁あふれこわいですね・・これも\varphi(p^{2r+1})=(p-1)p^{2r}で剰余とれそう。

q=2r+1となるように調整して2個の公式の左辺同士をかけると\frac{( (up+v)!)_p}{(v!)_p}がでてきて、vはp-1以下で、jp+vもjpもpq以下なので、pq以下の(x!)_p( (up+v)!)_pが求められることになります!良いですね!\alpha_j, \beta_jがうまく計算できれば実用的そうです。

コード(実際にはP,betanum,alphanumが馬鹿でかくなるのであまり使えませんね・・)
+...
public static long ff(long n, int p, int q, int[] fact, int[] ifact)
{
	long u = n/p;
	int v = (int)(n%p);
	int r = (q-1)/2;
	long P = 1;
	for(int j = 1;j <= q;j++)P *= p;
 
	long ret = fact[v];
	for(int j = 1;j <= r;j++){
		long betanum = u;
		long betaden = j;
		for(int i = 1;i <= r;i++){
			if(i != j){
				betanum *= u*u-i*i;
				betaden *= j*j-i*i;
			}
		}
		long d = betanum / betaden;
		long po = Modulo.pow(fact[j*p], Math.abs(d), P);
		if(d < 0)po = Modulo.invl(po, P);
		ret = ret * po%P;
 
	}
	if(p == 2)ret = -ret;
	for(int j = 1;j <= q-1;j++){
		long cell = (long)fact[j*p+v]%P
				*ifact[j*p]%P
				*ifact[v]%P;
		long alphanum = u;
		long alphaden = j;
		for(int i = 1;i <= q-1;i++){
			if(i != j){
				alphanum *= u-i;
				alphaden *= j-i;
			}
		}
		long d = alphanum / alphaden;
		long po = Modulo.pow(cell, Math.abs(d), P);
		if(d < 0)po = Modulo.invl(po, P);
		ret = ret * po%P;
	}
	return (ret+P)%P;
}
 

おまけ

\binom{n}{2}=\frac{n(n-1)}{2}\ (\text{mod}\ m)でnがでかくてmがやばい場合、n,n-1のいずれかは偶数なので、先に偶数の方を2で割ってしまえば幸せになります。
この手の式(\frac{n(n+1)}{2}, \frac{n(n+1)(2n+1)}{6}, \frac{n(n-1)(n-2)}{3})には同じ考え方が使えます。ProjectEulerの制約サイズが最近大きくてここによくはまるので自戒も込めて。
mがintをこえてさらにやばい場合はlongの範囲で収まる例の乗算をすればよいですね。
+...
public static long c2(long n, int mod)
{
	if(n % 2 == 0){
		return (n/2%mod)*((n-1)%mod)%mod;
	}else{
		return (n%mod)*((n-1)/2%mod)%mod;
	}
} 

あとがき

さて、いろいろ書きなぐってきましたが、これだけあれば競プロで出てくる二項係数には困らないんじゃないでしょうか。
ただこれだけあっても素早く求まらない二項係数はいくらでもあります。そのときは恥ずかしがらずに教えて下さいね。