杜教筛

我也不知现在我学这个干什么但是真的闲着没事干又不想卷whk所以就来了。

首先前置芝士懵逼钨丝反演我之前已经差不多写过一个,但是有些结论建议看看oiwiki什么的反正我也没写。

首先杜教筛大家知道是个亚线性复杂度求积性函数前缀和的东西。于是现在我们有一个积性函数\(f(i)\),要求它的前缀和\(S(n)=\sum_{i=1}^nf(i)\),范围随便给,大概给个1e9差不多了。

先不管别的有的没的,我们先随便另外找一个积性函数\(g\),那么有:

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

\[=\sum_{i=1}^n\sum_{j=1}^{\lfloor\frac ni\rfloor}g(i)f(j) \]

\[=\sum_{i=1}^ng(i)\sum_{j=1}^{\lfloor\frac ni\rfloor}f(j) \]

\[=\sum_{i=1}^ng(i)S(\lfloor\frac ni\rfloor) \]

于是我们把\(g(1)S(n)\)这一项扔到左边,就有了如下的柿子:

\[g(1)S(n)=\sum_{d|n}f(d)g(\frac nd)-\sum_{i=2}^ng(i)S(\lfloor\frac ni\rfloor) \]

发现我们要的\(S(n)\)就在左边,于是这就是杜教筛的基本原理。也就是说,只要我们可以快速求出\(f,g\)的卷积,就可以通过递归的方式来筛出前缀和。

时间复杂度是\(O(n^{\frac 34})\),如果预处理出前\(n^{\frac 23}\)项的前缀和那就是\(O(n^{\frac 23})\)

接下来是几个例题:

P4213 【模板】杜教筛

大概是求\(\mu\)\(φ\)的前缀和。也就是说我们要构造一个函数\(g\)使得它比较好算。

根据我们的经验,我们有:

\[\mu * I=\epsilon \]

\[φ * I=\operatorname{id} \]

其中\(\epsilon=[n=1],\operatorname{id}=n,I=1\)

于是就很显然了。还有一个小优化是记忆化,即用hash存下筛过的所有数。我这份板子是用的unordered_map,虽然跑的比手写hash慢不止一点半点但是好写。

void get(){
	miu[1]=1;phi[1]=1;
	for(int i=2;i<=5000000;i++){
		if(!v[i]){
			p[++p[0]]=i;miu[i]=-1;phi[i]=i-1;
		}
		for(int j=1;j<=p[0]&&i*p[j]<=5000000;j++){
			v[i*p[j]]=true;
			if(i%p[j]==0){
				phi[i*p[j]]=phi[i]*p[j];break;
			}
			phi[i*p[j]]=phi[i]*(p[j]-1);
			miu[i*p[j]]=-miu[i];
		}
	}
	for(int i=1;i<=5000000;i++){
		miu[i]+=miu[i-1];phi[i]+=phi[i-1];
	}
}
int getphi(int x){
	if(x<=5000000)return phi[x];
	if(ans1[x]!=0)return ans1[x];
	int ans=(x+1)*x/2;//卷积是n
	for(int l=2,r;l<=x;l=r+1){
		r=x/(x/l);
		ans-=(r-l+1)*getphi(x/l);
	}
	ans1[x]=ans;
	return ans;
}
int getmiu(int x){
	if(x<=5000000)return miu[x];
	if(ans2[x]!=0)return ans2[x];
	int ans=1;//卷积是1
	for(int l=2,r;l<=x;l=r+1){
		r=x/(x/l);
		ans-=(r-l+1)*getmiu(x/l);
	}
	ans2[x]=ans;
	return ans;
}

第二个:求\(\sum_{i=1}^nφ(i^2)\)

首先设\(f(i)=φ(i^2)=i\cdotφ(i)\)

通常来说我们可以把它的卷积形式拆开来观察:

\[\sum_{d|n}d\cdotφ(d) \]

我们好像可以整点东西把它消掉。

\[\sum_{d|n}d\cdotφ(d)\cdot\frac nd\cdot I(\frac nd) \]

\[=n\cdot (φ * I) \]

\[=n^2 \]

而且我们有:

\[\sum_{i=1}^ni^2=\frac {n(n+1)(2n+1)}6 \]

于是结束了。代码:

int getf(int x){
	if(x<=200000)return f[x];
	if(mp[x]!=0)return mp[x];
	int ans=1ll*x*(x+1)%mod*(2*x+1)%mod*inv%mod;
	for(int l=2,r;l<=x;l=r+1){
		r=x/(x/l);
		ans=(ans-1ll*(r-l+1)*(l+r)/2%mod*getf(x/l)%mod+mod)%mod;
	}
	mp[x]=ans;
	return ans;
}
posted @ 2022-09-03 11:46  gtm1514  阅读(8)  评论(0编辑  收藏  举报