2021C37 mod4 競プロ(nCk modP)

今年の東大数学で2021C37を4で割った余りを求める問題が出題されましたので、これに関連してnCk modPについてまとめておきたいとお思います(既存のものの紹介なのでNeuesは少ないです)。

まず競プロでお馴染みのp=10^9+7の場合です。
これはフェルマーの小定理より、a/b modpを求める場合Pが素数でb|pならbの逆元がb^(P-2)ということから、これをかけてあげれば簡単に求めることができます。

問題となるのはPが素数でない場合です。一つの方法として、パスカルの三角形の性質を利用して、n、kが小さいところから段階的に求めていくという方法があります。これはPの性質によらず使える方法で、2021C37くらいなら競プロでも余裕で通せます(大体O(n^2)くらい)。
ただ、これはnが大きくなると時間がかかってしまうのでもう少し速いアルゴが欲しいものです。

私が考えた一つの方法は以下のようになります。
実際に2021C37を出すには、(2021*2021*,,,,,*1985)/(37!)を計算すれば良いです。そのため、1~37、1985~2021を高速で素因数分解したら(vplとかに収納)、それを用いて分母の素因数分解、分子の素因数分解を出します。totalの素因数分解は、分子の素因数の指数ー分母の素因数の指数　をそれぞれのprimeについて計算すれば良いのでそれほどかかりません。そしたら、あとは順に掛けつつmod4で適宜あまりに変換する、という計算をすればわりかし早く出せます。計算量は0(NlogN)くらいかな。これは悪くはないですが、どうやらもっと速く出す方法があるみたいです。

以下の方法は日本語でサーチしても余りヒットするものではないので、あまり知られていないものと思われます。海外向けの競プロの大会でこれを使う問題があったようなので、それで少し話題になったようです。問題のurlは以下です。
http://www.lightoj.com/volume_showproblem.php?problem=1318
Pが素数pを使ってP=p^qとかける場合を考えます。
実はこの時、Andrew Granville's Theoremと呼ばれるものが存在し、これを用いることでO(log^2n+Plog^2P)くらいで高速計算ができるようです。以下、これを用いてnCr mod p^qを出す方法を説明します。参考文献の論文:"Binomical coefficients modulo prime powers" by Andrew Graville

まずn=k+rとすると、k=k0+k1*p+,,,,,+kd*p^dとかけます。そして、Kj(j>=0)を、floor(k/(p^j)) modp^qと定義し、Kdまで計算しておきます。同様にNj、Rjを定義し、それぞれN0~Nd、R0~Rdまで計算しておきます。
k=k0+k1*p+,,,,,+kd*p^dより、k=k0+k1*p+,,,,,+kd*p^d+0*p^(d+1)+o*p^(d+2),,,となります(係数をkiとする)。同様の展開をnに対しても行います(係数をkjとする)。kiが0となり、その後も永遠に0が続くようなiと、niが0となり、その後も永遠に0が続くようなiの最大値をとり(これをzとする)、そこまでそれぞれki、niを計算し2D配列にでも収容しておきます。
この時、ejをki>niとなるようなi(i>=j)の個数と定義します。累積和とかを使えばe0~ezは効率的にもとまります。
これらを使えば、nCk mod p^qは以下のように書けるというのがこの定理の主張です。
p^(e0)*((+-1)^(e(q-1)))* (N0!/(K0!*R0!)) * (N0!/(K0!*R0!)) *,,,,,, * (Nd!/(Kd!*Rd!))
ここで(+-1)の部分ですが、((p==2)&&(q>=3))の時にここが-1、elseで1となるようです。
この定理の証明を探しているのですがなかなか見つかりません。ご存知でしたら教えてください。

これはかなり高速なので2021C37を比較的少ないステップ数で出せます。

nCk mod p^q(p=prime)がわかるというのはかなり意義があります。というのは、以下のようにしてPが一般の場合に拡張できるからです。
P=(p0^a0)*(p1^a1)*,,,,のように書けるとします。その時、それぞれのprime p0, p1,,, に対して上の公式を使ってmodがpo^a0の時、p1^a1の時、、、
についてnCr を出すことができます。
簡略化のためにP=(p0^a0)*(p1^a1)でp0^a0=x、p1^a1=yとしますと、
nCk≡x1(modx)
nCk≡y1(mody)
となるx1、y1を出すことができます。
この時、xとyが互いに素の場合、
z≡x1(modx), z≡y1(mody)となるzが0からx*yの間に存在することが中国剰余定理よりわかるので、それを求めてあげれば解決します。
このプロセスを一般化してあげることでnCk mod Pを高速で求めることができることがわかりました。