BZOJ 1951 [Sdoi2010] 古代猪文 数论 Lucas定理

题目大意：给定N,G，求

N,G<=1,000,000,000

不妨设

那么如果求出来x，剩下的就是一步快速幂。问题转化为如何求x。

设p=999911659
我们发现p是一个质数且G不可能为p的倍数，根据费马小定理有

即在mod p的意义下G^(p-1)=1
可以得到
其中

我们还注意到，由于N,k很大，所以不能递推计算，只能用阶乘公式计算，这时需要用到Lucas定理。

然而计算组合数要用到除法，也就要用到逆元，p-1即999911658不是质数，求不了逆元，怎么办？

999911658=2*3*4679*35617

可以把其分解成4个质数，计算得到同余方程组，中国剩余定理合并即可。

#include <cstdio>
#include <cmath>
#define MOD 999911659
#define mod 999911658
using namespace std;
typedef long long LL;
LL N,G;
LL prime[]={2,3,4679,35617},X[4],inv[40000],fac[40000];
LL C(LL n,LL m,LL p) {//Lucas theorem
    if(n<m) return 0;
    if(n<p && m<p) return fac[n]*inv[m]%p*inv[n-m]%p;
    return C(n%p,m%p,p)*C(n/p,m/p,p)%p;
}
LL solve(LL p) {
    //initialize
    fac[0]=1;
    for(int i=1;i<p;i++) fac[i]=fac[i-1]*i%p;
    inv[p-1]=p-1;
    for(int i=p-2;~i;i--) inv[i]=inv[i+1]*(i+1)%p;
    //calculate
    LL lim=sqrt(N),ans=0;
    for(int i=1;i<=lim;i++) {
        if(N%i) continue;
        ans=(ans+C(N,i,p))%p;
        if(N/i==i) continue;
        ans=(ans+C(N,N/i,p))%p;
    }
    return ans;
}
void egcd(LL a,LL b,LL& d,LL& x,LL& y) {
    if(!b) { d=a; x=1; y=0; return ; }
    egcd(b,a%b,d,y,x); y-=x*(a/b);
    return ;
}
LL CRT() {//Chinese_Remainder_Theorem
    LL ans=0;
    for(int i=0;i<4;i++) {
        LL d,x,y;
        egcd(mod/prime[i],prime[i],d,x,y);
        LL cach=(x%mod*(mod/prime[i])%mod+mod)%mod;
        ans+=cach*X[i]%mod;
        ans%=mod;
    }
    return ans;
}
LL f_pow(LL x,LL y) {
    LL ans=1;
    while(y) {
        if(y&1) ans=ans*x%MOD;
        x=x*x%MOD;
        y>>=1;
    }
    return ans;
}
int main() {
    scanf("%lld%lld",&N,&G);
    if(!(G%MOD)) {
        printf("0\n");
        return 0;
    }
    for(int i=0;i<4;i++) X[i]=solve(prime[i]);
    printf("%lld\n",f_pow(G,CRT()));
    return 0;
}