数论【Sdoi2010】 古代猪文

题目大意：
给定N和G

求G^x mod 999911659 的值

题目分析：
step1：
观察第一个式子，k可以在sqrt（N）的时间内枚举出来；

step2：
但是这个x求出来必然时巨大的（因为N是1e9级别的），
根据费马小定理，a^(p-1)≡1（mod p）(gcd(a,p)==1)
因为999911659是质数，所以当G不是它的倍数的时候必然是与它互质的，那么G^x除以若干个G^999911658在模p的意义下是完全相同的，问题转化成了G^（x mod 999911658）；

step3：
接下来可以考虑组合数的问题了，组合数要是递推来算时间和空间都接受不了，那就只能采用阶乘公式来算了，C（n，m）=n！/（m!*（n-m）!）；

step4：
算组合数需要预处理出N以内的阶乘和它们的逆元，但是N很大，时间和空间都接受不了，需要用到卢卡斯定理：C（n，m）=C（n/p，m/p）*C（n mod p，m mod p）（p=999911658）；

step5：可以发现即使这样做，上一步骤中的p仍然很大，但是因为这个p不是质数，999911658=2* 3* 4679* 356717，可以用中国剩余定理把我们要求的x拆成几个同余方程，再合并出最终的x（step4中的东西就要算4遍，但是是可以接受的）；

step6：
最后一步轻松的快速幂。

时间复杂度：O（4*sqrt（N）+log（N））
（这道题时间复杂度我不太会算，大概是这样，枚举k要sqrt（N）的时间，求逆元（中国剩余定理）和快速幂是log级别的，卢卡斯定理那个比较玄学，但是应该不会花费太多时间，总之稳稳的能过就是）

注意事项：
1、卢卡斯定理最好用循环写不要递归（其实我也是照着版子打的，但是我的某同学用递归写wa了）；
2、当G是999911659的倍数是要特判，直接输出0；
3、最好用long long类型存储变量；
4、做预处理的时候其实可以发现（p-1）！≡p-1（mod p），并且（p-1）！的逆元也恰好等于p-1，所以在预处理的时候可以不必求一次逆元，直接赋值往回递推也可以。

代码如下：

#include<cstdio>
#include<cstdlib>
#include<algorithm>
#include<iostream>
#include<string>
#include<cstring>
#include<cmath>
//999911658=2*3*4679*35617
#define LL long long
#define mod 999911659
using namespace std;
LL N,G,e;
LL x[4];
LL pre[4][40000];
LL ni[4][40000];
LL zs[4]={2,3,4679,35617};
LL gcd(LL a,LL b,LL &x,LL &y)//求逆元
{
    if(b==0) {x=1;y=0;}
    else     {gcd(b,a%b,y,x); y-=(a/b)*x;}
}
LL lucas(LL p,LL q,LL o)//卢卡斯定理
{
    LL sum=1;
    while(p && q)
    {
        LL x=p%zs[o],y=q%zs[o];
        if(y>x) return 0;
        sum=(sum*(pre[o][x]*(ni[o][x-y]*ni[o][y])%zs[o])%zs[o])%zs[o];
        p/=zs[o];
        q/=zs[o];
    }
    return sum;
}
int main()
{
    scanf("%lld%lld",&N,&G);
    if(G%mod==0) {printf("0\n");return 0;}
    for(int o=0;o<4;o++)//预处理
    {
        pre[o][0]=1;ni[o][0]=1; 
        for(int i=1;i<=zs[o];i++)
        pre[o][i]=(pre[o][i-1]*i)%zs[o];
        ni[o][zs[o]-1]=zs[o]-1;
        for(int i=zs[o]-2;i>=1;i--)
        ni[o][i]=(ni[o][i+1]*(i+1))%zs[o];
        for(LL i=1;i<=sqrt(N);i++)//统计组合数
        if(N%i==0)
        {
            x[o]=(x[o]+lucas(N,i,o))%zs[o];
            if(i*i<N) x[o]=(x[o]+lucas(N,N/i,o))%zs[o];
        }
    };
    for(int i=0;i<4;i++)//中国剩余定理
    {
        LL o,p;
        LL M=(mod-1)/zs[i];
        gcd(M,zs[i],o,p);
        e=(e+x[i]*o*M)%(mod-1);
    }
    while(e<=0) e+=mod-1;
    LL ans=1,a=G;
    while(e)
    {
        if(e&1) ans=(ans*a)%mod;
        a=(a*a)%mod;
        e>>=1;
    }
    printf("%lld\n",ans);
    return 0;
}