[bzoj1951][SDOI2010]古代猪文

题目描述:

计算G^(Sigma(C(n,i)) % P， 其中 i | n

题目解答：

为叙述方便起见，我们令sigma()为W
显然W是一个非常大的数，如果暴力直接计算不是明智的选择。
我们这里使用费马小定理：
对于一个素数p：a^(p-1) = 1| (a, p) = 1
对于G^W, 确定一个常数k，有W = k * (p-1) + M
所以G^W = (G^(k(p-1)) G^M
显然G^W就等于G^M，所以问题就变成了求出组合数C(n,i)%(P-1)的值。
计算组合数的时候，有两种方法：
1) 暴力。因为组合数公式中有阶乘，所以用O(n)的时间暴力求
2) 用空间换时间，预处理出fact(n).
显然，如果直接算，我们需要预处理出1000000000这么多的fact，空间不允许。
所以我们介绍lucas定理：
C(n, m, p) = C(n / p, m / p, p) * C(n % p, m % p, p)
这样我们就有了高效的方法求组合数。
但是lucas定理的应用有一个局限：p必须是素数。
我们观察到：p-1并不是一个素数。
这里介绍一个分解因数的软件:在linux下直接调用factor就可以分解因数，非常的方便，而且在NOI Linux中会预装此软件。
分解因数，我们就可以应用中国剩余定理来合并线性方程。
这里给出中国剩余定理：
x = a1(mod b1)
x = a2(mod b2)
.....
.....
x = an(mod bn)
那么，对于pi(bi)的剩余系，令P = pi(bi)
x = sigma(a1 * (P/b1) * inv(P/b1, b1))
其中，inv()为逆元，可以使用费马小定理求解（条件：P为素数）
这样我们就解决了这个问题。
这个问题有一个陷阱：
当G == P时，(G, P) = 1, 费马小定理或者说欧拉定理在这里不适用，所以我们应该进行特判。
另外，在数论程序编写时要特别注意模与乘法的问题。
比如：x = a * b
最好写成：
(x = (ll) a * b % p) %= p
比较保险。
下面给出代码。

代码

#include <bits/stdc++.h>
using namespace std;
#define ll long long
const int P0=999911659;
const int P1=999911658; 
int G, N, M[4]; 
int t[4] = {35617, 4679, 3, 2};
int fac[66666];
int pow(int a, int b, int p) { //calc a^b % p
    int ans = 1;
    for(int i = b; i; i >>= 1, a = (ll)a*a%p) {
        if(i & 1) ans = (ll)ans * a % p;
    }
    return ans;
}
/*int pow(int x, int y, int P) {
    int ret = 1;
    while(y) {
        if(y&1) ret = (ll)ret * x % P;
        x = (ll)x*x%P;
        y >>= 1;
    }
    return ret;
}*/
int inv(int x, int P){return pow(x, P-2, P);}
int C(int n, int m, int P) { //calc (n, m) % t[x]
    if (n < m) return 0;
    return (ll)fac[n] * inv(fac[n-m], P) % P * inv(fac[m], P)%P;
}
int lucas(int n, int m, int P) {
    if(!n && !m) return 1;
    return (ll)lucas(n/P, m/P, P) * C(n%P, m%P, P) % P;
}
int main() {
    scanf("%d %d", &N, &G);
    if(G == P0) {
        puts("0");
        return 0;
    }
    for(int i = fac[0]=1; i <= t[0]; i++) fac[i] = (ll)fac[i-1]*i % P1;
    int tot = 0;
    for(int k = 0; k < 4; k++) {
        int now = 0;
        for(int i = 1; i * i <= N; i++) { if(N % i == 0) {
            (now += lucas(N, i, t[k])) %= t[k];
            if(N / i != i) {
                (now += lucas(N, N/i, t[k])) %= t[k];
            }
        }}
        (tot += (ll)now * (P1/t[k]) % P1 * inv(P1/t[k], t[k]) % P1) %= P1;
    }
    printf("%d\n", pow(G, tot, P0));
    return 0;
}