Powerful number 筛略解

Powerful number 筛可以用来求出一类积性函数的前缀和,最快可以达到根号复杂度。


  • 以下字母大都代表整数;
  • 以下 \(p\) 代表质数;
  • 以下 \(*\) 代表狄利克雷卷积;
  • 以下 \(a/b=\lfloor{a\over b}\rfloor\)
  • 以下任何不够严谨的地方都请忽略

定义

定义一个 Powerful number 为任意一个质因子的次数都不小于 2 的数,如 \(72=2^3\times 3^2\)\(4000=2^5\times5^3\)。Powerful number 有如下性质:

  • 一个 Powerful number 一定可以表示为 \(a^2b^3\) 的形式,如 \(72=3^2\times2^3\)\(4000=2^{(2+3)}\times5^3=2^2\times(2\times5)^3=2^2\times10^3\)
  • \(n\) 以内的 Powerful number 个数是 \(O(\sqrt n)\) 的:

    \[\int_1^{\sqrt n}\sqrt[3]{n\over x^2}dx=O(\sqrt n) \]

算法

假设有一个积性函数 \(f\),我们希望快速求出 \(f\) 的前缀和。我们尝试找到另一个积性函数 \(g\) 使得 \(g(p)=f(p)\),且可以较快求出 \(g\) 的前缀和。

再找到一个积性函数 \(h\),使得 \(f=g*h\)

显然 \(f(p)=g(1)h(p)+g(p)h(1)=h(p)+g(p)=h(p)+f(p)\),所以 \(h(p)=0\)。又因为 \(h\) 为积性函数,所以所有非 Powerful number 的数的 \(h\) 值都为 0。

\[\begin{aligned}\sum\limits_{i=1}^n f(i)&=\sum\limits_{i=1}^n g*h(i)\\ &=\sum\limits_{i=1}^n\sum\limits_{d|i}g(d)h(i/d)\\ &=\sum\limits_{i=1}^n h(i)\sum\limits_{j=1}^{n/i}g(j)\end{aligned} \]

因此枚举 \(\sqrt n\) 个 Powerful number 的 \(h\) 值并求出对应的 \(g\) 的前缀和便能求出 \(f\) 的前缀和。

难点在于如何构造 \(g,h\),下面举两个例子。

例 1

积性函数 \(f(p^q)=p^k\)\(k\) 为常数)。

构造 \(g(n)=n^k\),显然 \(f(p)=g(p)=p^k\)

考虑找到 \(g^{-1}\) 使 \(g^{-1}*g=\epsilon\),则 \(h=f*g^{-1}\)

\[\begin{aligned}g*g^{-1}(p^q)=&[p^q=1]=[q=0]\\ =&\sum\limits_{i=0}^qg(p^i)g^{-1}(p^{q-i})\\ =&\sum\limits_{i=0}^q (p^k)^ig^{-1}(p^{q-i})\end{aligned} \]

假如 \(q>0\),还可以推到:

\[\begin{aligned}g*g^{-1}(p^q)&=0\\&=g*g^{-1}(p^{q-1})\cdot p^k+g^{-1}(p^q)\end{aligned} \]

\(q=0\) 开始手玩,很容易发现 \(g^{-1}(p^0)=1\)\(g^{-1}(p^1)=-p^k\)\(g^{-1}(p^q)=0(q>1)\)。手玩一下 \(f*g^{-1}\) 可得 \(h=p^k-p^{2k}\)

例 2

