[BZOJ2627][洛谷P4464]JZPKIL(莫比乌斯反演+Pollard-rho+伯努利数)

题面

https://www.luogu.com.cn/problem/P4464

题解

前置知识

• Pollard-rho https://www.luogu.com.cn/problemnew/solution/P4718

• Miller-rabin https://www.cnblogs.com/xh092113/p/12288424.html

• 莫比乌斯反演：《具体数学：计算机科学基础》4.9节

  以及它的基础应用： https://www.cnblogs.com/xh092113/p/12268356.html

• 伯努利数：《具体数学：计算机科学基础》6.5节

本题要求({sum_{i=1}^{n}}(i,n)^x[i,n]^y)。上手即反演：

[={sumlimits_{i=1}^{n}}(i,n)^{x-y}i^yn^y ]

[=n^y{sumlimits_{i=1}^{n}}i^y(i,n)^{x-y} ]

[=n^y{sumlimits_{i=1}^{n}}i^y{sumlimits_{d|(i,n)}}g(d) ]

其中，

[g(d)={sumlimits_{p|d}}{mu(p)}({frac{d}{p}})^{x-y} ]

[原式=n^y{sumlimits_{d|n}}g(d){sumlimits_{i=1}^{n}}i^y[d|i] ]

[=n^y{sumlimits_{d|n}}g(d)d^y{sumlimits_{i=1}^{frac{n}{d}}}i^y ]

利用伯努利数，有

[{sumlimits_{i=1}^{frac{n}{d}}}i^y={frac{1}{y+1}}{sumlimits_{i=0}^{y}}{ binom{y+1}{i}}B_i(frac{n}{d})^{y+1-i} ]

所以原式=

[n^y*{frac{1}{y+1}}*{sumlimits_{i=0}^{y}}{ binom{y+1}{i}}B_i{sumlimits_{d|n}}g(d)d^y(frac{n}{d})^{y+1-i} ]

[=n^y*{frac{1}{y+1}}*{sumlimits_{i=0}^{y}}{ binom{y+1}{i}}B_in^{y+1-i}{sumlimits_{d|n}}g(d)d^{i-1} ]

前面两项均为常数；由于伯努利数满足

[B_i=egin{cases} 1(i=0) \ 1-{sumlimits_{j=0}^{i-1}}B_j{ binom{i+1}{j}}{frac{1}{i+1}} end{cases} ]

的递推式，所以所有的(B_i)可以(O(y^2))预处理出来，也变成了常数。所以问题的瓶颈就落在了计算({sum_{d|n}}g(d)d^{i-1})上。由于前面已经循环过i=0到y,并且还有T组数据，所以这里必须在最多(O(log n))内解决。不过，我们能够看出(g(d)d^{i-1})是一个积性函数。所以，设它是(h(d))，并设n的素因数分解式是({prod_{i=1}^{k}}p_k^{alpha_{k}})，就可以完成下列转化：

[{sumlimits_{d|n}}h(d) ]

[={prodlimits_{l=1}^{k}}{sumlimits_{j=0}^{alpha_l}}h(p_l^{j}) ]

现在我们再拆开h。

[h(d)=g(d)d^{i-1} ]

[={sumlimits_{e|d}}{mu(e)}(frac{d}{e})^{x-y}d^{i-1} ]

对于质数p和自然数j，发现

[h(p^j)=egin{cases} 1(j=0) \ p^{j(x-y+i-1)}-p^{j(x-y+i-1)+y-x} end{cases}(j{geq}1) ]

而且，(x-y+i-1)和(y-x)都处于([-3001,6000])的范围之中，这意味着我们可以对于所有的n的质因子p，预处理出p的-3001到6000次方，这样的预处理总时间是(O(Tlog n*y))，可以接受。在j>1的时候，只需要乘上(p^{x-y+i-1})，就可以从(h(p^{j-1}))转移过来了。

于是，计算({sum_{d|n}}h(d))的时间复杂度降为了(O({sum}{alpha_i}))即(O(log n))。

过程中，需要对n做素因数分拆；由于n过大，所以需要使用Pollard-rho配合Miller-rabin实现，其时间复杂度约为(O(n^{frac{1}{4}}))。

总时间复杂度(O(Tn^{frac{1}{4}}+Tylog n+y^2))。

代码

#include<bits/stdc++.h>

using namespace std;

#define ll long long
#define ld long double
#define rg register
#define M 3000
#define Mod 1000000007

namespace ModCalc{
	inline void Inc(ll &x,ll y,ll mod){
		x += y;if(x >= mod)x -= mod;
	}
	
	inline void Dec(ll &x,ll y,ll mod){
		x -= y;if(x < 0)x += mod;
	}
	
	inline void Adjust(ll &x,ll mod){
		x = (x % mod + mod) % mod;
	}
	
	inline void Tms(ll &x,ll y,ll mod){ //O(1)快速乘 
		x = x * y - (ll)((ld)x * y / mod + 1e-3) * mod;
		Adjust(x,mod);
	}
	
	inline ll Add(ll x,ll y,ll mod){
		Inc(x,y,mod);return x;
	}
	
	inline ll Sub(ll x,ll y,ll mod){
		Dec(x,y,mod);return x;
	}
	
	inline ll Check(ll x,ll mod){
		Adjust(x,mod);return x;
	}
	
	inline ll Mul(ll x,ll y,ll mod){
		Tms(x,y,mod);return x;
	}
}
using namespace ModCalc;

