古代猪文 HYSBZ - 1951

题意：

计算：

\[G^{\sum_{d|N}{C(N,d)}}\% p \]

分析：

当 \(G\) 和 \(p\) 不互质时：答案为 \(0\)；
当 \(G\) 和 \(p\) 互质时：
根据欧拉降幂：

\[G^{\sum_{d|N}{C(N,d)}}\% p=G^{\sum_{d|N}{C(N,d)}\%(p-1)}\% p \]

问题转化为求大组合数取模。
  大组合取模而且模数并非素数，可以用拓展 \(Lucas\) 定理，但是模数 \((p-1)\) 远远大于 \(1e5\)，而且拓展 \(Lucas\) 的复杂度比较高，加上要用很多次（因为 \(N\) 的因子有很多），肯定会超时（亲身尝试）。拓展 \(Lucas\) 计算时，主要处理了 \(C(n,m)\% p^t\) 这个问题。但此处，对 \((p-1)\) 质因子分解发现，其因子有：\(2,3,4679,35617\)，其幂均为 \(1\)。那么我们直接处理 \(C(n,m)\%p\)，是很简单的。借助拓展 \(Lucas\) 定理的思想，分别用 \(4\) 个因子当模数来处理组合数，得到 \(4\) 个同余方程，然后用 \(CRT\) 进行合并答案，最后快速幂即可。

代码：

#include <bits/stdc++.h>
using namespace std;
typedef long long ll;
const int maxn=4e4+5;
ll p=999911659;
int fac[5]={0,2,3,4679,35617},r[6];
ll a[maxn];
ll power(ll a,ll b,ll mod)
{
    ll res=1;
    a%=mod;
    while(b)
    {
        if(b&1)
            res=res*a%mod;
        a=a*a%mod;
        b>>=1;
    }
    return res%mod;
}
ll exgcd(ll a,ll b,ll &x,ll &y)
{
    if(b==0)
    {
        x=1,y=0;
        return a;
    }
    ll gcd=exgcd(b,a%b,x,y);
    ll tmp=x;
    x=y;
    y=tmp-(a/b)*y;
    return gcd;
}
ll get_inv(ll a,ll b)
{
    ll x,y;
    ll gcd=exgcd(a,b,x,y);
    x=(x%b+b)%b;
    return x;
}
void init(ll mod)//预处理出阶乘取模
{
    a[0]=1;//0!=1，注意，一开始写错了
    for(int i=1;i<=mod;i++)
        a[i]=a[i-1]*1LL*i%mod;
}
ll CRT()
{
    ll M=1,ans=0;
    for(int i=1;i<=4;i++)
        M*=fac[i];
    for(int i=1;i<=4;i++)
    {
        ll t=M/fac[i];
        ll x=get_inv(t,1LL*fac[i]);
        ans=(ans+t*x*r[i])%M;
    }
    return ans;
}
ll C(ll n,ll m,ll mod)
{
    if(m>n) return 0;
    m=min(m,n-m);
    ll tn=a[n];
    ll tm=a[m]*a[n-m]%mod;
    return tn*get_inv(tm,mod)%mod;
}
ll Lucas(ll n,ll m,ll mod)
{
    return m?C(n%mod,m%mod,mod)*Lucas(n/mod,m/mod,mod)%mod:1;
}
ll solve(ll n,ll m)
{
    ll res=0;
    for(int j=1;j<=4;j++)
    {
        init(1LL*fac[j]);
        for(ll i=1;i*i<=n;i++)
        {
            if(n%i==0)
            {
                r[j]=(r[j]+Lucas(n,i,fac[j]))%fac[j];
                if((n/i)!=i)
                    r[j]=(r[j]+Lucas(n,n/i,fac[j]))%fac[j];
            }
        }
    }
    ll ans=CRT();
    return power(m,ans,p);
}
int main()
{
    ll n,g;
    scanf("%lld%lld",&n,&g);
    if(__gcd(g,p)>1)
    {
        printf("0\n");
        return 0;
    }
    printf("%lld\n",solve(n,g));
    return 0;
}