【bzoj1951】【古代猪文】Lucas定理+欧拉定理+孙子定理

（上不了p站我要死了，当然是游戏原画啊）

Description
（题面倒是很有趣，就是太长了）
题意：
一个朝代流传的猪文文字恰好为N的k分之一，其中k是N的一个正约数（可以是1和N）。不过具体是哪k分之一，以及k是多少，由于历史过于久远，已经无从考证了。考虑到所有可能的k。显然当k等于某个定值时，该朝的猪文文字个数为N / k。然而从N个文字中保留下N / k个的情况也是相当多的。如果所有可能的k的所有情况数加起来为P的话，那么他研究古代文字的代价将会是G的P次方。 现在他想知道研究古代文字的代价除以999911659的余数是多少。
Input
有且仅有一行：两个数N、G，用一个空格分开。
Output
有且仅有一行：一个数，表示答案除以999911659的余数。
Sample Input
4 2
Sample Output
2048
HINT
10%的数据中，1 <= N <= 50；
20%的数据中，1 <= N <= 1000；
40%的数据中，1 <= N <= 100000；
100%的数据中，1 <= G <= 1000000000，1 <= N <= 1000000000。

若留下来的文字个数为x，则可能的情况有C(n,x)个。用o(sqrt(n))的复杂度来算出所有的可能情况。这样就可以统计出指数了。

但是这道题的精髓是取模上。
设指数为x，则由欧拉定理得G^x≡G^(x%φ(M)) (mod M)，gcd(G,M)==1
M=999911659，但是却惊讶的发现φ(M)=999911658，并不是一个质数。这下就GG了，因为题目的数据范围很明显要用Lucas定理，但是Lucas需要模数是质数。

这里就需要用上Lucas的一个拓展了。虽然模数不是质数不能直接用Lucas，但是如果模数可以分解成若干个不同的质数（每个质数最多出现一次）的乘积，就可以用孙子定理（中国剩余定理）来解决。
具体操作：
999911658=2*3*4679*35617
则可得出指数x需满足：x≡a1(mod 2), x≡a2(mod 3), x≡a3(mod 4679), x≡a4(mod 35617)
所以就用孙子定理（发现基本的孙子定理比拓展的孙子定理好些多了）。

还有一点需要注意：
欧拉定理的使用条件是gcd(G,M)==1，如果G==M则要输出0。
（就在这里wa了）

AC代码

#include<cstdio>
#include<cstring>
#include<algorithm>
using namespace std;
#define ll long long 
#ifdef WIN32
#define RIN "%I64d"
#else
#define RIN "%lld"
#endif

const ll mod=999911659;

ll g[5]={0,2,3,4679,35617},sum[5];
ll jiec[36000],niy[36000];
ll n,G;

void exgcd(ll a,ll b,ll &x,ll &y){
    if(b==0){
        x=1,y=0;return;
    }
    ll x0,y0;
    exgcd(b,a%b,x0,y0);
    x=y0;
    y=x0-(a/b)*y0;
}
ll inverse(ll a,ll b){
    ll x,y;
    exgcd(a,b,x,y);
    return (x%b+b)%b;
}
void init(int x){
    memset(jiec,0,sizeof(jiec));
    memset(niy,0,sizeof(niy));
    jiec[0]=niy[0]=1;
    for(int i=1;i<x;i++) jiec[i]=jiec[i-1]*i%x;
    niy[x-1]=inverse(jiec[x-1],x);
    for(int i=x-2;i>=1;i--) niy[i]=niy[i+1]*(i+1)%x;
}
ll comb(ll a,ll b,ll x){
    return jiec[a]*niy[b]%x*niy[a-b]%x;
}
ll lucas(ll a,ll b,ll x){
    if(a<b) return 0;
    if(b==0) return 1;
    if(a<x&&b<x) return comb(a,b,x);
    return lucas(a/x,b/x,x)*lucas(a%x,b%x,x)%x;
}
ll power(ll a,ll b){
    ll rt=1;
    for(;b;b>>=1,a=(a*a)%mod)
        if(b&1) rt=(rt*a)%mod;
    return rt;
}
int main(){
    scanf(RIN,&n);
    scanf(RIN,&G);
    if(G==mod){
        printf("0\n");
        return 0;
    }
    for(int k=1;k<=4;k++){
        init(g[k]);
        for(int i=1;i*i<=n;i++){
            if(n%i) continue;
            sum[k]=(sum[k]+lucas(n,i,g[k]))%g[k];
            if(i*i!=n) sum[k]=(sum[k]+lucas(n,n/i,g[k]))%g[k];
        }
    }
    ll M=mod-1,Mi,Ri;
    for(int i=1;i<=4;i++){
        Mi=M/g[i];
        Ri=inverse(Mi,g[i]);
        sum[0]=(sum[0]+sum[i]*Mi%M*Ri%M)%M;
    }
    printf(RIN"\n",power(G,sum[0]));
    return 0;
}