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

题面

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

题解

前置知识

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

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

\[=n^y{\sum\limits_{i=1}^{n}}i^y(i,n)^{x-y} \]

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

其中，

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

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

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

利用伯努利数，有

\[{\sum\limits_{i=1}^{\frac{n}{d}}}i^y={\frac{1}{y+1}}{\sum\limits_{i=0}^{y}}{\tbinom{y+1}{i}}B_i(\frac{n}{d})^{y+1-i} \]

所以原式=

\[n^y*{\frac{1}{y+1}}*{\sum\limits_{i=0}^{y}}{\tbinom{y+1}{i}}B_i{\sum\limits_{d|n}}g(d)d^y(\frac{n}{d})^{y+1-i} \]

\[=n^y*{\frac{1}{y+1}}*{\sum\limits_{i=0}^{y}}{\tbinom{y+1}{i}}B_in^{y+1-i}{\sum\limits_{d|n}}g(d)d^{i-1} \]

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

\[B_i=\begin{cases} 1(i=0) \\ 1-{\sum\limits_{j=0}^{i-1}}B_j{\tbinom{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}}\)，就可以完成下列转化：

\[{\sum\limits_{d|n}}h(d) \]

\[={\prod\limits_{l=1}^{k}}{\sum\limits_{j=0}^{\alpha_l}}h(p_l^{j}) \]

现在我们再拆开h。

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

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

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

\[h(p^j)=\begin{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(T\log 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}}+Ty\log 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;
}
posted @ 2020-02-11 21:13  coder66  阅读(382)  评论(0)  编辑  收藏  举报