题意即求下式:
m∑x|nn!(x!)nx(nx!)m∑x|nn!(x!)nx(nx!)m∑x|nn!(x!)nx(nx!)
m∑x|nn!(x!)nx(nx!)
根据欧拉定理: 对于互质的正整数 a 和 n ,有 aφ(n) ≡ 1 mod n 可知,求指数式子模1e9-402的值然后快速幂即可。
O(sqrt(n))分解质因数求值,注意到n非常大,无法预处理阶乘及逆元。
定理: (p-1)!%p=-1
对于 n!,将p提取出并记下p的指数k,fac(n)=(-1)^(n/p) * (n%p)! * fac(n/p)
做除法时求逆元并判断分子和分母中p的指数能否相消,不能为0,能则为fac(n)*inv( fac(x)^(n/x) )*fac(n/x)
(注:fac(n)表示n的阶乘,inv(n)表示n的逆元)
代码:
#include<algorithm>
#include<iostream>
#include<cstdlib>
#include<cstring>
#include<cstdio>
#include<cmath>
#include<queue>
#include<stack>
#include<map>
#include<set>
#define ll long long
#define mod 999999599
#define Mod 999999598
#define N 200200
using namespace std;
/*
2 1
13 1
5281 1
7283 1
*/
int n,m,cas;
int MOD[5];
ll ret[5],ans;
ll pre[N][5],inv[N][5];
int num[N][5];
ll qm(ll a,ll b,ll md)
{ll ret=1,t=a%md;while(b){if(b&1)ret=ret*t%md;t=t*t%md;b>>=1;}return ret;
}
void ex_gcd(int a,int b,int &x,int &y)
{if(b==0){x=1,y=0;return ;}ex_gcd(b,a%b,x,y);int t=x;x=y,y=t-a/b*y;
}
void init()
{MOD[1]=2,MOD[2]=13;MOD[3]=5281,MOD[4]=7283;pre[0][1]=pre[0][2]=pre[0][3]=pre[0][4]=1;inv[0][1]=inv[0][2]=inv[0][3]=inv[0][4]=1;for(int i=1;i<=4;i++)for(int j=1;j<MOD[i];j++){int t=j,cnt=0;while(t%MOD[i]==0)t/=MOD[i],cnt++;pre[j][i]=pre[j-1][i]*t%MOD[i];num[j][i]=num[j-1][i]+cnt;inv[j][i]=qm(pre[j][i],MOD[i]-2,MOD[i]);}
}
ll fac(int x,int p)
{if(x<MOD[p])return pre[x][p];ll t=x/MOD[p],ret=1;ret=((t&1)?-1:1);return ret*pre[x%MOD[p]][p]%MOD[p]*fac(x/MOD[p],p)%MOD[p];
}
int con(int x,int p)
{if(x<MOD[p])return 0;int ret=x/MOD[p];return con(x/MOD[p],p)+ret;
}
void cal(int x)
{//ret[1]=1;for(int i=1;i<=4;i++){//if(n>=MOD[i])continue;ll tmp,txp,now;tmp=fac(x,i),now=con(x,i);tmp=qm(tmp,n/x,MOD[i]),now=now*(n/x);tmp=qm(tmp,MOD[i]-2,MOD[i]),now=-now;txp=fac(n/x,i),txp=qm(txp,MOD[i]-2,MOD[i]);tmp=tmp*txp%MOD[i],now-=con(n/x,i);tmp=tmp*fac(n,i)%MOD[i],now+=con(n,i);if(now!=0)tmp=0;ret[i]=(ret[i]+tmp)%MOD[i];}
}
void CRT()
{for(int i=1;i<=4;i++){ll M=Mod/MOD[i];//int Mo,y;ex_gcd(M,MOD[i],Mo,y);ll Mo=qm(M,MOD[i]-2,MOD[i]);ans=(ans+M*Mo%Mod*ret[i]%Mod)%Mod;}ans=ans+Mod;ans=qm(m,ans,mod);
}
void solve()
{int x=sqrt(n);for(int i=1;i<=x;i++){if(n%i==0){int y=n/i;cal(i);if(y!=i)cal(y);}}CRT();printf("%lld\n",ans);
}
int main()
{scanf("%d",&cas);init();while(cas--){memset(ret,0,sizeof(ret));scanf("%d%d",&n,&m);ans=0;solve();}
}