P4464 JZPKIL 题解
又是一道独立（基本上是）做出的黑，好耶！注：下文为了简洁除法用 / 代替分数。
前置知识：伯努利数求自然数幂和。
伯努利数：\(B_0=1,\sum\limits_{i=0}^n\dbinom{n+1}{i}B_i=0(n\ge1)\)，于是可以 \(O(n^2)\) 预处理前 \(n\) 个伯努利数。
有结论：\(\sum\limits_{i=1}^{n-1}i^k=\dfrac{1}{k+1}\sum\limits_{i=1}^{k+1}\dbinom{k+1}{i}B_{k+1-i}n^i\)
于是 \(s_k(n)=\sum\limits_{i=1}^{n}i^k=\dfrac{1}{k+1}\sum\limits_{i=1}^{k+1}\dbinom{k+1}{i}B_{k+1-i}n^i+n^k\)，\(s_k\) 是一个关于 \(n\) 的 \(k+1\) 次多项式。于是我们可以 \(O(N^2)\) 预处理出 \(k\le N\) 的 \(s_k\) 的每项系数。
推式子：
\( \begin{aligned} &\quad\sum\limits_{i=1}^n (i,n)^x((ni)/(i,n))^y=n^y\sum\limits_{i=1}^n (i,n)^{x-y}i^y\\ &=n^y\sum\limits_{d\mid n}d^{x-y}\sum\limits_{i=1}^n i^y[(i,n)=d]=n^y\sum\limits_{d\mid n}d^{x-y}\sum\limits_{i=1}^{n/d} (id)^y[(i,n/d)=1]\\ &=n^y\sum\limits_{d\mid n}d^x\sum\limits_{i=1}^{n/d} i^y\sum\limits_{D\mid (i,n/d)}\mu(D)=n^y\sum\limits_{d\mid n}d^x\sum\limits_{D\mid n/d}\mu(D)D^y\sum\limits_{i=1}^{n/Dd} i^y\\ &=n^y\sum\limits_{D\mid n}\mu(D)D^y\sum\limits_{d\mid n/D}(Dd/D)^xs(n/dD)=n^y\sum\limits_{D\mid n}\mu(D)D^{y-x}\sum\limits_{Dd\mid n}(Dd)^xs(n/Dd) \end{aligned} \)
其中 \(s(x)=\sum\limits_{i=1}^x i^y\)，就是自然数幂和。接下来我们套路性的令 \(T=Dd\)：
\(n^y\sum\limits_{D\mid n}\mu(D)D^{y-x}\sum\limits_{Dd\mid n}(Dd)^xs(n/Dd)=n^y\sum\limits_{T\mid n}T^xs(n/T)\sum\limits_{D\mid T} \mu(D)D^{y-x}\)
我们可以把 \(s\) 用前置知识拆成多项式，令 \(s(n)=\sum\limits_{i=1}^{y+1}a_in^i\)，则
\(n^y\sum\limits_{T\mid n}T^xs(n/T)\sum\limits_{D\mid T} \mu(D)D^{y-x}=n^y\sum\limits_{i=1}^{y+1}a_i\sum\limits_{T\mid n}T^x(n/T)^{i}\sum\limits_{D\mid T} \mu(D)D^{y-x}=n^y\sum\limits_{i=1}^{y+1}a_in^i\sum\limits_{T\mid n}T^{x-i}\sum\limits_{D\mid T} \mu(D)D^{y-x}\)
注意到 \(a_in^i\) 后面很积性啊，有结论：若 \(f\) 为积性函数，\(n\) 质因数分解为：\(n=\prod\limits_{i=1}^kp_i^{\alpha_i}\)，则：
\(g(n)=\sum\limits_{d\mid n}f(d)=\prod\limits_{i=1}^k\left(\sum\limits_{j=0}^{\alpha_i} f(p_i^j)\right)\)，也是个积性。
考虑素数幂处的值（可以用 FFT 中的点值来理解下面）
令 \(f(n)=\sum\limits_{d\mid n} \mu(d)d^{y-x}\) 为积性函数，则对每个 \(n\) 的素数幂因子，\(f(p^\alpha)\) 是可以快速预处理的。
\(\sum\limits_{T\mid n}T^{x-i}\sum\limits_{D\mid T} \mu(D)D^{y-x}=\sum\limits_{T\mid n}T^{x-i}f(T)=\prod\limits_{i=1}^k\left(\sum\limits_{j=0}^{\alpha_i} p_i^{j(x-i)}f(p_i^j)\right)\)
这就很好算啦！先 Pollard-Rho 分解质因子，而后枚举原式中的 \(i\)，对于每个 \(n\) 的素因子计算后面那坨即可，注意 \(f\) 不受原式 \(i\) 的影响，于是只用预处理一次即可，当然这是常数优化。
复杂度：\(O(T(n^{1/4}+yw(n)\log n))\)，其中 \(w(n)\) 为 \(n\) 的质因子个数，为 \(O(\log n)\) 级别。
几乎没怎么优化就跑得很快，目前谷 rk3。代码：
#include<bits/stdc++.h>
#define P pair<int,int>
#define fi first
#define se second
#define LL long long
#define bll __int128
#define fr(x) freopen(#x".in","r",stdin);freopen(#x".out","w",stdout);
using namespace std;
mt19937 rnd(time(0));
const int N=3e3+5,M=N-5,mod=1e9+7;
int T,B[N],jc[N],inv[N],I[N],a[N][N],x,y,pw,ans;LL n;
inline int md(int x){return x>=mod?x-mod:x;}
inline int ksm(int x,int p){int s=1;for(;p;(p&1)&&(s=1ll*s*x%mod),x=1ll*x*x%mod,p>>=1);return s;}
inline int C(int n,int m){return 1ll*jc[n]*inv[m]%mod*inv[n-m]%mod;}
int t,c[95],f[95][95],g[95][95];LL PP[95];
namespace PRHO
{
#define mytz __builtin_ctzll
#define Abs(x) ((x)>0?(x):-(x))
inline LL gcd(LL a,LL b)
{
LL az=mytz(a),bz=mytz(b),z=min(az,bz),diff;b>>=bz;
while(a) a>>=az,diff=a-b,az=mytz(diff),b=min(a,b),a=Abs(diff);return b<<z;
}
inline LL ksm(LL x,LL p,LL mod){LL s=1;for(;p;(p&1)&&(s=(bll)s*x%mod),x=(bll)x*x%mod,p>>=1);return s;}
const LL pr[]={2,3,5,7,11,13,17,19,23,29,31,37};
inline bool check(LL a,LL p)
{
LL d=a-1,t=0;while(~d&1) d>>=1,t++;LL now=ksm(p,d,a);
if(now==1||now==0||now==a-1) return 1;
for(int i=0;i<t;i++)
{
now=(bll)now*now%a;
if(now==1) return 0;
if(now==a-1&&i!=t-1) return 1;
}
return 0;
}
inline bool pd(LL x)
{
if(x==1) return 0;
for(LL i:pr)
{
if(x==i) return 1;
if(x%i==0||!check(x,i)) return 0;
}return 1;
}
#define f(x,c,n) (((bll)(x)*(x)+(c))%(n))
inline LL Find(LL x)
{
LL t1=1ll*rnd()*rnd()%(x-1)+1,c=1ll*rnd()*rnd()%(x-1)+1,t2=f(t1,c,x),d,mul=1;
for(int i=1;;i<<=1,t1=t2,mul=1)
{
for(int j=1;j<=i;j++)
{
t2=f(t2,c,x);
mul=(bll)mul*Abs(t1-t2)%x;
if(j%127==0){d=gcd(mul,x);if(d>1) return d;}
}d=gcd(mul,x);
if(d>1) return d;
}
}
void po(LL x)
{
if(x==1) return;
if(pd(x)) return PP[++t]=x,void();LL num=Find(x);
while(x%num==0) x/=num;po(x),po(num);
}
inline void bk(LL x)
{
t=c[1]=0;po(x);sort(PP+1,PP+1+t);t=unique(PP+1,PP+1+t)-PP-1;
for(int i=1;i<=t;c[++i]=0) while(x%PP[i]==0) c[i]++,x/=PP[i];
}
}using PRHO::bk;
inline void sol1()
{
for(int i=1;i<=t;i++) PP[i]%=mod;
for(int i=1;i<=t;i++)
{
int w=ksm(PP[i],pw);
for(int j=f[i][0]=1,s=w;j<=c[i];j++,s=1ll*s*w%mod)
f[i][j]=f[i][j-1],(j<=1)&&(f[i][j]=md(f[i][j]-s+mod));
}
}
inline int sol2(int i)
{
int PW=(x-i+mod-1)%(mod-1),ans=1;
for(int i=1,w=ksm(PP[1],PW);i<=t;w=ksm(PP[++i],PW))
for(int j=0,s=1;j<=c[i];j++,s=1ll*s*w%mod) g[i][j]=1ll*f[i][j]*s%mod;
for(int i=1,s=0;i<=t;i++,ans=1ll*ans*s%mod,s=0) for(int j=0;j<=c[i];j++) s=md(s+g[i][j]);
return ans;
}
int main()
{
ios::sync_with_stdio(0);cin.tie(0);cout.tie(0);cin>>T;
for(int i=jc[0]=1;i<=M+1;i++) jc[i]=1ll*jc[i-1]*i%mod;
inv[M+1]=ksm(jc[M+1],mod-2);for(int i=M;i>=0;i--) inv[i]=1ll*inv[i+1]*(i+1)%mod;
I[1]=1;for(int i=2;i<=M+1;i++) I[i]=mod-1ll*(mod/i)*I[mod%i]%mod;
for(int i=B[0]=1;i<=M;i++)
{
for(int j=0;j<i;j++) B[i]=(B[i]+1ll*B[j]*C(i+1,j))%mod;
B[i]=md(mod-1ll*I[i+1]*B[i]%mod);
}
for(int i=0;i<=M;a[i][i]++,i++) for(int j=0;j<=i;j++) a[i][i+1-j]=1ll*B[j]*C(i+1,j)%mod*I[i+1]%mod;
while(T--)
{
cin>>n>>x>>y;bk(n);pw=(y-x+mod-1)%(mod-1);sol1();ans=0;
for(int i=1;i<=y+1;i++) ans=(ans+1ll*a[y][i]*ksm(n%mod,i)%mod*sol2(i))%mod;
cout<<(1ll*ans*ksm(n%mod,y)%mod)<<"\n";
}
return 0;
}
