ydc的题解 orzz
因为这道题才去学习了伯努利数 可以看这里
#include<cstdio>
#include<cstdlib>
#include<algorithm>
using namespace std;
typedef long long ll;
namespace RHO{
int prime[9]={2,3,5,7,11,13,17,19,23};
unsigned long long RR;
int Pcnt;
long long p[70];
inline long long R(long long Mo){return (RR+=4179340454199820289ll)%Mo;}
inline long long mul(long long x,long long y,long long Mo){
long long tmp=(x*y-(long long)((long double)x/Mo*y+0.5)*Mo);
return tmp<0?tmp+Mo:tmp;
}
inline long long Pow(long long x,long long k,long long Mo){
long long ans=1;
for(;k;k>>=1,x=mul(x,x,Mo))if(k&1)ans=mul(ans,x,Mo);
return ans;
}
inline bool MR(long long n){
if(n<=1)return false;
for(int i=0;i<9;i++)if(n==prime[i])return true;
long long d=n-1;int tmp=0;
while((d&1)==0)d>>=1,tmp++;
for(int i=0;i<9;i++){
long long x=Pow(prime[i],d,n),p=x;
for(int j=1;j<=tmp;j++){
x=mul(x,x,n);
if((x==1)&&(p!=1)&&(p!=n-1))return false;
p=x;
}
if(x!=1)return false;
}
return true;
}
inline long long f(long long x,long long c,long long Mo){return (mul(x,x,Mo)+c)%Mo;}
inline long long gcd(long long x,long long y){return x==0?y:gcd(y%x,x);}
inline long long check(long long c,long long n){
long long x=R(n),y=f(x,c,n),p=n;
while((x!=y)&&((p==n)||(p==1))){
if(x>y)p=gcd(n,x-y);
else p=gcd(n,y-x);
x=f(x,c,n);y=f(f(y,c,n),c,n);
}
return p;
}
inline void rho(long long n){
if(n<=1)return;
if(MR(n)){p[++Pcnt]=n;return;}
while(true){
long long tmp=check(R(n-1)+1,n);
if(tmp!=n && tmp!=1){rho(tmp),rho(n/tmp);return;}
}
}
}
const int M=3005;
const int P=1e9+7;
int C[M][M],inv[M],B[M];
int T[M];
inline void Pre(int m){
inv[1]=1; for (int i=2;i<=m+1;i++) inv[i]=(ll)(P-P/i)*inv[P%i]%P;
C[0][0]=1;
for (int i=1;i<=m+1;i++){
C[i][0]=1; for (int j=1;j<=i;j++) C[i][j]=(C[i-1][j-1]+C[i-1][j])%P;
}
B[0]=1;
for (int i=1;i<=m;i++){
ll ret=0;
for (int j=0;j<i;j++) ret+=(ll)C[i+1][j]*B[j]%P;
ret%=P;
B[i]=(P-ret*inv[i+1]%P)%P;
}
}
ll n; int x,y;
int Pcnt; ll p[70]; int q[70];
inline ll Pow(ll a,int b){
a%=P; ll ret=1;
for (;b;b>>=1,a=a*a%P) if (b&1) ret=ret*a%P;
return ret;
}
inline ll Inv(ll a){ return Pow(a,P-2); }
inline ll calc(ll p,int t,int x,int k){
if (t==1) return (Pow(p,x)+Pow(p,k))%P;
if (t==0) return 1;
if (x==k) return Pow(p,t*k)*(t+1)%P;
ll tt=x>k?Pow(p,x-k):Inv(Pow(p,k-x));
return Pow(p,t*k)*(Pow(tt,t+1)-1)%P*Inv(tt-1)%P;
/*ll ret=0;
for (int a=0;a<=t;a++)
ret+=Pow(p,a*x+(t-a)*k);
return ret%P;*/
}
inline ll Solve(int K){
ll ret=1;
for (int i=1;i<=Pcnt;i++){
ll tmp=(calc(p[i],q[i],x,K)+P-Pow(p[i],y)*calc(p[i],q[i]-1,x,K)%P)%P;
ret=ret*tmp%P;
}
return ret;
}
int main(){
int T;
freopen("t.in","r",stdin);
freopen("t.out","w",stdout);
Pre(3000); scanf("%d",&T);
while (T--){
scanf("%lld%d%d",&n,&x,&y);
RHO::Pcnt=0; RHO::rho(n);
Pcnt=RHO::Pcnt; for (int i=1;i<=Pcnt;i++) p[i]=RHO::p[i];
ll tmp=n;
sort(p+1,p+Pcnt+1); Pcnt=unique(p+1,p+Pcnt+1)-p-1;
for (int i=1;i<=Pcnt;i++){
q[i]=0;
while (tmp%p[i]==0) tmp/=p[i],q[i]++;
}
ll ans=0;
for (int k=1;k<=y+1;k++){
ll tk=(ll)inv[y+1]*C[y+1][k]%P*B[y+1-k]%P;
if (k==y) tk=(P-tk)%P;
ans+=tk*Solve(k)%P;
}
ans%=P; ans=ans*Pow(n,y)%P;
printf("%lld\n",ans);
}
return 0;
}