积性函数 \(f(p^q)=p\oplus q\)。(LOJ6053 简单的函数

显然 \(f(p)=\begin{cases}p-1\quad&(p\not=2)\\p+1&(p=2)\end{cases}\)

构造 \(g(n)=\begin{cases}\varphi(n)\quad&(2\nmid n)\\3\varphi(n)\quad&(2\mid n)\end{cases}\)

下面主要解决两个问题:\(g^{-1}\)\(g\) 的前缀和。

用类似于例 1 的方法推导一番可得:\(q\) 足够大时,\(g*g^{-1}(p^q)=\begin{cases} g^{-1}(p^q)-g^{-1}(p^{q-1})+p\cdot g*g^{-1}(p^{q-1})\quad(p\not =2)\\ g^{-1}(p^q)+g^{-1}(p^{q-1})+p\cdot g*g^{-1}(p^{q-1})\quad(p=2)\end{cases}\)
\(g^{-1}(p^q)=\begin{cases} g^{-1}(p^{q-1})\quad&(p\not =2)\\ -g^{-1}(p^{q-1})\quad&(p=2)\end{cases}\)

所以:

  • \(g^{-1}(1)=1\)
  • \(g^{-1}(p^q)=\begin{cases}-p+1\quad&(p\not =2)\\(-1)^q(p+1)\quad&(p=2)\end{cases}\quad(q>0)\)

然后推导 \(g\) 的前缀和:\(\sum\limits_{i=1}^{n}g(i)=\sum\limits_{i=1}^n\varphi(i)+2\sum\limits_{i=1}^{n/2}\varphi(2i)\)

\(S'(n)=\sum\limits_{i=1}^{n}\varphi(2i)\)。则当 \(2\mid n\) 时:

\[\begin{aligned}S'(n)=&\sum\limits_{i=1}^n\varphi(2i)\\ =&\sum\limits_{i=1}^{n/2}(\varphi(2\cdot (2i-1))+\varphi(2\cdot 2i))\\ =&\sum\limits_{i=1}^{n/2}(\varphi(2i-1)+2\varphi(2i))\\ =&\sum\limits_{i=1}^{n/2}\varphi(2i)+\sum\limits_{i=1}^{n/2}(\varphi(2i-1)+\varphi(2i))\\ =&S'(n/2)+\sum\limits_{i=1}^{n}\varphi(i)\end{aligned}\]

假如 \(2\nmid n\)\(S'(n)=S'((n-1)/2)+\sum\limits_{i=1}^{n-1}\varphi(i)+\varphi(2n)=S'(n/2)+\sum\limits_{i=1}^n\varphi(i)\)(完 全 一 致)。

因此在杜教筛出 \(O(\sqrt n)\) 个点的 \(\varphi\) 的前缀和后,再处理出 \(O(\sqrt n)\)\(S'\) 的值。\(\sum\limits_{i=1}^{n}g(i)=\sum\limits_{i=1}^n\varphi(i)+2S'(n/2)\)

于是我们得到了 \(g\) 的前缀和,并且可以通过暴力求出 \(f*g^{-1}(p^q)\) 得到 \(h(p^q)\) 的值,进一步在搜索时同步求出 \(h(x)\)

代码:

/**********
Author: WLBKR5
Problem: loj 6053
Name: 简单的函数 
Source: uploaded by fjzzq2002
Algorithm: powerful number 
Date: 2020/06/02
Statue: accepted
Submission: loj.ac/submission/823642
**********/
#include<bits/stdc++.h>
using namespace std;
#define ll long long
ll getint(){
	ll ans=0,f=1;
	char c=getchar();
	while(c<'0'||c>'9'){
		if(c=='-')f=-1;
		c=getchar();
	}
	while(c>='0'&&c<='9'){
		ans=ans*10-'0'+c;
		c=getchar();
	}
	return ans*f;
}
const int mod=1e9+7,N=1e7+10;
ll n;int sq;
bool boo[N+10];
int pri[1000100],phi[N+10],cnt=0;
void sieve(int n){
	//线性筛 
	boo[1]=1;
	phi[1]=1;
	for(int i=2;i<=n;i++){
		if(!boo[i]){
			pri[cnt++]=i;
			phi[i]=i-1;
		}
		for(int j=0;j<cnt&&pri[j]*i<=n;++j){
			boo[pri[j]*i]=1;
			if(i%pri[j]==0){
				phi[pri[j]*i]=phi[i]*pri[j];
				break;
			}else{
				phi[pri[j]*i]=phi[i]*(pri[j]-1);
			}
		}
	}
	for(int i=1;i<=n;i++){
		phi[i]=(phi[i]+phi[i-1])%mod;
	}
}
int calc(ll x){
	//x*(x+1)/2
	int ans=x%mod*((x+1ll)%mod)%mod;//此处处理不好可能炸 long long 
	if(ans&1)ans+=mod;
	return ans>>1;
}
ll s[10010];//sum_1^{n/i} g(i)
int s_1[200010],s_2[200010];
ll ans=0;
int h[30100][66],d[30100];
int& s_(ll x){
	//S'
	return x<=200000?s_1[x]:s_2[n/x];
}
ll sumphi(ll x){
	//phi 的前缀和 
	return x<=N?phi[x]:s[n/x];
}
ll sumg(ll x){
	//g 的前缀和 
	return (sumphi(x)+2*s_(x/2))%mod;
}
void dfs(ll x,ll v,ll p){
	//枚举 x 
	//v=h(x) 
	//h(x)*\sum_1^{n/x}g(j)
	ans=(ans+v*1ll*sumg(n/x)%mod)%mod;
	for(int i=p;i<cnt;i++){
		if(x*1ll*pri[i]*pri[i]>n)break;
		ll y=x*pri[i];
		for(int j=1;y<=n;++j,y*=pri[i]){
			if(j>d[i]){
				++d[i];
				if(pri[i]==2){
					h[i][j]=((pri[i]^j)+(j%2?mod-pri[i]-1ll:pri[i]+1ll))%mod;
					for(int k=1;k<j;k++){
						h[i][j]=(h[i][j]+
							(pri[i]^k)*((j-k)%2?mod-pri[i]-1ll:pri[i]+1ll)%mod)%mod;
					}
				}else{
					h[i][j]=((pri[i]^j)+1ll-pri[i]+mod)%mod;
					for(int k=1;k<j;k++){
						h[i][j]=(h[i][j]+(pri[i]^k)*(mod-pri[i]+1ll)%mod)%mod;
					}
				}
			}
			if(!h[i][j])continue;
			dfs(y,v*1ll*h[i][j]%mod,i+1);
		}
	}
}

signed main(){
	n=getint();
	sieve(N);//预处理质数、phi 的前缀和 
	for(int i=1;i<=cnt&&pri[i]*1ll*pri[i]<=n;i++){
		h[i][0]=1;
	}
	int u=1;while(n/u>=N)++u;
	for(int i=u;i>=1;--i){
		//杜教筛筛出 phi 的前缀和 
		//phi*1=Id
		ll n_=n/i;
		s[i]=calc(n_);
		for(ll l=2,r=0;l<=n_;l=r+1){
			r=n_/(n_/l);
			s[i]=(s[i]-(r-l+1ll)%mod*1ll*sumphi(n_/l)%mod+mod)%mod;
		}
	}
	vector<ll>v;
	for(ll l=1,r=0;l<=n;l=r+1){
		r=n/(n/l);
		v.push_back(n/l);
	}
	for(int i=v.size()-1;i>=0;i--){
		//算出 S' 的前缀和 
		s_(v[i])=(s_(v[i]/2)+sumphi(v[i]))%mod;
	}
	dfs(1,1,0);
	cout<<(ans%mod+mod)%mod;
	return 0;
}

(原发布于 csdn

知识共享许可协议
若文章内无特别说明,公开文章采用知识共享署名-相同方式共享 4.0 国际许可协议进行许可。
posted @ 2020-10-30 12:33  破壁人五号  阅读(249)  评论(0编辑  收藏  举报