bzoj1951 [Sdoi2010]古代猪文（Lucas定理+CRT+逆元）

题目就是求

$$G^{\sum\limits_{k|n}C_n^k}\%P,P=999911659$$
根据欧拉定理，我们就是要求 $\sum\limits_{k|n}C_n^k\%\phi(p)$
因为P是个素数，所以 $phi(p)=p-1$
n很大，所以计算这个组合数求模要用Lucas定理来算

Lucas定理：portal
计算$\binom{n}{m} \mod p$，p是个质数，

$n= n_k * p ^ k + n_{k-1} * p^{k-1}+ ... + n_2 * p^2 + n_1 * p + n_0$

$m=m_k * p ^ k +m_{k-1} * p^{k-1}+ ... +m_2 * p^2 +m_1 * p+m_0$

则$\binom{n}{m} = \prod\limits_{i=0}^k \binom{n_i}{m_i}$

实现：这个形式很像多项式啊变量为p，n和m迭代/=p然后算C(n%p,m%p)就行了，预处理阶乘，逆元之后，复杂度大概$O(log_pn)$

然而我们不是质数怎么办！中国剩余定理！
999911658=2*3*4679*35617
次数都是1次诶
那就不用扩展Lucas了，分别计算一个余数然后用CRT合并一下就好啦。

注意一个小细节：G=999911659!,此时不互质，要特殊处理，也就是输出0直接qaq

#include <cstdio>
#include <cstring>
#include <algorithm>
#include <cmath>
using namespace std;
#define ll long long
#define inf 0x3f3f3f3f
#define N 36000
#define mo 999911659
inline char gc(){
    static char buf[1<<16],*S,*T;
    if(S==T){T=(S=buf)+fread(buf,1,1<<16,stdin);if(T==S) return EOF;}
    return *S++;
}
inline int read(){
    int x=0,f=1;char ch=gc();
    while(ch<'0'||ch>'9'){if(ch=='-')f=-1;ch=gc();}
    while(ch>='0'&&ch<='9') x=x*10+ch-'0',ch=gc();
    return x*f;
}
int fac[N],prime[4]={2,3,4679,35617},inv[N],r[4];
inline int ksm(ll x,int k,int mod){
    ll res=1;for(;k;k>>=1,x=x*x%mod) if(k&1) res=res*x%mod;return res;
}
inline int C(int n,int m,int mod){
    if(m>n) return 0;
    return fac[n]*inv[m]%mod*inv[n-m]%mod;
}
inline int Lucas(int n,int m,int mod){
    int res=1;
    for(;m;n/=mod,m/=mod) res*=C(n%mod,m%mod,mod),res%=mod;return res;
}
inline int solve(int n,int mod){
    fac[0]=1;inv[1]=1;inv[0]=1;
    for(int i=2;i<mod;++i) inv[i]=inv[mod%i]*(mod-mod/i)%mod;
    for(int i=1;i<mod;++i) fac[i]=(ll)fac[i-1]*i%mod;
    for(int i=1;i<mod;++i) inv[i]=(ll)inv[i-1]*inv[i]%mod;
    int x=1,res=0;
    for(x;x*x<n;++x) if(n%x==0) res+=Lucas(n,x,mod),res+=Lucas(n,n/x,mod),res%=mod;
    if(x*x==n) res+=Lucas(n,x,mod),res%=mod;
    return res;
}
inline int CRT(int mod){
    int res=0;
    for(int i=0;i<4;++i){
        ll a=mod/prime[i];int b=prime[i];
        res+=r[i]*a*ksm(a,b-2,b)%mod;res%=mod;
    }return res;
}
int main(){
//  freopen("a.in","r",stdin);
    int n=read(),G=read();if(G==mo){puts("0");return 0;}
    for(int i=0;i<4;++i) r[i]=solve(n,prime[i]);
    printf("%d\n",ksm(G,CRT(mo-1),mo));
    return 0;
}