[egin{aligned}
ans&=sumlimits_{i=1}^ngcd(i,n)^xoperatorname{lcm}(i,n)^y\
&=sumlimits_{i=1}^n(in)^ygcd(i,n)^{x-y}\
&=n^ysumlimits_{d|n}d^xsumlimits_{t|frac nd}mu(t)t^ysumlimits_{i=1}^{frac n{dt}}i^y\
&=frac{n^y}{y+1}sumlimits_{d|n}d^xsumlimits_{t|frac nd}mu(t)t^ysumlimits_{k=0}^y{y+1choose k}B^+_k(frac n{dt})^{y-k+1}\
&=sumlimits_{k=0}^y{y+1choose k}B^+_k(id_x imes(mucdot id_y) imes id_{y-k+1})(n)
end{aligned}
]
(id_x imes(mucdot id_y) imes id_{y-k+1})是一个积性函数，因此我们只需要关心(id_x imes(mucdot id_y) imes id_{y-k+1}(p^e))的值。
我们知道(forall ein[2,+infty),(mucdot id_y)(p^e)=0)，那么只需要对((mucdot id_y)(1),(mucdot id_y)(p))的情况分类讨论，剩下的暴力卷积即可。
#include<map>
#include<cctype>
#include<cstdio>
#include<random>
#include<cstring>
#include<algorithm>
using i64=long long;
const int N=3007,P=1000000007;
i64 fac[N],ifac[N],B[N];
i64 read(){i64 x=0;char c=getchar();while(isspace(c))c=getchar();while(isdigit(c))(x*=10)+=c&15,c=getchar();return x;}
void inc(i64&a,i64 b){a+=b-P,a+=a>>63&P;}
void dec(i64&a,i64 b){a-=b,a+=a>>63&P;}
i64 pow(i64 a,i64 b){i64 r=1;for(a%=P;b;b>>=1,a=a*a%P)if(b&1)r=r*a%P;return r;}
i64 C(int n,int m){return fac[n]*ifac[m]%P*ifac[n-m]%P;}
void init(int n)
{
fac[0]=ifac[0]=B[0]=1;
for(int i=1;i<=n;++i) ifac[i]=pow(fac[i]=i*fac[i-1]%P,P-2);
for(int i=1;i<=n;++i)
if(i==1||~i&1)
{
for(int k=0;k<i;++k) dec(B[i],C(i+1,k)*B[k]%P);
B[i]=B[i]*pow(C(i+1,i),P-2)%P;
}
B[1]=P-B[1];
}
namespace Factoring
{
std::map<i64,int>cnt;std::mt19937 rng(19260817);
i64 mul(i64 a,i64 b,i64 p){return (__int128)a*b%p;}
i64 add(i64 a,i64 b,i64 p){return a+=b-p,a+=a>>63&p;}
i64 pow(i64 a,i64 b,i64 p){i64 r=1;for(;b;b>>=1,a=mul(a,a,p))if(b&1)r=mul(a,r,p);return r;}
int Judge(i64 x)
{
static int p[7]={2,3,5,7,11,61,24151};
for(int i=0;i<7;++i) if(x==p[i]) return 1; else if(!(x%p[i])||pow(p[i],x-1,x)!=1) return 0;
i64 phi=x-1;while(~phi&1) phi>>=1;
for(int i=0;i<7;++i)
{
i64 y=pow(p[i],phi,x);
if(y==1) continue;
while(y!=1&&y!=x-1) y=mul(y,y,x);
if(y==1) return 0;
}
return 1;
}
i64 Find(i64 x)
{
i64 c=rng()%(x-1)+1,las=0,now=0;
for(int len=1;;las=now,len<<=1)
{
i64 prod=1,d;
for(int i=1;i<=len;++i) if(now=mul(now,now,x)+c,prod=mul(prod,llabs(now-las),x),!(i&127)&&(d=std::__gcd(prod,x))!=1) return d;
if((d=std::__gcd(prod,x))!=1) return d;
}
}
void Factor(i64 x,int f)
{
if(x==1) return ;
if(Judge(x)) return cnt[x]+=f,void();
i64 p=Find(x);int c=0;
while(!(x%p))x/=p,++c;
Factor(p,f*c),Factor(x,f);
}
}using Factoring::cnt;
void solve()
{
static i64 f[N],g[N];
i64 n=read(),x=read(),y=read(),ans=0;
cnt.clear(),Factoring::Factor(n,1);
for(int k=0;k<=y;++k)
{
if(!B[k]) continue;
i64 t=C(y+1,k)*B[k]%P,z=y+1-k;
for(auto[p,e]:cnt)
{
i64 t1=0,t2=0,u=f[1]=pow(p,x),v=g[1]=pow(p,z);f[0]=g[0]=1;
for(int j=2;j<=e;++j) f[j]=u*f[j-1]%P,g[j]=v*g[j-1]%P;
for(int j=0;j<=e;++j) (t1+=f[j]*g[e-j])%=P;
for(int j=0;j<e;++j) (t2+=f[j]*g[e-j-1])%=P;
(t2*=P-pow(p,y))%=P,(t*=(t1+t2))%=P;
}
inc(ans,t);
}
printf("%lld
",ans*pow(n,y)%P*pow(y+1,P-2)%P);
}
int main()
{
init(3001);
for(int t=read();t;--t) solve();
}