杜教筛

简述

杜教筛是一种小于线性时间复杂度的求积性数论函数前缀和的算法，它利用了积性函数的性质和数论分块将前缀和进行递推。实际上很多数论相关的题目可以在线性时间内跑过，杜教筛也是一种优化的方法，它的基础建立于狄利克雷卷积&莫比乌斯反演之上。

积性函数

定义

我们将对于(forall a,b,aperp b)均有(f(ab)=f(a)f(b))的函数(f)称为积性函数。
像一些常见的数论函数其实都是积性函数，如：(mu,epsilon,varphi,operatorname{id},operatorname{d})等。

性质

若(f,g)为积性函数，则如下都为积性函数：

[f^i(n)\ f(n^i)\ f imes g\ f*g]

杜教筛

设积性函数(f)，另有(s(n)=sum_{i=1}^n f(i))。为了加速，我们需要构造出(s(n))关于(s(lfloorfrac{n}{i} floor))的递推式。
设有(g)为积性函数，(h=f*g)，且(h,g)的前缀和比较方便求出。

[egin{align*} sum_{i=1}^n h(i)&=sum_{i=1}^nsum_{d|i}f(i)g(frac{i}{d})\ &=sum_{d=1}^n g(d)sum_{i=1}^n f(lfloorfrac{i}{d} floor)[d|i]\ &=sum_{d=1}^n g(d)sum_{i=1}^{lfloorfrac{n}{d} floor} f(i)\ &=sum_{d=1}^n g(d)s(lfloorfrac{n}{d} floor)\ Leftrightarrow g(1)s(n)&=sum_{i=1}^n h(i)-sum_{j=2}^n g(j)s(lfloorfrac{n}{j} floor) end{align*}]

基于(h,g)的前缀和比较方便求出，就可以用数论分块来算出(g(j)s(lfloorfrac{n}{j} floor))以求解(s(n))。
(g)的选择是一个略有技巧性的事情，需要对于数论函数的性质有一定的认识和熟练度才方便设出能让式子变得最舒服的函数。

时间复杂度

直接计算时间复杂度为(Θ(n^{frac{3}{4}}))；预处理(Θ(n^{frac{2}{3}}))个前缀和时，时间复杂度为(Θ(n^{frac{2}{3}}))。
证明略去。

例题

1. 求(mu)的前缀和

设此时(g=1)，由( exttt{Dirichlet})卷积可知(h=mu*1=epsilon)，则：

[s(n)=1-sum_{j=2}^n s(lfloorfrac{n}{j} floor) ]

接下来的事就是依上所述，预处理一部分后递推一部分。

2. 求(varphi)的前缀和

有两种方法：直接根据(operatorname{id}=varphi*1)做或根据(varphi=mu*operatorname{id})在(mu)的前缀和基础上做。

第一种

[s_n=frac{(1+n)n}{2} - sum_{j=2}^n s(lfloorfrac{n}{j} floor) ]

第二种

直接套(varphi=mu*operatorname{id})：

[sum_{i=1}^nvarphi(i)=sum_{d=1}^n mu(d)sum_{j=1}^{lfloorfrac{n}{d} floor} j=sum_{d=1}^n mu(d) frac{(1+lfloorfrac{n}{d} floor)lfloorfrac{n}{d} floor}{2} ]

由于(lfloorfrac{n}{d} floor)的取值最多只有(2sqrt{n})种，所以存下了这些(mu)的前缀和之后再数论分块时间复杂度是(Θ(sqrt{n}))，总时间复杂度还是处理(mu)的前缀和的(Θ(n^{frac{2}{3}}))。

luoguP4213 【模板】杜教筛（Sum）

[click]

#include <cstdio>
#include <cctype>
#include <cmath>
#include <map>
using std::map;
typedef long long ll;
const int maxn=1e7+10;
int prim[maxn],s[maxn];
int tot,lim;
bool vis[maxn];
map<int,int> S;

int read()
{
	int res=0;
	char ch=getchar();
	while(!isdigit(ch))
		ch=getchar();
	while(isdigit(ch))
		res=res*10+ch-'0',ch=getchar();
	return res;
}
void prework(int n)
{
	vis[0]=vis[1]=1;
	s[1]=1;
	for (int i=2;i<=n;i++)
	{
		if (!vis[i])
			prim[++tot]=i,s[i]=s[i-1]-1;
		else
			s[i]+=s[i-1];
		for (int j=1;j<=tot&&(ll)i*prim[j]<=n;j++)
		{
			vis[i*prim[j]]=1;
			if (i%prim[j]==0)
				break;
			s[i*prim[j]]=-(s[i]-s[i-1]);
		}
	}
}
int solve(int n)
{
	if (n<=lim)
		return s[n];
	else
	{
		map<int,int>::iterator p=S.find(n);
		if (p!=S.end())
			return (*p).second;
	}
	int res=1; 
	for (int l=2,r=1;l<=n;l=r+1)
	{
		r=n/(n/l);
		res-=(r-l+1)*solve(n/l);
		if (r==n)
			break; 
	}
	return S[n]=res;
}
int main()
{
	int T=read();
	lim=pow((1LL<<31)-1, 2.0/3.0);
	prework(lim);
	while(T--)
	{
		int n=read();
		ll ans1=0;
		int ans2=solve(n);
		for (int l=1,r=1;l<=n;l=r+1)
		{
			r=n/(n/l);
			ans1+=(ll)(1+n/l)*(n/l)/2*(solve(r)-solve(l-1));
			if (r==n)
				break;
		}
		printf("%lld %d
",ans1,ans2);
	}
	return 0;
}

参考

• OI WIKI