杜教筛学习笔记

时间:2024-01-27 14:10:13

杜教筛学习笔记

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

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

应用

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

方法

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

\[\sum\limits_{i=1}^n (f\ast g)(i)=\sum\limits_{i=1}^n g(i)\cdot S(\lfloor\frac{n}{i}\rfloor) \]

我们给个证明:

\[\begin{aligned} \sum\limits_{i=1}^n (f\ast g)(i) \\=\sum\limits_{i=1}^n\sum\limits_{d\mid i} f(d)\cdot g(\frac{i}{d})\\=\sum\limits_{d=1}^n\sum\limits_{i=1}^{\lfloor \frac{n}{d}\rfloor} g(d)\cdot f(i)\\=\sum\limits_{d=1}^n g(d)\cdot\sum\limits_{i=1}^{\lfloor \frac{n}{d}\rfloor} f(i)\\=\sum\limits_{d=1}^n g(d)\cdot S(\lfloor \frac{n}{d}\rfloor)\\=\sum\limits_{i=1}^n g(i)\cdot S(\lfloor \frac{n}{i}\rfloor)\\&\square \end{aligned}\]

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

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

\[\begin{aligned} \sum\limits_{i=1}^n (f\ast g)(i)=\sum\limits_{i=1}^n g(i)\cdot S(\lfloor\frac{n}{i}\rfloor)\\\sum\limits_{i=1}^n (f\ast g)(i)=g(1)\cdot S(n)+\sum\limits_{i=2}^n g(i)\cdot S(\lfloor\frac{n}{i}\rfloor)\\g(1)\cdot S(n)=\sum\limits_{i=1}^n (f\ast g)(i)-\sum\limits_{i=2}^n g(i)\cdot S(\lfloor\frac{n}{i}\rfloor) \end{aligned}\]

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

例子

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

求莫比乌斯函数前缀和

\(S(n)=\sum\limits_{i=1}^n \mu(i)\ \ \ \ (n<2^{31})\)

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

\[\begin{aligned} g(1)\cdot S(n)=\sum\limits_{i=1}^n (f\ast g)(i)-\sum\limits_{i=2}^n g(i)\cdot S(\lfloor\frac{n}{i}\rfloor)\\1\cdot S(n)=\sum\limits_{i=1}^n \varepsilon(i)-\sum\limits_{i=2}^n 1\cdot S(\lfloor\frac{n}{i}\rfloor)\\S(n)=1-\sum\limits_{i=2}^n S(\lfloor\frac{n}{i}\rfloor) \end{aligned}\]

根据数论分块的知识,\(\lfloor\frac{n}{i}\rfloor\) 的取值有 \(O(\sqrt{n})\) 个。那么时间复杂度为 \(O(\sum\limits_{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\n",pre(read()));
	}
	return 0;
}

求欧拉函数前缀和

\(S(n)=\sum\limits_{i=1}^n \varphi(i)\ \ \ \ (n<2^{31})\)

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

\[\begin{aligned} g(1)\cdot S(n)=\sum\limits_{i=1}^n (f\ast g)(i)-\sum\limits_{i=2}^n g(i)\cdot S(\lfloor\frac{n}{i}\rfloor)\\1\cdot S(n)=\sum\limits_{i=1}^n id(i)-\sum\limits_{i=2}^n 1\cdot S(\lfloor\frac{n}{i}\rfloor)\\S(n)=\frac{n^2+n}{2}-\sum\limits_{i=2}^n S(\lfloor\frac{n}{i}\rfloor) \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\n",pre(read()));
	}
	return 0;
}

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

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

\(S(n)=\sum\limits_{i=1}^n \mu^2(i)\ \ \ \ (n<2^{31})\)

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

\[\begin{aligned} g(1)\cdot S(n)=\sum\limits_{i=1}^n (f\ast g)(i)-\sum\limits_{i=2}^n g(i)\cdot S(\lfloor\frac{n}{i}\rfloor)\\1\cdot S(n)=\sum\limits_{i=1}^n \mu(i)-\sum\limits_{i=2}^n 1\cdot S(\lfloor\frac{n}{i}\rfloor)\\S(n)=\sum\limits_{i=1}^n \mu(i)-\sum\limits_{i=2}^n S(\lfloor\frac{n}{i}\rfloor) \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!!!
//#pragma GCC optimize("Ofast","-funroll-loops","-fdelete-null-pointer-checks")
//#pragma GCC target("ssse3","sse3","sse2","sse","avx2","avx")
#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=4e6+10;

int n,mu[maxn],p[maxn],cnt,s[maxn],s2[maxn];
bool is[maxn];
map<long,long> mp,mp2;

void sieve_mu(int n) {
	fill(is,is+n+1,1);
	mu[1]=1;
	
	for(int i=2;i<=n;i++) {
		if(is[i]) {
			cnt++;
			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];
		}
	}
	s[0]=s2[0]=0;
	for(int i=1;i<=n;i++){
		s[i]=s[i-1]+mu[i];
		s2[i]=s2[i-1]+mu[i]*mu[i];
	}
}

int premu(int n) {
	if(n<=maxn-10) {
		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)*premu(n/i);
	}
	return mp[n]=ret;
}

int pre(int n) {
	if(n<=maxn-10) {
		return s2[n];
	}
	if(mp2[n]) {
		return mp2[n];
	}
	int j=2,ret=premu(n);
	for(int i=2;i<=n;i=j+1) {
		j=(n/(n/i));
		
		ret-=(j-i+1)*pre(n/i);
	}
	return mp2[n]=ret;
}

signed main() {
	//freopen(".in","r",stdin);
	//freopen(".out","w",stdout);
	
	sieve_mu(maxn-10);
	int t=read();
	
	while(t--) {
		printf("%lld\n",pre(read()));
	}

	//fclose(stdin);
	//fclose(stdout);
	return 0;
}

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

\(S(n)=\sum\limits_{i=1}^n (\mu(i))^2\ \ \ \ (n<2^{31})\)

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

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

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

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

\[\begin{aligned} g(1)\cdot S(n)=\sum\limits_{i=1}^n (f\ast g)(i)-\sum\limits_{i=2}^n g(i)\cdot S(\lfloor\frac{n}{i}\rfloor)\\1\cdot S(n)=\sum\limits_{i=1}^n 1-\sum\limits_{i=2}^{\sqrt{n}} 1\cdot S(\lfloor\frac{n}{i^2}\rfloor)\\S(n)=n-\sum\limits_{i=2}^{\sqrt{n}} S(\lfloor\frac{n}{i^2}\rfloor) \end{aligned}\]

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

相关练习题:

  1. 模板(简单版)
  2. 模板(强化版)
  3. 完全平方数