杜教筛小结

前置芝士

狄利克雷卷积

进入正题

现在要求一个函数 $f$ 的前缀和 $S$。

但是往往连 $O(n)$ 的算法都满足不了毒瘤的出题人，所以需要一个更快的做法，于是时间复杂度 $O(n^{\frac{2}{3}})$ 的杜教筛出现了。

具体做法：

找两个积性函数 $g,h$，满足 $f\ast g=h$。

那么有：
$$\begin{gathered} h(i)=\sum _{d|i}f(d)\times g(\frac{i}{d}) \\ \sum _{i=1}^{n}h(i)=\sum _{i=1}^n\sum _{d|i}f(d)\times g(\frac{i}{d}) \\ \sum _{i=1}^{n}h(i)=\sum _{i=1}^{n}g(i)\sum _{j=1}^{\lfloor \frac{n}{i}\rfloor }f(j) \\ \sum _{i=1}^{n}h(i)=\sum _{i=1}^{n}g(i)S(\lfloor \frac{n}{i}\rfloor ) \end{gathered}$$

把右边的第一项提出来，就有：
$$\begin{gathered} \sum _{i=1}^{n}h(i)=S(n)g(1)+\sum _{i=2}^{n}g(i)S(\lfloor \frac{n}{i}\rfloor ) \\ S(n)g(1)=\sum _{i=1}^{n}h(i)-\sum _{i=2}^{n}g(i)S(\lfloor \frac{n}{i}\rfloor ) \end{gathered}$$

那么如果 $h$ 和 $g$ 的前缀和很好求，那么用一个除法分块加递归就可以做了。

需要注意的是，递归里面需要对每一次求出来的前缀和记忆化，这个可以用哈希或map实现，这样时间复杂度就是 $O(n^{\frac{3}{4}})$。

你可能会问，这样递归下去求出来的前缀和数量不会很多吗？事实上只有 $2\sqrt{n}$ 个，也就是不同的 $\lfloor \frac{n}{i}\rfloor$ 的个数，因为每一个求到的前缀和都可以表示为 $\lfloor \frac{n}{x_1{x}_2...}\rfloor$ 的形式。

假如线性筛预处理 $S$ 的一部分，那么上面的杜教筛的时间复杂度就可以更加优秀，比如预处理前 $K$ 个，根据一些我不会的微积分知识可以算出复杂度会变成 $O(K+\frac{n}{\sqrt{K}})$，取 $K=n^{\frac{2}{3}}$ 时有最优复杂度，即 $O(n^{\frac{2}{3}})$。

实例

拿筛 $\mu$ 和筛 $\phi$ 做例子。

筛莫比乌斯函数的前缀

设 $f=\mu$。

因为 $\mu$ 有一个很优秀的性质——${\sum }_{d|n}\mu (d)=[n=1]$，写成狄利克雷卷积的形式就是： $f\ast I=ϵ$。（不认识这些函数的话请将鼠标滚轮上滑，点开狄利克雷卷积那篇博客）

那么此时使函数 $g=I,h=ϵ$，就可以做杜教筛了。

把上面的柿子拿来，就是：
$$S(n)g(1)=\sum _{i=1}^{n}h(i)-\sum _{i=2}^{n}g(i)S(\lfloor \frac{n}{i}\rfloor )$$

显然等于：
$$S(n)=1-\sum _{i=2}^{n}1\times S(\lfloor \frac{n}{i}\rfloor )$$

代码实现：

#define maxn 6000010
int t;
ll n;
ll mu[maxn];
int prime[maxn],tt=0;
bool v[maxn];
void work()//预处理mu的前6000000位的前缀和
{
	mu[1]=1ll;
	for(int i=2;i<=maxn-10;i++)
	{
		if(!v[i])prime[++tt]=i,mu[i]=-1;
		for(int j=1;j<=tt&&i*prime[j]<=maxn-10;j++)
		{
			v[i*prime[j]]=true;
			if(i%prime[j]==0)
			{
				mu[i*prime[j]]=0;
				break;
			}
			mu[i*prime[j]]=-mu[i];
		}
	}
	for(int i=2;i<=maxn-10;i++)
	mu[i]+=mu[i-1];
}
map<int,int> mapmu;
int solvemu(int x)
{
	if(x<=maxn-10)return mu[x];//假如在6000000位之前
	if(mapmu[x])return mapmu[x];//假如之前算过
	int re=1,l=2,r;
	while(l<=x)//除法分块
	{
		r=x/(x/l);
		re-=(r-l+1)*solvemu(x/l);
		l=r+1;
	}
	mapmu[x]=re;//记忆化
	return re;
}

筛欧拉函数的前缀和

因为欧拉函数有这样一个优秀的性质：${\sum }_{d|n}\phi (d)=n$，所以有：$\phi \ast I=id$。

于是使$g=I$，$h=id$，用上面的公式套下来就是：
$$S(n)=\sum _{i=1}^{n}h(i)-\sum _{i=2}^{n}1\times S(\lfloor \frac{n}{i}\rfloor )$$

$h$ 显然是个等差数列，可以 $O(1)$ 求。

代码和上面类似，就不贴了。

顺手补上模板题传送门。

代码如下：

#include <cstdio>
#include <cstring>
#include <map>
#include <algorithm>
using namespace std;
#define ll long long
#define maxn 4000010
#define uint unsigned int

