蒟蒻惊叹于一道小小的数论题竟能涉及这么多知识点!不过,掌握了这些知识点,拿下这道题也并非难事。
题意一行就能写下来:
给定\(N,G\),求\(G^{\sum \limits _{d|N}C(N,d)}(\mod999911659)\)
乍一看,指数这么大,要怎么处理好呢?上费马小定理。
平时用费马小定理求逆元用多了,\(a^{p-2}\equiv inv(a)(\mod p)\),搞得蒟蒻差点忘了它原本的样子\(a^{p-1}=1(\mod p)\),那原式的指数\(\sum \limits _{d|N}C(N,d)\)就可以直接在\(\mod999911658\)意义下求了。
首先枚举\(d|N\),这个还好办,直接\(O(\sqrt{N})\)把约数筛出来就好了,注意特判\(N\)为完全平方数的情况。
紧接着我们就碰到了组合数取模。如果模数是个质数还好办,直接上卢卡斯定理。但是\(999911658\)明显不是,把它放到计算器里点一下fact会出来\(2×3×4679×35617\)。
只能用扩展卢卡斯了。所谓扩展卢卡斯,其实就是对于一个非质数的模数,把它质因数分解,求出组合数模每一个质因数的结果。这样等于说得到了一个同余方程组,中国剩余定理还原答案即可。注意对于每一个质因数都要先重新打一遍阶乘表,再对每一个\(N\)的约数算组合数求和取模。
这样指数就求出来了。最后快速幂得到答案。
#include<cstdio>
#include<cmath>
#define LL long long
#define R register LL
#define add(a,b,p) a=(a+b)%p
const LL YL=999911659,N=4e4,pr[4]={2,3,4679,35617};
LL f[N],fac[N],res[4],x,y,t;
LL qpow(R b,R k,R p){//快速幂求b^k%p
R a=1;
for(;k;k>>=1,b=b*b%p)
if(k&1)a=a*b%p;
return a;
}
LL C(R n,R m,R p){//组合数,注意特判
return n<m?0:fac[n]*qpow(fac[m]*fac[n-m],p-2,p)%p;
}
LL lucas(R n,R m,R p){//同样注意特判
return m?C(n%p,m%p,p)*lucas(n/p,m/p,p)%p:1;
}
void exgcd(R a,R b){
if(b)exgcd(b,a%b),t=x,x=y,y=t-a/b*y;
}
int main(){
fac[0]=1;
R n,g,m,p=0,i,j,ans=0;
scanf("%lld%lld",&n,&g);
if(!(g%YL))return puts("0"),0;//又一次注意特判
m=sqrt(n);
for(i=1;i<=m;++i)
if(!(n%i))f[++p]=i,f[++p]=n/i;
p-=f[p-1]==f[p];//还是注意特判,防止被算两次
for(i=0;i<4;++i){
for(j=1;j<pr[i];++j)
fac[j]=fac[j-1]*j%pr[i];
for(j=1;j<=p;++j)
add(res[i],lucas(n,f[j],pr[i]),pr[i]);
//下面中国剩余定理,注意x变成正数
x=1;y=0;exgcd(m=(YL-1)/pr[i],pr[i]);
add(ans,res[i]*(x%pr[i]+pr[i])*m,(YL-1));
}
printf("%lld\n",qpow(g,ans,YL));
return 0;
}