题目链接：https://www.luogu.com.cn/problem/P4464
简记$gcd(x,y)=(x,y)$。
推式子：
$sum_{i=1}^{n}{(i,n)^xlcm(i,n)^y}$
$=sum_{i=1}^{n}{(i,n)^{x-y}(in)^y}$
$=n^ysum_{d|n}d^{x-y}sum_{i}i^y[(i,n)=d]$
$=n^ysum_{d|n}{d^{x-y}sum_{i=1}^{frac{n}{d}}{(id)^y[(i,frac{n}{d})=1]}}$
$=n^ysum_{d|n}{d^xsum_{i=1}^{frac{n}{d}}{i^y[(i,frac{n}{d})=1]}}$
$=n^ysum_{d|n}{d^xsum_{i=1}^{frac{n}{d}}{i^ysum_{k|(i,frac{n}{d})}{mu(k)}}}$
改变最后的求和式的顺序，注意k的取值要整除$frac{n}{d}$。
$=n^ysum_{d|n}{d^xsum_{k|frac{n}{d}}{mu(k)}sum_{i=1}^{frac{n}{dk}}{(ik)^y}}$
$=n^ysum_{d|n}{d^xsum_{k|frac{n}{d}}{mu(k)k^y}sum_{i=1}^{frac{n}{dk}}{i^y}}$
若我们把最后的和式看成是一个y+1次的多项式，它的第i项系数为$t_i$，那么
$=n^ysum_{j=1}^{y+1}{t_i}sum_{d|n}{d^xsum_{k|frac{n}{d}}{mu(k)k^y}(frac{n}{dk})^j}$
（注意最后一项的变化）
注意到后面所有的形式都是积性函数的卷积，因此我们尝试一个质数的若干次方的贡献，最后将其相乘就得到某一个对应的j的答案。
于是对于固定的$d=p^q,n=p^w$，最后一个和式的贡献为$(p^{w-q})^j-p^y(p^{w-q-1})^j$（分别考虑$d=1$和$d=p$）。特别地，若$d=n$，那么贡献为1。如果我们能快速得到所有$t_i$，那么就能在$O(n^{0.25}logn+ylog^3n)$的复杂度完成一次询问。
考虑如何得到$t_i$。一种方法是使用插值，这里讨论使用伯努利数的方法（应用，证明看https://www.luogu.com.cn/blog/chihik/bo-nu-li-shu）。
递归定义：$sum_{i=0}^{n}{ binom{n+1}{i}B_i=[n=0]}$
令$S_m(n)=sum_{i=0}^{n-1}{i^m}$，那么$S_m(n)=frac{1}{m+1}sum_{i=0}^{m}{ binom{m+1}{i}B_in^{m+1-i}}$
于是对于上面的y+1次多项式，有$S_y(frac{n}{dk})=sum_{i=0}^{frac{n}{dk}-1}{i^y}=frac{1}{y+1}sum_{i=0}^{y}{ binom{y+1}{i}{B_i(frac{n}{dk})^{y+1-i}}}$，所以$t_{y+1-i}=frac{1}{y+1} binom{y+1}{i}B_i$。（它的零次项系数为0！（这是感叹号，不是阶乘））
但我们注意到少了一项$(frac{n}{dk})^y$，因此$t_y$需要增加1（对应了$i=1$）。
1 #include<bits/stdc++.h> 2 using namespace std; 3 typedef long long int ll; 4 typedef long double ld; 5 namespace MATH 6 { 7 inline ll mul(ll x,ll y,ll mod) 8 { 9 return ((x*y-(ll)((ld)x*y/mod)*mod)%mod+mod)%mod; 10 } 11 inline ll qpow(ll x,ll y,ll mod) 12 { 13 ll ans=1,base=x; 14 while(y) 15 { 16 if(y&1) 17 ans=mul(ans,base,mod); 18 base=mul(base,base,mod); 19 y>>=1; 20 } 21 return ans; 22 } 23 const int prime[10]={2,3,5,7,11,13,17,19,23,29}; 24 const int LEN=10; 25 inline bool miller(ll x) 26 { 27 for(int i=0;i<LEN;++i) 28 { 29 if(x<prime[i]) 30 return false; 31 else if(x==prime[i]) 32 return true; 33 ll now=x-1,y=qpow(prime[i],x-1,x); 34 if(y!=1) 35 return false; 36 while(now%2==0&&y==1) 37 { 38 now>>=1; 39 y=qpow(prime[i],now,x); 40 if(y!=1&&y!=x-1) 41 return false; 42 } 43 } 44 return true; 45 } 46 ll gcd(ll x,ll y) 47 { 48 return x%y==0?y:gcd(y,x%y); 49 } 50 inline ll f(ll x,ll y,ll mod) 51 { 52 return (mul(x,x,mod)+y)%mod; 53 } 54 inline ll get(ll n,ll c,int steps) 55 { 56 if(!(n&1)) 57 return 2; 58 ll x=2,y=2,d=1; 59 while(true) 60 { 61 ll tmpx=x,tmpy=y; 62 for(int i=1;i<=steps;++i) 63 { 64 x=f(x,c,n); 65 y=f(y,c,n); 66 y=f(y,c,n); 67 d=mul(d,((y-x)%n+n)%n,n); 68 } 69 d=gcd(d,n); 70 if(d==1) 71 continue; 72 if(d!=n) 73 return d; 74 x=tmpx,y=tmpy; 75 for(int i=1;i<=steps;++i) 76 { 77 x=f(x,c,n); 78 y=f(y,c,n); 79 y=f(y,c,n); 80 d=gcd(((y-x)%n+n)%n,n); 81 if(d!=1&&d!=n) 82 return d; 83 } 84 return 0; 85 } 86 } 87 vector<ll>wait; 88 void pollard(ll n) 89 { 90 if(miller(n)) 91 { 92 wait.push_back(n); 93 return; 94 } 95 ll now=0,c=1,g=pow(n,0.1)+1; 96 while(!now) 97 now=get(n,++c,g); 98 pollard(now),pollard(n/now); 99 } 100 } 101 const ll mod=1000000007; 102 ll C[3005][3005],B[3005]; 103 inline ll qpow(ll x,ll y) 104 { 105 ll ans=1,base=x%mod; 106 while(y) 107 { 108 if(y&1) 109 ans=ans*base%mod; 110 base=base*base%mod; 111 y>>=1; 112 } 113 return ans; 114 } 115 inline void init() 116 { 117 for(int i=0;i<=3004;++i) 118 { 119 C[i][0]=1; 120 for(int j=1;j<=i;++j) 121 C[i][j]=(C[i-1][j-1]+C[i-1][j])%mod; 122 } 123 B[0]=1; 124 for(int i=1;i<=3002;++i) 125 { 126 for(int j=0;j<i;++j) 127 B[i]=(B[i]+B[j]*C[i+1][j])%mod; 128 B[i]=(B[i]*qpow(-i-1,mod-2)%mod+mod)%mod; 129 } 130 } 131 ll tmp[3005]; 132 inline void solve() 133 { 134 ll n,x,y; 135 cin>>n>>x>>y; 136 MATH::wait.clear(); 137 MATH::pollard(n); 138 vector<ll>wait=MATH::wait; 139 map<ll,int>vis; 140 for(ll x:wait) 141 ++vis[x]; 142 ll ans=0; 143 for(int i=0;i<=y;++i) 144 tmp[y+1-i]=C[y+1][i]*B[i]%mod*qpow(y+1,mod-2)%mod; 145 tmp[y]+=1; 146 for(int i=1;i<=y+1;++i) 147 { 148 ll s=1; 149 for(auto A:vis) 150 { 151 ll p=A.first,k=A.second; 152 ll g=0; 153 ll now=1,delta=qpow(p,x); 154 ll G=qpow(p,y); 155 for(int j=0;j<k;++j,now=now*delta%mod) 156 g=(g+now*(qpow(qpow(p,k-j),i)-G*qpow(qpow(p,k-j-1),i)%mod)%mod)%mod; 157 g=(g+now%mod)%mod; 158 s=s*g%mod; 159 } 160 ans=(ans+tmp[i]*s)%mod; 161 } 162 ans=ans*qpow(n%mod,y)%mod; 163 ans=(ans%mod+mod)%mod; 164 cout<<ans<<endl; 165 } 166 int main() 167 { 168 ios::sync_with_stdio(false); 169 init(); 170 int T; 171 cin>>T; 172 while(T--) 173 solve(); 174 return 0; 175 }