luoguP2480 [SDOI2010]古代猪文

题意

考虑所求即为：(G^{sumlimits_{d|n}C_n^d}\%999911659)。

发现系数很大，先用欧拉定理化简系数：(G^{sumlimits_{d|n}C_n^d\%999911658}\%999911659)。

实际上我们只用求(sumlimits_{d|n}C_n^d\%999911658)，之后快速幂即可。

发现(999911658)不是个质数，没办法用Lucas定理求组合数，于是考虑拆开(999911658)，发现为(2,3,4679,35617)。

于是对模意义下这四个数分别求(sumlimits_{d|n}C_n^d)，假设第(i)个求出的为(a_i)

发现我们得到了四个形如(xequiv a_ipmod{p_i})的方程，用中国剩余定理合并即可得到答案（这其实就是exLucas的简化版）。

code:

#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
const ll mod1=999911659;
const ll mod2=999911658;
const int maxs=40010;
ll n,m;
ll a[5],prime[]={0,2,3,4679,35617};
ll fac[maxs][5],inv[maxs][5];
inline ll power(ll x,ll k,ll mod)
{
	ll res=1;
	while(k)
	{
		if(k&1)res=res*x%mod;
		x=x*x%mod;k>>=1;
	}
	return res;
}
inline ll C(ll n,ll m,ll op)
{
	if(m>n)return 0;
	return fac[n][op]*inv[n-m][op]%prime[op]*inv[m][op]%prime[op];
}
inline ll Lucas(ll n,ll m,ll op)
{
	if(!m)return 1;
	return C(n%prime[op],m%prime[op],op)*Lucas(n/prime[op],m/prime[op],op)%prime[op];
}
void exgcd(ll a,ll b,ll& x,ll& y)
{
	if(!b){x=1,y=0;return;}
	exgcd(b,a%b,x,y);
	ll z=x;x=y,y=z-(a/b)*y;
}
inline ll CRT()
{
	ll res=0;
	for(int i=1;i<=4;i++)
	{
		ll x,y,M=mod2/prime[i];
		exgcd(M,prime[i],x,y);
		x=(x%prime[i]+prime[i])%prime[i];
		res=(res+a[i]*x%mod2*M%mod2)%mod2;
	}
	return res;
}
int main()
{
	scanf("%lld%lld",&n,&m);
	if(m==mod1){puts("0");return 0;}
	for(int i=1;i<=4;i++)
	{
		fac[0][i]=1;
		for(int j=1;j<prime[i];j++)fac[j][i]=fac[j-1][i]*j%prime[i];
		inv[prime[i]-1][i]=power(fac[prime[i]-1][i],prime[i]-2,prime[i]);
		for(int j=prime[i]-1;j;j--)inv[j-1][i]=inv[j][i]*j%prime[i];
	}
	for(ll i=1;i*i<=n;i++)
	{
		if(n%i)continue;
		for(int j=1;j<=4;j++)
		{
			a[j]=(a[j]+Lucas(n,i,j))%prime[j];
			if(i*i!=n)a[j]=(a[j]+Lucas(n,n/i,j))%prime[j];
		}
	}
	printf("%lld",power(m,CRT(),mod1));
	return 0;
}