杜教筛学习笔记 杜教筛学习笔记

接着 Dirichlet 卷积,继续学习杜教筛。

本文中一些函数的定义见此博文

应用

通过杜教筛,我们可以快速求出某一数论函数 (f) 的前缀和,即,我们可以在低于线性的时间复杂度内求出 (S(n)=sumlimits_{i=1}^n f(i))

方法

杜教筛主要运用一个公式,通过这个公式我们建立了 (S(n))(S(frac{n}{i})) 的关系式,已知两个数论函数 (f,g)(S(n)=sumlimits_{i=1}^n f(i)),则有:

[sumlimits_{i=1}^n (fast g)(i)=sumlimits_{i=1}^n g(i)cdot S(lfloorfrac{n}{i} floor) ]

我们给个证明:

[egin{aligned} sumlimits_{i=1}^n (fast g)(i) \=sumlimits_{i=1}^nsumlimits_{dmid i} f(d)cdot g(frac{i}{d})\=sumlimits_{d=1}^nsumlimits_{i=1}^{lfloor frac{n}{d} floor} g(d)cdot f(i)\=sumlimits_{d=1}^n g(d)cdotsumlimits_{i=1}^{lfloor frac{n}{d} floor} f(i)\=sumlimits_{d=1}^n g(d)cdot S(lfloor frac{n}{d} floor)\=sumlimits_{i=1}^n g(i)cdot S(lfloor frac{n}{i} floor)\&square end{aligned}]

证明的关键就在于中间对于 (d,i) 枚举顺序的变化,这个证明值得细品。

那么我们有了这个公式有什么用呢?我们再来导一下:

[egin{aligned} sumlimits_{i=1}^n (fast g)(i)=sumlimits_{i=1}^n g(i)cdot S(lfloorfrac{n}{i} floor)\sumlimits_{i=1}^n (fast g)(i)=g(1)cdot S(n)+sumlimits_{i=2}^n g(i)cdot S(lfloorfrac{n}{i} floor)\g(1)cdot S(n)=sumlimits_{i=1}^n (fast g)(i)-sumlimits_{i=2}^n g(i)cdot S(lfloorfrac{n}{i} floor) end{aligned}]

我们可以轻松地求出 (S(n)) 的值。而 (S(lfloorfrac{n}{i} floor)) 可以通过数论分块来解决。由此看来,杜教筛的时间复杂度为低于线性复杂度。

例子

看了理论,可能还有点晕,我们来看几个例子。

求莫比乌斯函数前缀和

(S(n)=sumlimits_{i=1}^n mu(i) (n<2^{31}))

我们根据公式来求,我们令 (f=mu,g=1),由 Dirichlet 卷积知识,我们知道 (muast 1=varepsilon)。那么:

[egin{aligned} g(1)cdot S(n)=sumlimits_{i=1}^n (fast g)(i)-sumlimits_{i=2}^n g(i)cdot S(lfloorfrac{n}{i} floor)\1cdot S(n)=sumlimits_{i=1}^n varepsilon(i)-sumlimits_{i=2}^n 1cdot S(lfloorfrac{n}{i} floor)\S(n)=1-sumlimits_{i=2}^n S(lfloorfrac{n}{i} floor) end{aligned}]

根据数论分块的知识,(lfloorfrac{n}{i} floor) 的取值有 (O(sqrt{n})) 个。那么时间复杂度为 (O(sumlimits_{i=1}2^{frac{1}{2^i}}) approx O(n^{frac{3}{4}}))

这明显会超时,我们需要进一步优化。其实可以这样想,一部分前缀和我们先预处理出来,另一部分再通过杜教筛来算。我们设预处理前 (k) 个,后 (n-k) 个用杜教筛做,那么此时时间复杂度是 (O(k+(n-k)^{frac{3}{4}})),由基本不等式得,(k=(n-k)^{frac{3}{4}}) 原式取最小值。我们算一下发现 (k=n^{frac{2}{3}}) 是不错的选择,此时复杂度大约为 (O(n^{frac{2}{3}}))

//Don't act like a loser.
//You can only use the code for studying or finding mistakes
//Or,you'll be punished by Sakyamuni!!!
#include<bits/stdc++.h>
#define int long long
using namespace std;

int read() {
	char ch=getchar();
	int f=1,x=0;
	while(ch<'0'||ch>'9') {
		if(ch=='-')
			f=-1;
		ch=getchar();
	}
	while(ch>='0'&&ch<='9') {
		x=x*10+ch-'0';
		ch=getchar();
	}
	return f*x;
}

const int maxn=4e7+10;

int n,mu[maxn],p[maxn/10],cnt,s[maxn];
bool is[maxn];
map<long,long> mp;//记忆化搜索

void sieve_mu(int n) {
	fill(is,is+n+1,1);
	mu[1]=1;
	
	for(int i=2;i<=n;i++) {
		if(is[i]) {
			p[++cnt]=i;
			mu[i]=-1;
		}
		for(int j=1;j<=cnt;j++) {
			if(p[j]*i>n) {
				break;
			}
			
			is[p[j]*i]=0;
			if(i%p[j]==0) {
				mu[i*p[j]]=0;
				break;
			}
			mu[i*p[j]]=-mu[i];
		}
	}
	for(int i=1;i<=n;i++){
		s[i]=s[i-1]+mu[i];
	}
}

int pre(int n) {
	if(n<=4e7) {
		return s[n];
	}
	if(mp[n]) {
		return mp[n];
	}
	int j=2,ret=1;
	for(int i=2;i<=n;i=j+1) {
		j=(n/(n/i));//数论分块
		
		ret-=(j-i+1)*pre(n/i);
	}
	return mp[n]=ret;
}

