【BZOJ1951】古代猪文（SDOI2010）-数论大集合

测试地址：古代猪文
做法：本题需要用到的数论知识有：费马小定理，枚举约数，中国剩余定理/合并模线性方程，扩展欧几里得，Lucas定理。可以说是一道很复杂的数论题了。
注意到题目要求的式子是：

G∑k|nCn/knmodp$$G^{\sum _{k|n}{C}_n^{n/k}}\mathrm{mod}p$$

注意到p=999911659$p=999911659$是一个质数，所以根据费马小定理，答案等同于：

G∑k|nCn/knmod(p−1)modp$$G^{\sum _{k|n}{C}_n^{n/k}\mathrm{mod}(p-1)}\mathrm{mod}p$$

所以我们只需计算里面那个和式，然后快速幂算出最后的答案即可。枚举约数k$k$是O(n−−√)$O(\sqrt{n})$的，而计算组合数的话，首先要注意p−1$p-1$不是质数，它质因数分解的结果为2×3×4679×35617$2\times 3\times 4679\times 35617$，我们没有办法，只能分别求出结果对这四个质数取模的取值，然后用中国剩余定理或者合并模线性方程来算出最后的答案，这里我用的是合并模线性方程的做法，即如有两个同余方程：M≡r1(moda1)$M\equiv r_1(\mathrm{mod}{a}_1)$，M≡r2(moda2)$M\equiv r_2(\mathrm{mod}{a}_2)$，那么它们合并之后的结果为：M≡r1+a1x0(moda1a2gcd(a1,a2))$M\equiv r_1+a_1{x}_0(\mathrm{mod}\frac{a_1{a}_2}{gcd(a_1,a_2)})$，其中x0$x_0$为方程a1x1−a2x2=r2−r1$a_1{x}_1-a_2{x}_2=r_2-r_1$的一个解，可以用扩展欧几里得算法求出。
而求组合数对一个质数取模的结果可以用Lucas定理来解决，即Cmn%p=Cm/pn/pCm%pn%p%p$C_n^m\mathrm{\%}p=C_{n/p}^{m/p}{C}_{n\mathrm{\%}p}^{m\mathrm{\%}p}\mathrm{\%}p$。这样我们就解决了这道题，时间复杂度为O(n−−√logn)$O(\sqrt{n}\log n)$。
还有一点要注意，有一个点是G=p$G=p$的，而若是你算出来的幂数模p−1$p-1$的结果是0$0$，那么你会输出1$1$，而正确答案是0$0$，要特判一下。
我傻逼的地方：一大堆公式记错，例如线性求逆元的状态转移方程是inv(i)=(p−p/i)×inv(p%i)%p$inv(i)=(p-p/i)\times inv(p\mathrm{\%}i)\mathrm{\%}p$，以及合并模线性方程的公式也记错了……看来以后要经常背……
以下是本人代码：

#include <bits/stdc++.h>
using namespace std;
typedef long long ll;
int p;
ll n,g,mod[10]={999911659,2,3,4679,35617};
ll fac[40010],inv[40010],fi[40010],ans=0,tot,lasta=1,x;

ll C(ll n,ll m)
{
    if (n<m) return 0;
    if (n>=mod[p]||m>=mod[p])
        return C(n/mod[p],m/mod[p])*C(n%mod[p],m%mod[p])%mod[p];
    return fac[n]*fi[m]%mod[p]*fi[n-m]%mod[p];
}

ll power(ll a,ll b)
{
    ll s=1,ss=a;
    while(b)
    {
        if (b&1) s=(s*ss)%mod[0];
        ss=(ss*ss)%mod[0];b>>=1;
    }
    return s;
}

ll exgcd(ll a,ll b)
{
    ll x0=1,y0=0,x1=0,y1=1;
    while(b)
    {
        ll tmp,q;
        q=a/b;
        tmp=x0,x0=x1,x1=tmp-q*x1;
        tmp=y0,y0=y1,y1=tmp-q*y1;
        tmp=a,a=b,b=tmp%b;
    }
    return x0;
}

int main()
{
    scanf("%lld%lld",&n,&g);

    for(p=1;p<=4;p++)
    {
        fac[0]=fac[1]=inv[0]=inv[1]=fi[0]=fi[1]=1;
        for(ll i=2;i<=mod[p];i++)
        {
            fac[i]=fac[i-1]*i%mod[p];
            inv[i]=(mod[p]-mod[p]/i)*inv[mod[p]%i]%mod[p];
            fi[i]=fi[i-1]*inv[i]%mod[p];
        }
        tot=0;
        for(ll i=1;i*i<=n;i++)
            if (n%i==0)
            {
                tot+=C(n,i);
                if (i*i!=n) tot+=C(n,n/i);
            }
        x=exgcd(lasta,mod[p]);
        x=(x*(tot-ans)%mod[p]+mod[p])%mod[p];
        ans=(ans+x*lasta)%(lasta*mod[p]);
        lasta*=mod[p];
    }
    ans=power(g,ans);
    if (g==mod[0]) printf("0");
    else printf("%lld",ans);

    return 0;
}