浅(kou)谈(hu)杜教筛

引入

让我们从一道例题开始
这道题目需要我们求出\(\varphi(n)\)\(\mu(n)\)这两个函数的前缀和
显然,这可以通过线性筛\(\mathcal O(n)\)完成,但这道题目的\(n\)的范围是\(\le 2^{31}\),线性筛无法完成,我们需要一种更高效的算法————杜教筛
test


test

前置知识

1.线性筛
2.积性函数:如果函数\(f(n)\)满足当\(gcd(n,m)=1\)\(f(nm)=f(n)f(m)\),那么我们称\(f(n)\)为积性函数。
3.完全积性函数:就是积性函数去掉那个互质的条件。
4.艾弗森约定:对于布尔型变量\(x\)\([x]\)表示\(x\)为真时取\(1\),为假时取\(0\)
5.一些常见积性函数

  • 元函数\(\epsilon(n)=[n=1]\),是完全积性的。
  • 幂函数\(id_k(n)=n^k\),是完全积性的。
  • 单位函数\(I(n)=id_0(n)\)
  • 约数幂函数\(\sigma_k(n)=\sum_{d|n}d^k\),是积性的。
  • 欧拉函数\(\varphi(n)=\sum_{d≤n}[gcd(d,n)=1]\),是积性的.
  • 莫比乌斯函数\(μ(n)\)满足\(\sum_{d|n}μ(d)=[n=1]=\epsilon(n)\)是积性的

6.狄利克雷卷积
\(f,g\)是两个数论函数,它们的狄利克雷卷积是:\((f \otimes g)(n)=\sum_{d|n} f(d) g(\frac nd)\)
性质:满足交换律,结合律


概述

杜教筛,因由圈内著名的杜老师(杜瑜皓)引进而得名,可以以\(\mathcal O(n^{\frac 23})\)的时间复杂度求解一类积性函数的前缀和
假设\(f(n)\)是我们需要求出的积性函数,令

\[S(n)=\sum_{i=1}^{n} f(i) \]

我们再寻找一个积性函数\(g(n)\),令\(h(n)=(f \otimes g)(n)\),那么

\[\sum_{i=1}^{n} h(i) =\sum_{i=1}^{n} \sum_{d|i} f(d)g(\frac id) \]

\[=\sum_{d=1}^{n} \sum_{i=1}^{\lceil \frac nd \rceil} g(d)f(i) \]

\[=\sum_{d=1}^{n} g(d) \sum_{i=1}^{\lceil \frac nd \rceil} f(i) \]

\[=\sum_{d=1}^{n} g(d)S(\lceil \frac nd \rceil) \]

移项后得到

\[g(1)S(n)=\sum_{i=1}^{n} h(i)-\sum_{d=2}^{n}g(d)S(\lceil \frac nd \rceil) \]

那么,如果我们能够较快地求出\(g(n)\)\(h(n)\)的前缀和,那么减号右边的柿子可以整除分块后递归处理,我们就能够在低于线性的复杂度下完成对\(S(n)\)的求解
可以证明直接递归并使用\(unordered\_map\)或哈希表记忆化后该算法的复杂度是\(\mathcal O(n^{\frac 34})\)
如果在此基础上,先线性筛求出前\(n^{\frac 23}\)\(f\)的前缀和,那么我们就可以进一步优化到\(\mathcal O(n^{\frac 23})\)
这就是杜教筛的核心思想


例题

洛谷P4213【模板】杜教筛(Sum)
就是引入中的那道题。。。
sol.
根据经典卷积柿子\(\mu \otimes I=\epsilon,\varphi \otimes I=id_1\)
\(\epsilon,I,id_1\)的前缀和都可以\(\mathcal O(1)\)求出
于是直接上模板就行了

#include<bits/stdc++.h>
using namespace std;
const int N=5000000;
typedef long long ll;
typedef unsigned int uint;
int T;
uint n;
ll phi[N+10],mu[N+10];
unordered_map<uint,ll> ans1,ans2;
inline ll sumphi(uint m){ 
	if(m<=N) return phi[m];
	if(ans1.count(m)) return ans1[m];
	ll ret=1ll*m*(m+1)/2;uint r;
	for(uint l=2;l<=m;l=r+1){
		r=m/(m/l);
		ret-=1ll*(r-l+1)*sumphi(m/l);
	}
	return ans1[m]=ret;
}
inline ll summu(uint m){ 
	if(m<=N) return mu[m];
	if(ans2.count(m)) return ans2[m];
	ll ret=1;uint r;
	for(uint l=2;l<=m;l=r+1){
		r=m/(m/l);
		ret-=1ll*(r-l+1)*summu(m/l);
	}
	return ans2[m]=ret;
}
int p[N+10],f[N+10],pcnt;
inline void init(){
	mu[1]=phi[1]=1;
	for(int i=2;i<=N;++i){
		if(!f[i]) p[++pcnt]=i,phi[i]=i-1,mu[i]=-1;
		for(int j=1;j<=pcnt&&i*p[j]<=N;++j){
			int d=i*p[j];f[d]=1;
			if(i%p[j]) mu[d]=-mu[i],phi[d]=phi[i]*phi[p[j]];
			else{
				mu[d]=0;phi[d]=phi[i]*p[j];
				break;
			}
		}
	}
	for(int i=2;i<=N;++i) phi[i]+=phi[i-1],mu[i]+=mu[i-1];
} 
int main(){
	scanf("%d",&T);
	init();
	while(T--){
		scanf("%u",&n);
		printf("%lld %lld\n",sumphi(n),summu(n));
	}
	return 0;
}

