nCr mod mの求め方

これは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以上の平方数で割り切れない)のとき
値(階乗+逆元) (n<10^7, r<10^7, n,r<m, m:素数)
値(逆元) (r<10^7, n,r<m, m:素数)
列挙(テーブル型) (n*r<10^7, m:任意)
値(行列冪乗) (r<100, n,m:任意)
値(素因数分解) (n<10^7, r,m:任意)
値(Lucas' Theorem) (m:素数,m<=1000, n,r:任意)
列挙(乗算型SegmentTree) (n<=10^6, r:すべて, m:任意)
値(Lucas' Theoremの拡張) (m=p^q:素数の累乗, m<10^7, n,r:任意)
値(Lucas' Theoremの拡張2) (pq<10^7)
おまけ
あとがき

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

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

コード

...

値(階乗+逆元) (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は素数でなければいけません。

コード

...

値(逆元) (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が大きくても大丈夫。

コード

...

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

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

コード

...

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&0&0&\hdots &0&0\\ 1&1&0&\hdots &0&0 \\ 0&1&1&\hdots &0&0 \\ \vdots &\vdots &\ddots &\vdots &\vdots \\ 0&0&0&\hdots &1&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}$ で繰り返し二乗法しようという話です。

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

コード

...

値(素因数分解) (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$ で書けるので・・

コード

...

値(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)$ と計算できます。

コード

...

列挙(乗算型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を表します。まっさらな状態からはじめて

nを素因数分解したものを足し, 1を素因数分解したものを引き、ルートの値を取得
n-1を素因数分解したものを足し, 2を素因数分解したものを引き、ルートの値を取得

ということを繰り返していくと列挙できます。結局TLEしたんですがー

...

値(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.

コード

...

値(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が馬鹿でかくなるのであまり使えませんね・・)

...

おまけ

$\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の範囲で収まる例の乗算をすればよいですね。

...

あとがき

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

「nCr mod mの求め方」をウィキ内検索

最終更新：2012年12月29日 05:00

uwicoder