超过 90% 的人都不知道!居然能这样用 O(√N) 次运算求出 N 个数之和!网友:长见识了 | Powerful number 筛

笔记 同时被 2 个专栏收录
10 篇文章 0 订阅

标题党

不过标题好像没有什么不符合事实的地方。

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


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

定义

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

  • 一个 Powerful number 一定可以表示为 a 2 b 3 a^2b^3 a2b3 的形式,如 72 = 3 2 × 2 3 72=3^2\times2^3 72=32×23 4000 = 2 ( 2 + 3 ) × 5 3 = 2 2 × ( 2 × 5 ) 3 = 2 2 × 1 0 3 4000=2^{(2+3)}\times5^3=2^2\times(2\times5)^3=2^2\times10^3 4000=2(2+3)×53=22×(2×5)3=22×103
  • n n n 以内的 Powerful number 个数是 O ( n ) O(\sqrt n) O(n ) 的: ∫ 1 n n x 2 3 d x = O ( n ) \int_1^{\sqrt n}\sqrt[3]{n\over x^2}dx=O(\sqrt n) 1n 3x2n dx=O(n )

算法

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

再找到一个积性函数 h h h,使得 f = g ∗ h f=g*h f=gh

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

∑ i = 1 n f ( i ) = ∑ i = 1 n g ∗ h ( i ) = ∑ i = 1 n ∑ d ∣ i g ( d ) h ( i / d ) = ∑ i = 1 n h ( i ) ∑ j = 1 n / i g ( j ) \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} i=1nf(i)=i=1ngh(i)=i=1ndig(d)h(i/d)=i=1nh(i)j=1n/ig(j)

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

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

例 1

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

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

考虑找到 g − 1 g^{-1} g1 使 g − 1 ∗ g = ϵ g^{-1}*g=\epsilon g1g=ϵ,则 h = f ∗ g − 1 h=f*g^{-1} h=fg1

g ∗ g − 1 ( p q ) = [ p q = 1 ] = [ q = 0 ] = ∑ i = 0 q g ( p i ) g − 1 ( p q − i ) = ∑ i = 0 q ( p k ) i g − 1 ( p q − i ) \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} gg1(pq)===[pq=1]=[q=0]i=0qg(pi)g1(pqi)i=0q(pk)ig1(pqi)

假如 q > 0 q>0 q>0,还可以推到:
g ∗ g − 1 ( p q ) = 0 = g ∗ g − 1 ( p q − 1 ) ⋅ p k + g − 1 ( p q ) \begin{aligned}g*g^{-1}(p^q)&=0\\&=g*g^{-1}(p^{q-1})\cdot p^k+g^{-1}(p^q)\end{aligned} gg1(pq)=0=gg1(pq1)pk+g1(pq)

q = 0 q=0 q=0 开始手玩,很容易发现 g − 1 ( p 0 ) = 1 g^{-1}(p^0)=1 g1(p0)=1 g − 1 ( p 1 ) = − p k g^{-1}(p^1)=-p^k g1(p1)=pk g − 1 ( p q ) = 0 ( q > 1 ) g^{-1}(p^q)=0(q>1) g1(pq)=0(q>1)。手玩一下 f ∗ g − 1 f*g^{-1} fg1 可得 h = p k − p 2 k h=p^k-p^{2k} h=pkp2k

例 2

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

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

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

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

用类似于例 1 的方法推导一番可得: q q q 足够大时, g ∗ g − 1 ( p q ) = { g − 1 ( p q ) − g − 1 ( p q − 1 ) + p ⋅ g ∗ g − 1 ( p q − 1 ) ( p ≠ 2 ) g − 1 ( p q ) + g − 1 ( p q − 1 ) + p ⋅ g ∗ g − 1 ( p q − 1 ) ( p = 2 ) 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} gg1(pq)={g1(pq)g1(pq1)+pgg1(pq1)(p=2)g1(pq)+g1(pq1)+pgg1(pq1)(p=2)
g − 1 ( p q ) = { g − 1 ( p q − 1 ) ( p ≠ 2 ) − g − 1 ( p q − 1 ) ( p = 2 ) 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} g1(pq)={g1(pq1)g1(pq1)(p=2)(p=2)

所以:

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

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

S ′ ( n ) = ∑ i = 1 n φ ( 2 i ) S'(n)=\sum\limits_{i=1}^{n}\varphi(2i) S(n)=i=1nφ(2i)。则当 2 ∣ n 2\mid n 2n 时: S ′ ( n ) = ∑ i = 1 n φ ( 2 i ) = ∑ i = 1 n / 2 ( φ ( 2 ⋅ ( 2 i − 1 ) ) + φ ( 2 ⋅ 2 i ) ) = ∑ i = 1 n / 2 ( φ ( 2 i − 1 ) + 2 φ ( 2 i ) ) = ∑ i = 1 n / 2 φ ( 2 i ) + ∑ i = 1 n / 2 ( φ ( 2 i − 1 ) + φ ( 2 i ) ) = S ′ ( n / 2 ) + ∑ i = 1 n φ ( i ) \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} S(n)=====i=1nφ(2i)i=1n/2(φ(2(2i1))+φ(22i))i=1n/2(φ(2i1)+2φ(2i))i=1n/2φ(2i)+i=1n/2(φ(2i1)+φ(2i))S(n/2)+i=1nφ(i)

假如 2 ∤ n 2\nmid n 2n S ′ ( n ) = S ′ ( ( n − 1 ) / 2 ) + ∑ i = 1 n − 1 φ ( i ) + φ ( 2 n ) = S ′ ( n / 2 ) + ∑ i = 1 n φ ( i ) 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) S(n)=S((n1)/2)+i=1n1φ(i)+φ(2n)=S(n/2)+i=1nφ(i)(完 全 一 致)。

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

于是我们得到了 g g g 的前缀和,并且可以通过暴力求出 f ∗ g − 1 ( p q ) f*g^{-1}(p^q) fg1(pq) 得到 h ( p q ) h(p^q) h(pq) 的值,进一步在搜索时同步求出 h ( x ) h(x) 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;
}

  • 0
    点赞
  • 0
    评论
  • 1
    收藏
  • 一键三连
    一键三连
  • 扫一扫,分享海报

©️2021 CSDN 皮肤主题: 深蓝海洋 设计师:CSDN官方博客 返回首页
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、C币套餐、付费专栏及课程。

余额充值