signed main() {
	sieve_mu(4e7);
	int t=read();
	
	while(t--) {
		printf("%lld
",pre(read()));
	}
	return 0;
}

求欧拉函数前缀和

(S(n)=sumlimits_{i=1}^n varphi(i) (n<2^{31}))

我们根据公式来求,我们令 (f=varphi,g=1),由 Dirichlet 卷积知识,我们知道 (varphiast 1=id)。那么:

[egin{aligned} g(1)cdot S(n)=sumlimits_{i=1}^n (fast g)(i)-sumlimits_{i=2}^n g(i)cdot S(lfloorfrac{n}{i} floor)\1cdot S(n)=sumlimits_{i=1}^n id(i)-sumlimits_{i=2}^n 1cdot S(lfloorfrac{n}{i} floor)\S(n)=frac{n^2+n}{2}-sumlimits_{i=2}^n S(lfloorfrac{n}{i} floor) end{aligned}]

和之前一样,我们还是采取一部分直接算,一部分杜教筛的想法,时间复杂度为 (O(n^{frac{2}{3}}))

//Don't act like a loser.
//You can only use the code for studying or finding mistakes
//Or,you'll be punished by Sakyamuni!!!>
#define int long long
using namespace std;

int read() {
	char ch=getchar();
	int f=1,x=0;
	while(ch<'0'||ch>'9') {
		if(ch=='-')
			f=-1;
		ch=getchar();
	}
	while(ch>='0'&&ch<='9') {
		x=x*10+ch-'0';
		ch=getchar();
	}
	return f*x;
}

const int maxn=1e7+10;

int n,phi[maxn],p[maxn/10],cnt,s[maxn];
bool is[maxn];
map<long,long> mp;

void sieve_euler(int n) {
	fill(is,is+n+1,1);
	phi[1]=1;
	
	for(int i=2;i<=n;i++) {
		if(is[i]) {
			p[++cnt]=i;
			phi[i]=i-1;
		}
		for(int j=1;j<=cnt;j++) {
			if(p[j]*i>n) {
				break;
			}
			
			is[p[j]*i]=0;
			if(i%p[j]==0) {
				phi[i*p[j]]=p[j]*phi[i];
				break;
			}
			phi[i*p[j]]=(p[j]-1)*phi[i];
		}
	}
	for(int i=1;i<=n;i++){
		s[i]=s[i-1]+phi[i];
	}
}

int pre(int n) {
	if(n<=1e7) {
		return s[n];
	}
	if(mp[n]) {
		return mp[n];
	}
	int j=2,ret=(n*(n+1)/2);
	for(int i=2;i<=n;i=j+1) {
		j=(n/(n/i));
		
		ret-=(j-i+1)*pre(n/i);
	}
	return mp[n]=ret;
}

signed main() {
	sieve_euler(1e7);
	int t=read();
	
	while(t--) {
		printf("%lld
",pre(read()));
	}
	return 0;
}

这两个代码结合起来就可以通过洛谷上的模板题。

求莫比乌斯函数平方前缀和

(S(n)=sumlimits_{i=1}^n mu^2(i) (n<2^{31}))

我们继续根据公式来求,我们令 (f=mu^2,g=1),由 Dirichlet 卷积知识,我们知道 (mu^2*1=mu)。那么:

[egin{aligned} g(1)cdot S(n)=sumlimits_{i=1}^n (fast g)(i)-sumlimits_{i=2}^n g(i)cdot S(lfloorfrac{n}{i} floor)\1cdot S(n)=sumlimits_{i=1}^n mu(i)-sumlimits_{i=2}^n 1cdot S(lfloorfrac{n}{i} floor)\S(n)=sumlimits_{i=1}^n mu(i)-sumlimits_{i=2}^n S(lfloorfrac{n}{i} floor) end{aligned}]

所以我们在求莫比乌斯函数平方和时,要求出莫比乌斯函数的前缀和,也是杜教筛,总时间复杂度 (O(n^{frac{2}{3}}))

求莫比乌斯函数值平方前缀和

(S(n)=sumlimits_{i=1}^n (mu(i))^2 (n<2^{31}))

这个题和之前的就有所不同了。我们令 (f(i)=mu(i)^2),但是要找一个方便计算的 (g) 才行。

我们选取函数 (g(x)=[x=k^2 kin N^+])。我们来算一下 (g)(f) 的狄利克雷卷积。我们会发现 (fast g=1)?!简要证明一下:

((fast g)(n)=sumlimits_{dmid n} g(d)cdot f(frac{n}{d}))。首先分类讨论,如果 (d) 是一个完全平方数,只有当 (d)(n) 的所有约数中最大的完全平方数时,乘积为 (1),否则 (mu(frac{n}{d})=0)。如果 (d) 不是完全平方数,那么 (g(d)=0)。得证。

那我们的计算就简化了不少。

[egin{aligned} g(1)cdot S(n)=sumlimits_{i=1}^n (fast g)(i)-sumlimits_{i=2}^n g(i)cdot S(lfloorfrac{n}{i} floor)\1cdot S(n)=sumlimits_{i=1}^n 1-sumlimits_{i=2}^{sqrt{n}} 1cdot S(lfloorfrac{n}{i^2} floor)\S(n)=n-sumlimits_{i=2}^{sqrt{n}} S(lfloorfrac{n}{i^2} floor) end{aligned}]

我们甚至不用数论分块就可以解决此题。采用同样的优化方法,时间复杂度还是约为 (O(n^{frac{2}{3}}))

相关练习题:

  1. 模板(简单版)
  2. 模板(强化版)
  3. 完全平方数
  4. Sum
  5. Lucas的数论
  6. 自己找