洛谷P3768 简单的数学题
题目要求

\[\sum_{i=1}^{n} \sum_{j=1}^{n}ij gcd(i,j) \]

sol.
因为

\[\sum_{d|n} \varphi(d)=n \]

\[原式=\sum_{i=1}^{n} \sum_{j=1}^{n} ij \sum_{d|n} \varphi(d) \]

\[=\sum_{d=1}^{n} \varphi(d) \sum_{d|i} i\sum_{d|j} j \]

\[=\sum_{d=1}^{n} \varphi(d)d^2 \sum_{i=1}^{\lceil \frac nd \rceil} i \sum_{j=1}^{\lceil \frac nd \rceil} j \]

\[=\sum_{d=1}^{n} \varphi(d)d^2 (\sum_{i=1}^{\lceil \frac nd \rceil} i)^2 \]

其中\((\sum_{i=1}^{\lceil \frac nd \rceil} i)^2\)可以数论分块后\(\mathcal O(\sqrt{n})\)求出
于是我们现在就只需要求\(f(x)=\varphi(x)x^2\)的前缀和了
因为

\[\sum_{d|n}\varphi(d)=n \]

又因为\(d*\frac nd=n\),于是

\[\sum_{d|n}\varphi(d)d^2 (\frac nd)^2 =n^3 \]

于是令\(g=n^2,h=n^3\)即可
显然\(g\)的前缀和\(=\frac {n(n+1)(2*n+1)}{6}\),\(h\)的前缀和\(=(\frac{n(n+1)}{2})^2\),可以\(\mathcal O(1)\)求出
于是上杜教筛即可

#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
const int N=5e6;
int phi[N+10],f[N+10],p[N+10],pcnt,mod;
ll n;
inline int add(int x,int y){return (x+y>=mod)?x+y-mod:x+y;}
inline int dec(int x,int y){return (x-y<0)?x-y+mod:x-y;}
inline void init(){
	phi[1]=1;
	for(int i=2;i<=N;++i){
		if(!f[i]) p[++pcnt]=i,phi[i]=i-1;
		for(int j=1;j<=pcnt&&i*p[j]<=N;++j){
			f[i*p[j]]=1;
			if(i%p[j]) phi[i*p[j]]=1ll*phi[i]*phi[p[j]]%mod;
			else{
				phi[i*p[j]]=1ll*phi[i]*p[j]%mod;
				break;
			}
		}
	}
	for(int i=1;i<=N;++i) f[i]=add(1ll*phi[i]*i%mod*i%mod,f[i-1]);
}
unordered_map<ll,int> ans;
inline int ksm(int x,int y){
	int ret=1;
	for(;y;x=1ll*x*x%mod,y>>=1) if(y&1) ret=1ll*ret*x%mod;
	return ret;
}
int iv2,iv6;
inline int sum(ll x){x%=mod;return 1ll*x*(x+1)%mod*iv2%mod;}
inline int sumpow(ll x){x%=mod;return 1ll*x*(x+1)%mod*(2*x+1)%mod*iv6%mod;}
inline int sumf(ll x){
	if(x<=N) return f[x];
	if(ans.count(x)) return ans[x];
	ll ret=1ll*sum(x)*sum(x)%mod,r;
	for(ll l=2;l<=x;l=r+1){
		r=x/(x/l);
		ret=dec(ret,1ll*dec(sumpow(r),sumpow(l-1))*sumf(x/l)%mod);
	}
	ans[x]=ret;
	return ret;
}
int main(){
//	freopen("lx.out","w",stdout);
	scanf("%d%lld",&mod,&n);
	init();
	int ret=0;iv2=ksm(2,mod-2);iv6=ksm(6,mod-2);ll r;
	for(ll l=1;l<=n;l=r+1){
		r=n/(n/l);
		ret=add(ret,1ll*dec(sumf(r),sumf(l-1))*sum(n/l)%mod*sum(n/l)%mod);
	}
	printf("%d\n",ret);
	return 0;
}

想必你也注意到了,我们的杜教筛往往筛的是一些与\(\varphi,\mu\)等积性函数相关的函数,于是我们这里给出一些常见的卷积,其他的柿子大部分都能由这些导出

\[\mu \otimes I=\epsilon \]

\[\phi \otimes I=id_1 \]

\[id_k \otimes I=\sigma_k \]

\[f \otimes \epsilon=f \]

原文地址:https://www.cnblogs.com/tqxboomzero/p/14039589.html