淦，为什么题面这么复杂啊
一句话题意：
给定整数(q,n(1leq q,n leq 10^9)),求出(q^{sum_{d|n} C_n^d}(mod 999911659))
根据欧拉定理的推论
当(a,p)互质时(a^b=a^{b\% varphi(p)}(mod p))
可知
[q^{sum_{d|n} C_n^d\%999911658}(mod 999911659)
]
现在我们只需要求解出(sum_{d|n}C_n^d\%999911658)即可
我们尝试分解质因数,发现(99911658=2*3*4679*35617),它的所有质因子的指数都是(1)
所有我们可以枚举(n)的约数(d),然后用卢卡斯定理
[C_n^d=C_{n\%p}^{d\%p}C_{n/p}^{d/p}(mod p)
]
求(C_n^d),分别计算出(sum_{d|n}C_n^d)对(2,3,4679,35617)这四个质数取模的结果分别记为(a_1,a_2,a_3,a_4)
可以得出线性同余方程组
[egin{cases}xmod 2=a_1\xmod 3=a_2\x mod 4679=a_3\
xmod35617=a_4end{cases}]
如果看着不习惯可以自己改一下格式,qwq
利用中国剩余定理就诶出最小非负整数解(x),然后利用快速幂求出(q^x)即为答案
/*
@ author:pyyyyyy
-----思路------
-----debug-------
n,q输入的时候反了
*/
#include<vector>
#include<bits/stdc++.h>
#define int long long
using namespace std;
const int mod=999911658;//这个模数是已经用过欧拉定理的推论的
const int N=50111;
int q,n;
int jc[N];
int a[5],b[5]={0,2,3,4679,35617};
void init(int p)
{
jc[0]=1;
for(int i=1;i<=p;++i) jc[i]=(jc[i-1]%p*i)%p;
}
int ksm(int a,int b,int p)
{
int ret=1;
while(b)
{
if(b&1) ret=(ret*a)%p;
a=(a*a)%p;
b>>=1;
}
return ret%p;
}
int exgcd(int a,int b,int &x,int &y)
{
if(b==0) {
x=1;y=0;return a;
}
int d=exgcd(b,a%b,x,y);
int z=x;x=y;y=z-y*(a/b);
return d;
}
int C(int n,int m,int p)
{
if(n<m) return 0;
return (jc[n]%p*ksm(jc[m],p-2,p)%p*ksm(jc[n-m],p-2,p)%p)%p;
}
int lucas(int n,int m,int p)
{
if(n<m) return 0;
if(!n) return 1;
return lucas(n/p,m/p,p)*C(n%p,m%p,p)%p;
}
int M[5];
int crt()
{
int ans=0;
int m=1;
for(int i=1;i<=4;++i) m*=b[i];
for(int i=1;i<=4;++i) M[i]=m/b[i];
for(int i=1;i<=4;++i)
{
int x,y,d;
d=exgcd(M[i],b[i],x,y);
x=(x%b[i]+b[i])%b[i];
ans+=(a[i]%mod*M[i]%mod*x%mod)%mod;
}
return ans%mod;
}
signed main()
{
// freopen(".in","r",stdin);
// freopen(".out","w",stdout);
// cin>>q>>n;
cin>>n>>q;
if(q%(mod+1)==0)
{
cout<<0;
return 0;
}
for(int k=1;k<=4;++k)//分别对应a_k
{
init(b[k]);//预处理出来当前模数下的阶乘
//枚举d|n
// for(int i=1;i<=10;++i) cout<<jc[i]<<' ';
// cout<<endl;
for(int i=1;i*i<=n;++i) {
if(n%i==0)
{
a[k]=(a[k]+lucas(n,i,b[k]))%b[k];
if(i*i!=n) a[k]=(a[k]+lucas(n,n/i,b[k]))%b[k];
// if(i!=sqrt(n)) a[k]+=(lucas(n,i,b[k])%b[k]+lucas(n,n/i,b[k])%b[k])%b[k];
// else if(i==sqrt(n)) a[k]+=lucas(n,i,b[k])%b[k];
// cout<<a[k]<<'
';
}
}
}
// for(int i=1;i<=4;++i) cout<<a[i]<<' ';
int x=crt();
// int x=0;
// for(int i=1;i<=4;i++)
// x=(x+a[i]*(mod/b[i])%mod*ksm(mod/b[i],b[i]-2,b[i]))%mod;
cout<<ksm(q,x,mod+1);
return 0;
}