[BZOJ1951][SDOI2010]古代猪文(Lucas+CRT)

简化题意

  给定$n,g(n,g\leq10^9)$求$g^{\displaystyle\sum_{d|n}C_{n}^{d}} \mod 999911659$

分析

  首先我们发现，指数的那个东西，我们可以通过$O(\sqrt n)$的时间来枚举出所有的约数，然后直接算组合数。
  但因为$n$非常大，所以算组合数的时候我们还需要$Lucas$定理来求组合数，即：
$C_{n}^{m}\mod p=C_{n\mod p}^{m \mod p}\times C_{n\div p}^{m\div p}\mod p$。
  之后我们就发现一个问题：这个指数非常大，我们无法直接求出然后进行快速幂，这时候我们就需要用欧拉降幂定理了，即：
$a^b\mod p=a^{b\mod \phi(p)+\phi(p)}\mod p$。
  又因为当$p$为质数时：
$a^{p-1}\equiv1\mod p,\phi(p)=p-1$。
  所以当$p$为质数时，我们可以得到：
$a^{b\mod \phi(p)+\phi(p)}=a^{b\mod \phi(p)}*a^{\phi(p)}=a^{b\mod p-1}*a^{p-1}=a^{b\mod p-1}$。
  又注意到在这个地方$mod-1$不是质数：$999911658=2*3*4679*35617$。
  所以我们要对4个质数分别计算出答案后通过CRT来合并，这样就做完了。

Code

#pragma GCC optimize(3)
#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
bool Finish_read;
template<class T>inline void read(T &x){Finish_read=0;x=0;int f=1;char ch=getchar();while(!isdigit(ch)){if(ch=='-')f=-1;if(ch==EOF)return;ch=getchar();}while(isdigit(ch))x=x*10+ch-'0',ch=getchar();x*=f;Finish_read=1;}
template<class T>inline void print(T x){if(x/10!=0)print(x/10);putchar(x%10+'0');}
template<class T>inline void writeln(T x){if(x<0)putchar('-');x=abs(x);print(x);putchar('\n');}
template<class T>inline void write(T x){if(x<0)putchar('-');x=abs(x);print(x);}
/*================Header Template==============*/
const int maxn=100005;
const ll mod[4]={2,3,4679,35617},modx=999911659;
namespace {
    inline ll Add(ll a,ll b,ll p) {
        ll res=a+b;
        if(res>=p)
            res-=p;
        return res;
    }
    inline ll Sub(ll a,ll b,ll p) {
        ll res=a-b;
        if(res<0)
            res+=p;
        return res;
    }
    inline ll Mul(ll a,ll b,ll p) {
        return a*b%p;
    }
    inline ll Pow(ll a,ll b,ll p) { 
        ll res=1;
        while(b) {
            if(b&1)
                res=Mul(res,a,p);
            a=Mul(a,a,p);
            b>>=1;
        }
        return res;
    }
    inline ll Inv(ll a,ll p) { // Assert p is prime
        return Pow(a,p-2,p);
    }
    inline ll Div(ll a,ll b,ll p) { // Assert p is prime
        return Mul(a,Inv(b,p),p);
    }
    inline ll Fix(ll a,ll p) {
        if(a<0)
            a+=p;
        if(a>=p)
            a-=p;
        return a;
    }
}
ll fac[4][maxn],ifac[4][maxn],n,g,sum[4];
inline ll C(ll a,ll b,int k) {
    if(a<b)
        return 0;
    return Mul(fac[k][a],Mul(ifac[k][b],ifac[k][a-b],mod[k]),mod[k]);
}
inline ll Lucas(ll a,ll b,int k) {
//  cout<<a<<" "<<b<<" "<<k<<endl;
    if(b==0)
        return 1;
    return Mul(C(a%mod[k],b%mod[k],k),Lucas(a/mod[k],b/mod[k],k),mod[k]);
}
inline void exgcd(ll a,ll b,ll &x,ll &y) {
    if(!b) {
        x=1,y=0;
        return;
    }
    exgcd(b,a%b,y,x);
    y-=a/b*x;
}
inline void Work(int x) {
    for(int i=0;i<4;++i)
        sum[i]=Add(sum[i],Lucas(n,x,i),mod[i]);
}
inline ll CRT() {
    int a1,b1,a2,b2,a,b,c;
    ll x,y;
    a1=mod[0],b1=sum[0];
    for(int i=1;i<4;++i) {
        a2=mod[i],b2=sum[i];
        a=a1,b=a2,c=b2-b1;
        exgcd(a,b,x,y);
        x=Fix(Mul(c,x,b),b);
        b1=b1+a1*x;
        a1=a1*b;
    }
    return b1;
}
int main() {
//  static int N=1e5;
    for(int k=0;k<4;++k) {
        fac[k][0]=1;
        for(int i=1;i<mod[k];++i)
            fac[k][i]=Mul(fac[k][i-1],i,mod[k]);
        ifac[k][mod[k]-1]=Inv(fac[k][mod[k]-1],mod[k]);
        for(int i=mod[k]-2;~i;--i)
            ifac[k][i]=Mul(ifac[k][i+1],i+1,mod[k]);
        for(int i=0;i<mod[k];++i)
            assert(Mul(fac[k][i],ifac[k][i],mod[k])==1);
    }
//  cerr<<Lucas(5,2,3)<<endl;
    read(n),read(g);
    if(g%modx==0)
        return 0*puts("0");
    g%=modx;
//  cout<<n<<" "<<g<<endl;
    for(int i=1;i*i<=n;++i)
        if(n%i==0) {
            Work(i);
            if(i*i!=n)
                Work(n/i);
        }
    printf("%lld\n",Pow(g,CRT(),modx));
}