コード

struct combLarge{
    typedef long int_type;
    vector<int_type> fact;
    vector<int_type> revFact;
    //in:nCrのn
    void init(int_type N){
        N++;
        fact.assign(N,0);
        revFact.assign(N,0);
        fact[0] = 1;
        for( int i = 1 ; i < N ; i++ ){
            fact[i] = fact[i-1] * i;
            fact[i] %= mod;
        }
        revFact[N-1] = powMod(fact[N-1] , mod - 2);
        for(int i = N-2 ; i >= 0 ; i-- ){
            revFact[i] = revFact[i+1] * (i+1);
            revFact[i] %= mod;
        }
    }
    //in:nCrのn，r out:nCr % mod
    int_type getC(int_type a , int_type b){
        return (((fact[a] * revFact[b]) % mod) * revFact[a-b] ) % mod;
    }
};

割り算の余剰

くそデカイ数とくそデカイ数の割り算の結果を大きい素数で割った余りを求まるのと，それぞれの大きい素数で割った余りの割り算をするのは意味が異なる（小数点以下云々のせいか？）ので，　逆元を求めて，掛け算に落とし込む必要があるらしい．

要するに次式にする．
$$\frac{a}{b}=a・b^{-1}$$
$b^{-1}$も無理数になることがあるので，$b^{-1}$をかけることと等価になるものを求めなければならない(?)．
フェルマーの小定理というのがあり（バカなので知らなかった)，aとpが互いに素のとき次式で表せる．

$$a^{p-1}\equiv1 \ mod \ p$$

これも知らなかったが，$x$と$y$を$n$で割った余剰が等しいとき次式でかけるらしい．
$$x \equiv y \ mod \ n $$

先ほどのフェルマーの小定理の式の両辺に$a^{-1}$をかけると，次式になる．
$$a^{p-2}\equiv a^{-1} \ mod \ p$$

求めたいaをbで割った答えのnで割った余りは，aにbをp-2乗したものを掛けたものをnで割った余りを求めればよい．

組み合わせ

nCrの定義は次式
$$nCr=\frac{nPr}{r!}=\frac{n!}{r!(n-r)!}=n!・(r!)^{-1}・(n-r)!^{-1}$$

求めたいのは余剰なので
$$nCr\equiv n!・(r!)^{-1}・(n-r)!^{-1}\equiv n!・(r!)^{p-2}・(n-r)!^{p-2} \ mod \ p$$

コードの説明

階乗とその逆元を$O(n)$で配列に格納し，nCrを$O(1)$で求められるようにしている．
考えた人は天才だと思った．

階乗とその逆元の格納

1.配列(vector)初期化 ($O(1)$?)

N++;
fact.assign(N,0);
revFact.assign(N,0);

入力nに対してnCrを求めたいので，ここでn++している．

2.$N!$までの値の余剰を求め配列Factに格納 ($O(N)$)

fact[0] = 1;
for( int i = 1 ; i < N ; i++ ){
    fact[i] = fact[i-1] * i;
    fact[i] %= mod;
}

3.$N!$の逆元を計算 ($O(logN)$)

revFact[N-1] = powMod(fact[N-1] , mod - 2);

long powMod( long a , long p ){
    if( p == 0 ) return 1;
    else if( p == 1 ) return a % mod;
    else if( p % 2 == 0 ){
        long tmp = powMod( a , p/2 );
        return (tmp * tmp) % mod;
    }
    else{
        long tmp = powMod( a , p/2 );
        return ((tmp%mod * tmp%mod) * a ) % mod;
    }
}

累乗の余剰を求めるのは，ビットシフトでの実装の方が早いが見た目がわかりやすいので再帰で実装している．

4.階乗の逆元の余剰を求めて配列revFactに格納 ($O(N)$)

for(int i = N-2 ; i >= 0 ; i-- ){
    revFact[i] = revFact[i+1] * (i+1);
    revFact[i] %= mod;
}

こういう解釈でいいのだろうか......

(n+1)!^{-1}=\frac{1}{1}・\frac{1}{2}・...・\frac{1}{n}・\frac{1}{n+1}\\
(n+1)!^{-1}・(n+1)=\frac{1}{1}・\frac{1}{2}・...・\frac{1}{n}・\frac{1}{n+1}・(n+1)\\
(n!)^{-1}=\frac{1}{1}・\frac{1}{2}・...・\frac{1}{n}

よって$((n+1)!)^{-1}・(n+1)=(n!)^{-1}$なので次式が成り立つ(?)

((n+1)!)^{p-2}・(n+1)\equiv((n+1)!)^{-1}・(n+1)\equiv(n!)^{-1}\equiv(n!)^{p-2} \ mod \ p

オーダーは$O(1+N+logN+N)=O(N)$となる

nCrを返す

あらかじめ階乗とその逆元を求めているので，O(1)でnCrを出力できる．

return (((fact[a] * revFact[b]) % mod) * revFact[a-b] ) % mod;

感想

フェルマーの小定理とかすっかり忘れていたのでいい復習になった気がする．
考えた人は天才だと思った．

ハムスターでも分かった組み合わせの余剰

コード

割り算の余剰

組み合わせ

コードの説明

階乗とその逆元の格納

nCrを返す

感想