BZOJ1951 [Sdoi2010]古代猪文

题目大意：给出n,G,求G^P%M，其中M=999911659，P=∑C(n,i) (i=1~n,n%i=0)。

题解：此神题也。注意到M是个素数，则G^P%M=G^(P%phi(M)+phi(M))%M=G^(P%(M-1)+M-1)%M，所以求P%(M-1)即可。错误！前式应满足条件P>=phi(M)即P>=M-1，如果P<M-1，可以随便举出一个反例。所以，判断P是否≥M-1，若不是直接求答案，若是就进行下一步。令人头疼的是M-1不是一个素数，求P%(M-1)就要用Lucas求∑C(n,i) (i=1~n,n%i=0)，其中求C需要乘法逆元，而面对模数M-1不是素数的情况不能直接求。于是，把M-1分解质因数，用Lucas(n,i,p)(p为质因数)求出几个模数后，解一波线性模方程组，用中国剩余定理即可。Lucas(a,b,p)=C(a%p,b%p,p)*Lucas(a/p,b/p,p)%p，C的求法见代码。最后的答案跑快速幂。

代码：

#include<cstdio>
#include<cstring>
#include<cstdlib>
#include<algorithm>
using namespace std;

#define LL long long
#define M 999911659
LL son[]={0,2,3,4679,35617};
#define maxs 36000
LL n,g,fac[maxs][5],p,mo[5];
LL pow_mod(LL a,LL b,LL p)
{
    LL tmp=a%p,ans=1;
    while (b) {if (b&1) ans=ans*tmp%p;tmp=tmp*tmp%p;b>>=1;}
    return ans;
}
void ext_gcd(LL a,LL b,LL& x,LL& y)
{
    if (!b) {x=a;y=0;}
    else ext_gcd(b,a%b,y,x),y-=a/b*x;
}
LL C(LL a,LL b,LL pos)
{
    if (b>a) return 0;if (b==0 || b==a) return 1;LL p=son[pos];
    return fac[a][pos]*pow_mod(fac[b][pos]*fac[a-b][pos],p-2,p)%p;
}
LL Lucas(LL a,LL b,LL pos)
{
    if (!b) return 1;LL p=son[pos];
    return (C(a%p,b%p,pos)*Lucas(a/p,b/p,pos)%p);
}
LL china(LL n,LL* a,LL* m)
{
    LL x,y,ans=0;
    for (int i=1;i<=n;i++)
    {
        LL w=(M-1)/m[i];
        ext_gcd(w,m[i],x,y);
        ans=(ans+w*x*a[i])%(M-1);
    }
    return (ans+(M-1))%(M-1);
}
LL C2(LL a,LL b)
{
    if (b>a) return 0;b=min(a-b,b);if (b>25) return M;
    LL ans=1;
    for (int i=1;i<=b;i++) {ans=ans*(a-i+1)/i;if (ans>=M-1) break;}
    if (ans>=M-1) return M;return ans;
}
int main()
{
    scanf("%lld%lld",&n,&g);
    for (int j=1;j<=4;j++) fac[0][j]=1;
    for (int j=1;j<=4;j++)
        for (int i=1;i<=son[j];i++)
            fac[i][j]=fac[i-1][j]*i%son[j];
    LL tmp=0;int ok=0;
    #define TMP if (tmp>=M-1) {ok=1;break;}
    for (LL i=1;i*i<=n;i++)
        if (!(n%i))
        {
            tmp+=C2(n,i);TMP
            if (i*i!=n) {tmp+=C2(n,n/i);TMP}
        }
    for (int j=1;j<=4;j++)
    {
        mo[j]=0;
        for (LL i=1;i*i<=n;i++)
            if (!(n%i))
            {
                mo[j]=(mo[j]+Lucas(n,i,j))%son[j];
                if (i*i!=n) mo[j]=(mo[j]+Lucas(n,n/i,j))%son[j];
            }
    }
    p=ok?china(4,mo,son)+M-1:tmp;
    printf("%lld",pow_mod(g,p,M)%M);
    return 0;
}