bzoj 1951 [Sdoi2010]古代猪文（数论知识）

【题目链接】

　　http://www.lydsy.com/JudgeOnline/problem.php?id=1951

【思路】

       一道优(e)秀(xin)的数论题。

       首先我们要求的是(G^sigma{ C(n,n/i),i|n })%P，即G^M %P，根据费马小定理G^(P-1) ≡1(mod P),我们要求的就是G^(M%(P-1)) %P。

       考虑C(n,i)%(P-1)，由于n i P都比较大所以不好求组合数。发现P-1可以分解质因数为2,3,4679,35617，将C(n,i)对每一个质因子取模，会得到一个形为x≡ ai(mod pi)的模线性方程组，可以用中国剩余定理确定x。对于C(n,i)%p，此时p比较小，我们可以用lucas定理求解。

       总的来说就是先用O(sqrt(n))的时间枚举约数，然后用lucas定理求出不同模数下的ai，最后联立方程组，中国剩余定理解。

       注意当(G,P)!=1的时候费马小定理不成立，此时答案为0。

       关于lucas的写法，a^p-2 %p是a在模p下的逆，因为a^(p-2) *a=a^(p-1)，由费马小定理得a^(p-1) %p=1，因此p必须满足为质数才能使用这种方法。

【定理】

1—  欧拉定理

当a与p互质时，a^phi(p) mod p=1 

费马小定理即欧拉定理在p为质数时的特例， a^(p-1)  ≡1 mod p

2—  Lucas定理

C(n,m)%p=C(n/p,m/p)*C(n%p,m%p)

【代码】

#include<cmath>
#include<cstdio>
#include<cstring>
#include<iostream>
using namespace std; 

typedef long long LL;
const int N = 36000;
const LL mod[5]={2,3,4679,35617,999911659}; 

LL fac[4][N],a[4],n,G;

void gcd(LL a,LL b,LL &x,LL& y) {
    if(!b) { x=1;y=0; }
    else gcd(b,a%b,y,x) , y-=x*(a/b);
}
LL pow(LL x,LL p,LL MOD) {
    LL ans=1;
    while(p) {
        if(p&1) ans=(ans*x)%MOD;
        x=(x*x)%MOD; p>>=1;
    }
    return ans;
}
LL C(LL n,LL m,int x) {
    if(n<m) return 0;
    return (fac[x][n]*pow(fac[x][n-m]*fac[x][m],mod[x]-2,mod[x]))%mod[x];
}
LL lucas(LL n,LL m,int x) {
    if(!m) return 1;
    return lucas(n/mod[x],m/mod[x],x)*C(n%mod[x],m%mod[x],x)%mod[x];
}
LL china() {
    LL M=1,d,x,y,ans=0;
    for(int i=0;i<4;i++) M*=mod[i];
    for(int i=0;i<4;i++) {
        d=M/mod[i];
        gcd(d,mod[i],x,y);
        ans=(ans+d*x*a[i])%M;
    }
    while(ans<=0) ans+=M;
    ans=(ans+M)%M;
    return ans;
}
void get_fac() {
    for(int i=0;i<4;i++) {
        fac[i][0]=1;
        for(int j=1;j<=mod[i];j++)
            fac[i][j]=(fac[i][j-1]*j)%mod[i];
    }
}
int main() {
    get_fac();
    scanf("%lld%lld",&n,&G);
    G%=mod[4];
    if(!G) { puts("0"); return 0; }
    for(int i=1;i*i<=n;i++) if(n%i==0) {
        LL tmp=n/i;
        for(int j=0;j<4;j++) {
            a[j]=(a[j]+lucas(n,tmp,j))%mod[j];
            if(tmp!=i) a[j]=(a[j]+lucas(n,i,j))%mod[j];
        }
    }
    printf("%lld",pow(G,china(),mod[4]));
    
    return 0;
}