int t;
uint n;
ll mu[maxn],phi[maxn];
int prime[maxn],tt=0;
bool v[maxn];
void work()
{
	mu[1]=phi[1]=1ll;
	for(int i=2;i<=maxn-10;i++)
	{
		if(!v[i])prime[++tt]=i,mu[i]=-1,phi[i]=i-1;
		for(int j=1;j<=tt&&i*prime[j]<=maxn-10;j++)
		{
			v[i*prime[j]]=true;
			if(i%prime[j]==0)
			{
				mu[i*prime[j]]=0;
				phi[i*prime[j]]=phi[i]*(ll)prime[j];
				break;
			}
			mu[i*prime[j]]=-mu[i];
			phi[i*prime[j]]=phi[i]*phi[prime[j]];
		}
	}
	for(int i=2;i<=maxn-10;i++)
	mu[i]+=mu[i-1],phi[i]+=phi[i-1];
}
map<uint,ll> mapphi;
int tot=0;
ll solvephi(uint x)
{
	if(x<=maxn-10)return phi[x];
	if(mapphi[x])return mapphi[x];
	ll re=(ll)(1+x)*x/2ll;
	uint l=2,r;
	while(l<=x)
	{
		r=x/(x/l);
		re-=(ll)(r-l+1)*solvephi(x/l);
		l=r+1;
	}
	mapphi[x]=re;
	return re;
}
map<int,ll>mapmu;
ll solvemu(uint x)
{
	if(x<=maxn-10)return mu[x];
	if(mapmu[x])return mapmu[x];
	ll re=1;
	uint l=2,r;
	while(l<=x)
	{
		r=x/(x/l);
		re-=(ll)(r-l+1)*solvemu(x/l);
		l=r+1;
	}
	mapmu[x]=re;
	return re;
}

int main()
{
	work();
	scanf("%d",&t);
	while(t--)
	{
		scanf("%u",&n);
		printf("%lld %lld\n",solvephi(n),solvemu(n));
	}
}

但事实上这道模板题还可以做的更优。

利用 $\phi$ 和 $\mu$ 之间的一个性质：$\phi =\mu \ast id$，那么 $\phi$ 的前缀和就可以写成：
$$\begin{gathered} \sum _{i=1}^n\phi (i)=\sum _{i=1}^n\sum _{d|i}\mu (d)\times \frac{i}{d} \\ =\sum _{i=1}^{n}i\sum _{j=1}^{\lfloor \frac{n}{i}\rfloor }\mu (j) \\ =\sum _{i=1}^{n}i\times S(\lfloor \frac{n}{i}\rfloor ) \end{gathered}$$

这样就可以用除法分块来做，转化成了求 $\mu$ 的前缀和，因为 $\mu$ 的前缀和我们有记忆化，所以此时跑起来飞快。

代码如下：

#include <cstdio>
#include <cstring>
#include <map>
#include <algorithm>
using namespace std;
#define maxn 5000010
#define uint unsigned int
#define ll long long

uint T,n;ll ans1,ans2;
int prime[maxn],t=0;
ll mu[maxn];
bool v[maxn];void work()
{
	mu[1]=1;
	for(int i=2;i<=maxn-10;i++)
	{
		if(!v[i])prime[++t]=i,mu[i]=-1;
		for(int j=1;j<=t&&i*prime[j]<=maxn-10;j++){
			v[i*prime[j]]=true;
			if(i%prime[j]==0)break;
			mu[i*prime[j]]=-mu[i];
		}
		mu[i]+=mu[i-1];
	}
}
map<uint,ll>S;
ll solvemu(uint x)
{
	if(x<=maxn-10)return mu[x];
	if(S.count(x))return S[x];
	ll re=1;
	for(uint l=2,r;l<=x;l=r+1)
	r=x/(x/l),re-=(r-l+1)*solvemu(x/l);
	S[x]=re; return re;
}
ll solvephi(uint x)
{
	ll re=0;
	for(uint l=1,r;l<=x;l=r+1)
	r=x/(x/l),re+=1ll*(l+r)*(r-l+1)/2*solvemu(x/l);
	return re;
}

int main()
{
	scanf("%u",&T);work();
	while(T--){
		scanf("%u",&n);
		ans2=solvemu(n);
		ans1=solvephi(n);
		printf("%lld %lld\n",ans1,ans2);
	}
}

但是，这还不是杜教筛的时间复杂度最玄学的地方。

比如说现在有一个函数 $f$，需要求下面这个柿子的值：
$$\sum _{i=1}^{n}f(\lfloor \frac{n}{i}\rfloor )$$

显然是一个除法分块。

但是如果 $n\le 10^{10}$，$f$ 无法预处理，那么就做不了了。

考虑杜教筛。

假如写除法分块，那么里面求值的语句是长这样子的：

ans+=(r-l+1)*f(n/l);

现在的问题是 $f(n/l)$ 求不了。

设 $S(n)={\sum }_{i=1}^{n}f(i)$，那么 $f(n/l)=S(n/l)-S(n/l-1)$。

没错，你肯定猜到了，这样转化一下之后，后面的 $S$ 是个前缀和，我们用杜教筛来搞他。

于是这个$f(n/l)$的问题就解决了。

不知道从哪里蹦出来的读者： 可是你这样的时间复杂度不是 $O(\sqrt{n}\times n^{\frac{2}{3}})$的吗？更慢了！

看起来是很慢的，但是实际上，由于我们对 $S(n)$ 求杜教筛后，$S(n/l)$ 的前缀和也都被求出来了，于是后面用到的时候复杂度就变成了 $O(1)$ 而不需要重新筛。

所以实际上，这样做的时间复杂度就是 $O(n^{\frac{2}{3}}+\sqrt{n})$ 的。

这个技巧比如这题就有用到。

另外一道题：

简单的数学题   题解