フェルマーの小定理

LayCurseさんの日記を見て面白そうなのでチャレンジした問題。

http://www.codeforces.com/contest/57/problem/C

①重複組み合わせ（Homogeneous Combination）と

②包除原理（Inclusion-exclusion principle）

の練習にと解いてみたが、見事にhogeりました。。

LayCurseさんのソースと睨めっこすること約2日、どうやら素数pと自然数aに対して

a^p = a (mod p)

が成り立つのではないか？？と気付きました。

これは「フェルマーの小定理」と呼ばれているようです。

別の書き方をすれば、

a^(p-1) = 1 (mod p)や

a^(p-2) = 1/a (mod p)

です。

上の問題は、この定理を用いて重複組み合わせを計算します。

普通、組み合わせはパスカルの三角形を利用してDPで出せますが、この問題の場合は抽出対象の個数が大きいためメモリーオーバーとなり、NGです。

フェルマーの小定理を利用すると、法の世界の割り算ができます。

普通に

nCk = n! / (n-k)!k!

を計算しましょう。

pが大きいので、累乗の計算は繰り返し二乗法で行いましょう。

以下ソースです。

view plainprint?

1. const int MOD = 1000000007;  
2. long long int f[200001];  
3. long long int rf[200001];  
4.   
5. long long int pw(long long int x, long long int y) {  
6.     long long ret = 1LL;  
7.   
8.     REP(mask, 50) {  
9.         if (y >> mask & 1)  
10.             ret = ret * x % MOD;  
11.         x = x * x % MOD;  
12.     }  
13.     return ret % MOD;  
14. }  
15.   
16. long long int homoComb(int x, int y) {  
17.     int n = x + y - 1;  
18.     int k = y;  
19.   
20.     rf[0] = f[0] = 1;  
21.     FOR (i, 1, n + 1) {  
22.         f[i] = f[i-1] * i % MOD;  
23.         rf[i] = pw(f[i], MOD-2);  
24.     }  
25.   
26.     long long int ret = f[n];  
27.     ret = ret * rf[k] % MOD;  
28.     ret = ret * rf[n-k] % MOD;  
29.   
30.     return ret;  
31. }  
32.   
33. int main() {  
34.     int n;  
35.   
36.     cin >> n;  
37.     memset(f, 0x00, sizeof(f));  
38.     memset(rf, 0x00, sizeof(rf));  
39.   
40.     long long int ret = homoComb(n, n) % MOD;  
41.     ret = 2 * ret % MOD;  
42.     ret -= n;  
43.     if (ret < 0) ret += MOD;  
44.   
45.     cout << ret << endl;  
46.   
47.     return 0;  
48. }  

const int MOD = 1000000007;
long long int f[200001];
long long int rf[200001];

long long int pw(long long int x, long long int y) {
    long long ret = 1LL;

    REP(mask, 50) {
        if (y >> mask & 1)
            ret = ret * x % MOD;
        x = x * x % MOD;
    }
    return ret % MOD;
}

long long int homoComb(int x, int y) {
    int n = x + y - 1;
    int k = y;

    rf[0] = f[0] = 1;
    FOR (i, 1, n + 1) {
        f[i] = f[i-1] * i % MOD;
        rf[i] = pw(f[i], MOD-2);
    }

    long long int ret = f[n];
    ret = ret * rf[k] % MOD;
    ret = ret * rf[n-k] % MOD;

    return ret;
}

int main() {
    int n;

    cin >> n;
    memset(f, 0x00, sizeof(f));
    memset(rf, 0x00, sizeof(rf));

    long long int ret = homoComb(n, n) % MOD;
    ret = 2 * ret % MOD;
    ret -= n;
    if (ret < 0) ret += MOD;

    cout << ret << endl;

    return 0;
}