inline ll read(){
	ll s = 0,ww = 1;
	char ch = getchar();
	while(ch < '0' || ch > '9'){if(ch == '-')ww = -1;ch = getchar();}
	while('0' <= ch && ch <= '9'){s = 10 * s + ch - '0';ch = getchar();}
	return s * ww;
}

inline ll power(ll a,ll n,ll mod){
	ll x = a,s = 1;
	while(n){
		if(n & 1)Tms(s,x,mod);
		Tms(x,x,mod);
		n >>= 1;
	}	
	return s;
}

ll B[M+5],b[M+5];
ll jc[M+5],iv[M+5]; //阶乘和阶乘的逆元 

inline void prepro(){ //预处理jc,iv 
	jc[0] = 1;
	for(rg ll i = 1;i <= M + 1;i++)jc[i] = jc[i-1] * i % Mod;
	iv[M+1] = power(jc[M+1],Mod - 2,Mod);
	for(rg ll i = M;i >= 0;i--)iv[i] = iv[i+1] * (i + 1) % Mod;
}

inline void calcB(){
	B[0] = 1;
	for(rg ll i = 1;i <= M;i++){
		ll cur = 0;
		for(rg ll j = 0;j < i;j++)Inc(cur,jc[i+1] * iv[j] % Mod * iv[i+1-j] % Mod * B[j] % Mod,Mod); 
		cur = cur * jc[i] % Mod * iv[i+1] % Mod;
		B[i] = Sub(1,cur,Mod);
	}
}

inline ll test(ll p,ll n){
	if(n == p)return 1;
	ll m = n - 1;
	while(!(m & 1))m >>= 1;
	ll x = power(p,m,n);
	if(x == 1)return 1;
	for(;m < n - 1;m <<= 1){
		ll y = Mul(x,x,n);
		if(y == 1)return x == n - 1; 
		x = y;
	}
	return 0;
}

ll pri[9] = {2,3,5,7,11,13,17,19,23};

inline ll MR(ll n){ //Miller-rabin
	if(n < 2)return 0;
	if(!(n & 1))return n == 2;
	for(rg ll i = 0;i < 9;i++)if(!test(pri[i],n))return 0;
	return 1;
}

inline ll f(ll x,ll c,ll mod){
	return Add(Mul(x,x,mod),c,mod);
}

inline ll gcd(ll a,ll b){
	return b ? gcd(b,a % b) : a;
}

inline ll PR(ll n){ //Pollard-rho
	ll c = 1ll * rand() % (n - 1) + 1,x = 0,y = 0,prod = 1;
	for(rg ll step = 1;;step <<= 1){
		for(rg ll i = 1;i <= step;i++){
			y = f(y,c,n);
			Tms(prod,abs(x-y),n);
			if(i % 127 == 0){
				ll d = gcd(prod,n);
				if(d > 1)return d;
			}
		}
		ll d = gcd(prod,n);
		if(d > 1)return d;
		x = y;
	}
}

vector<ll>fac; //PR算法分拆出的无序质因子 
vector<pair<ll,ll> >Pfac; //整理好后的质因子和对应的次数 

inline void divide(ll n){ //分拆n 
	if(n == 1)return;
	if(MR(n)){
		fac.push_back(n);
		return;
	}
	ll d = PR(n);
	while(d == n)d = PR(n);
	divide(d);
	divide(n / d);
}

ll pw[60+5][3*M+5]; //n的所有素因子的-3001~6000次方 

int main(){
	srand((int)time(NULL));
	prepro();
	calcB();
	ll T = read();
	while(T--){
		ll n = read(),x = read(),y = read();
		fac.resize(0);
		Pfac.resize(0);
		divide(n);
		sort(fac.begin(),fac.end());
		ll p = 0,v = 0;
		for(rg ll i = 0;i < fac.size();i++){
			if(fac[i] != p){
				if(p)Pfac.push_back(make_pair(p,v));
				p = fac[i],v = 1;
			}
			else v++;
		}
		if(p)Pfac.push_back(make_pair(p,v));
		for(rg ll i = 0;i < Pfac.size();i++){ //预处理pw 
			ll p = Pfac[i].first,vp = power(p % Mod,Mod - 2,Mod);
			pw[i][3001] = 1;
			for(rg ll j = 1;j <= 6000;j++)pw[i][j+3001] = pw[i][j+3000] * (p % Mod) % Mod;
			for(rg ll j = -1;j >= -3001;j--)pw[i][j+3001] = pw[i][j+3002] * vp % Mod;
		}
		ll ans = 0;
		for(rg ll i = 0;i <= y;i++){
			ll b = B[i] * jc[y] % Mod * iv[i] % Mod * iv[y+1-i] % Mod;
			b = b * power(n % Mod,y + 1 - i,Mod) % Mod;
			for(rg ll l = 0;l < Pfac.size();l++){
				ll s = 1,w = pw[l][x-y+i-1+3001];
				ll cur = Sub(1,pw[l][y-x+3001],Mod);
				for(rg ll j = 1;j <= Pfac[l].second;j++){
					cur = cur * w % Mod;
					Inc(s,cur,Mod);
				}
				b = b * s % Mod;
			}
			Inc(ans,b,Mod);
		}
		ans = ans * power(n % Mod,y,Mod) % Mod;
		cout << ans << endl;
	}
	return 0;
}