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;
}