ABC 034 C

昨日この問題を解きました。

abc034.contest.atcoder.jp

wとhの値が与えられます。(w+h)Cw % (10^9+7)の値を求めよ、という問題です。
wとhの値が小さければ動的計画法が使えるのですが、今回はwとhのどちらもmax10^5ということで使えません。
今回は(w+h)Cw = (w+h)!/(w!h!)という公式を用います。
w!h!という値はlong long int の範囲すら超える可能性があります。そこでフェルマーの小定理を利用します。
フェルマーの小定理はこんなかんじです。
aとpが互いに素のとき、
a^(p-1) % p = 1 % p
さらに、
a^(p-2) % p = a^(-1) % p

で、これを用いて、
(w+h)Cw % p = (w+h)!/(w!h!) % p
= (((w+h)! % p) * ( (w!h!) % p) ) % p
= (((w+h)! % p) * (((w!^(-1) % p) * (h!^(-1) % p)) % p)) % p
= (((w+h)! % p) * (((w!^(p - 2) % p) * (h!^(p - 2) % p)) % p)) % p
という感じに変形して計算します。

答えはこんな感じ。

    #include <cstdio>
    #include <iostream>
    #include <cmath>
    #include <cstring> 
    #include <sstream>
    #include <iostream>
    #include <algorithm>
    #include <cstdlib>
    #include <map>
    #include <queue>
    #include <utility>
    #include <vector>
    #include <set>
    #include <memory.h>
    #include <iomanip>
    #include <bitset>
    #include <list>
    #include <stack>
     
    using namespace std;
     
    #define pai 3.14159265359
    #define mod 1000000007
     
    int main()
    {
        int w, h;
        cin >> w >> h;
        w--;
        h--;
        // (w + h)Cw = (w + h)! / (w! * h!)を求める
        // ((w + h)! / (w! * h!)) % mod = (((w + h)! % mod) * (((w! * h!) ^ (-1)) % mod)) % mod
                                     // = (((w + h)! % mod) * ((((w!) ^ (-1)) % mod) * (((h!) ^ (-1)) % mod) % mod)) % mod
        // フェルマーの小定理より　modは素数だから ((w!) ^ (mod - 1)) % mod = 1 % mod よって ((w!) ^ (-1)) % mod = ((w!) ^ (mod - 2)) % mod
        long long int factw = 1;
        for(int i = 2; i <= w; i++) factw = (factw * i) % mod;
        long long int inverse_factw_mod = 1;
        long long int pow_factw_mod[50];
        // pow_factw_mod[i] = factw ^ (2^i) % mod
        pow_factw_mod[0] = factw % mod;
        for(int i = 1; i < 50; i++){
            pow_factw_mod[i] = (pow_factw_mod[i - 1] * pow_factw_mod[i - 1]) % mod;
        }
        long long int tmp = mod - 2;
        for(int i = 49; i >= 0; i--){
            if(tmp >= (long long int)pow(2, i)){
                inverse_factw_mod = (inverse_factw_mod * pow_factw_mod[i]) % mod;
                tmp -= (long long int)pow(2, i);
            }
        }
        long long int facth = 1;
        for(int i = 2; i <= h; i++) facth = (facth * i) % mod;
        long long int inverse_facth_mod = 1;
        long long int pow_facth_mod[50];
        // pow_factw_mod[i] = factw ^ (2^i) % mod
        pow_facth_mod[0] = facth % mod;
        for(int i = 1; i < 50; i++){
            pow_facth_mod[i] = (pow_facth_mod[i - 1] * pow_facth_mod[i - 1]) % mod;
        }
        tmp = mod - 2;
        for(int i = 49; i >= 0; i--){
            if(tmp >= (long long int)pow(2, i)){
                inverse_facth_mod = (inverse_facth_mod * pow_facth_mod[i]) % mod;
                tmp -= (long long int)pow(2, i);
            }
        }
        long long int fact_h_plus_w = 1;
        for(int i = 2; i <= h + w; i++) fact_h_plus_w = (fact_h_plus_w * i) % mod;
        // cout << fact_h_plus_w << " " << inverse_facth_mod << " " << inverse_factw_mod << endl;
        cout << (fact_h_plus_w * ((inverse_facth_mod * inverse_factw_mod) % mod)) % mod << endl;
        return 0;
    }

フェルマーの小定理、なんだか不思議です。a^(-1)のmodってとれなさそうなのに。