ICode9

精准搜索请尝试: 精确搜索
首页 > 其他分享> 文章详细

浅谈杜教筛

2021-11-01 06:31:16  阅读:176  来源: 互联网

标签:frac 浅谈 limits int sum varphi 杜教 include


今天本身是要学习莫比乌斯反演的,然后题单里面出现了一道需要用杜教筛的题,我还推到了只剩杜教筛的地方

于是迫不得已学习了这种神仙算法,真是迫不得已啊。。。。

杜教筛是一种在线性复杂度以下的求积性函数前缀和的高级算法,大概在\(O(n^{\frac{2}{3}})\)左右

前置知识

积性函数

积性函数:对于任意互质的整数\(a,b\)有\(f(ab)=f(a)f(b)\)则称\(f(x)\)的数论函数

完全积性函数:对于任意整数\(a,b\)有\(f(ab)=f(a)f(b)\)的数论函数

常见的积性函数:\(\varphi,\mu,\sigma,d\)

常见的完全积性函数:\(\epsilon,I,id\)

这里特殊解释一下\(\epsilon,I,id\)分别是什么意思:\(\epsilon(n)=[n=1],I(n)=1,id(n)=n\)

狄利克雷卷积

设\(f,g\)是两个数论函数,它们的狄利克雷卷积卷积是:\((f*g)(n)=\sum\limits_{d|n}f(d)g(\frac{n}{d})\)

性质:满足交换律,结合律

单位元:\(\epsilon\)(即\(f*\epsilon=f\))

结合狄利克雷卷积得到的几个性质:

\(\mu * I = \epsilon\)

\(\varphi * I = id\)

\(\mu * id = \varphi\)

然后一定要记住这些常用的卷积函数,关键时候要反应过来

杜教筛

设\(f\)为积性函数,我们现在要求\(\sum\limits_{i=1}^{n}f(i)\)

线性范围一下显然可以直接筛,但是一般题目要求范围为\(n\in [1,1e10]\)

这样我们需要复杂度在线性以下的算法

设\(S(n)=\sum\limits_{i=1}^{n}f(i)\)

那么再找一个积性函数\(g\),并把它和\(f\)卷起来,考虑这个函数的前缀和

\(\begin{aligned} &= \sum\limits_{i=1}^{n} \sum \limits _{d|i} f(d)g(\frac{i}{d}) \\ &= \sum \limits _{d=1}^{n} g(d)\sum\limits _{i=1}^{\lfloor \frac{n}{d}\rfloor } f(i) \\ &= \sum \limits _{d=1}^{n} g(d) S(\lfloor \frac{n}{d} \rfloor) \end{aligned}\)

考虑\(g(1)S(n)\)就等于

\(\sum\limits_{i=1}^{n}g(i)S(\frac{n}{i})-\sum\limits_{i=2}^{n}g(i)S(\frac{n}{i})\)

若不特殊说明,本文中所写的分数均是向下取整

那么前面的式子可以化简为\(\sum\limits_{i=1}^{n}(f*g)(i)\)

所以得到了\(\huge{大套路式}\)

\(g(1)S(n)=\sum\limits_{i=1}^{n}(f*g)(i)-\sum\limits_{i=2}^{n}g(i)S(\frac{n}{i})\)

这样的话我们剩下的工作就是找到一个合适的\(g\)函数

要求尽量是可以在复杂度\(O(1)orO(\sqrt{n})\)的情况下出\(\sum\limits_{i=1}^{n}(f*g)(i)\),否则无法递归求出解

找到\(g\)后剩下的就是数论分块加递归求解加记忆化了

下面放一个模板题

主要就是使用上述的

\(\varphi * I = id\)

以及

\(\mu * I = \epsilon\)

code
#include<cmath>
#include<cstdio>
#include<cstring>
#include<iostream>
#include<algorithm>
#include<unordered_map>
#define int long long
using namespace std;
const int NN=5000000;
namespace AE86{
	auto read=[](){
		int x=0,f=1;char ch=getchar();
		while(ch<'0'||ch>'9'){if(ch=='-')f=-1;ch=getchar();}
		while(ch>='0'&&ch<='9'){x=(x<<1)+(x<<3)+(ch^48);ch=getchar();}
		return x*f;
	};
}using namespace AE86;
int T,n,mu[NN+5],len,prime[NN+5];
int phi[NN+5];
auto getprime=[](){
	mu[1]=1; phi[1]=1;
	for(int i=2;i<=NN;i++){
		if(!phi[i]) prime[++len]=i,mu[i]=-1,phi[i]=i-1;
		for(int j=1;j<=len&&prime[j]*i<=NN;j++){
			if(i%prime[j]==0){
				phi[i*prime[j]]=phi[i]*prime[j];
				break;
			}
			mu[i*prime[j]]=-mu[i];
			phi[i*prime[j]]=phi[i]*(prime[j]-1);
		}
	}
	for(int i=1;i<=NN;i++) mu[i]+=mu[i-1],phi[i]+=phi[i-1];
};
unordered_map<int,int>smu;
unordered_map<int,int>sphi;
inline int getmu(int n){
	if(n<=NN) return mu[n];
	if(smu.find(n)!=smu.end()) return smu[n];
	int res=1,l=2,r;
	while(l<=n){
		r=n/(n/l);
		res-=(r-l+1)*getmu(n/l);
		l=r+1;
	} return smu[n]=res;
}
inline int getphi(int n){
	if(n<=NN) return phi[n];
	if(sphi.find(n)!=sphi.end()) return sphi[n];
	int res=(1+n)*n/2;int l=2,r;
	while(l<=n){
		r=n/(n/l);
		res-=(r-l+1)*getphi(n/l);
		l=r+1;
	} return sphi[n]=res;
}
namespace WSN{
	inline short main(){
		T=read(); getprime();
		while(T--){
			n=read();
			printf("%lld %lld\n",getphi(n),getmu(n));
		}
		return 0;
	}
}
signed main(){return WSN::main();}

接下来就是今天做莫比乌斯反演时碰见的杜教筛题

简单的数学题

\(\sum\limits_{i=1}^{n}\sum\limits_{j=1}^{n}ij\times gcd(i,j),(\mod{P})\)

\(\begin{aligned} & \sum\limits_{i=1}^{n}\sum\limits_{j=1}^{n}ij\times gcd(i,j)\\ &=\sum\limits_{d=1}^{n}\sum\limits_{i=1}^{\frac{n}{d}}\sum\limits_{j=1}^{\frac{n}{d}}id\times jd\times d[gcd(i,j)=1]\\ &=\sum\limits_{d=1}^{n}d^3\sum\limits_{i=1}^{\frac{n}{d}}\sum\limits_{j=1}^{\frac{n}{d}}i\times j[gcd(i,j)=1] \\ &=\sum\limits_{d=1}^{n}\sum\limits_{i=1}^{\frac{n}{d}}\sum\limits_{j=1}^{\frac{n}{d}}\sum\limits_{k|gcd(i,j)}\mu(k)ij \\ &=\sum\limits_{d=1}^{n}d^3\sum\limits_{k=1}^{\frac{n}{d}}\mu(k)k^2(1+2+...+\frac{\frac{n}{d}}{k})^2 \\ \end{aligned}\)

接下来我们令\(sum(n)=\sum\limits_{i=1}^{n}i\),\(T=kd\)

则有,

\(\begin{aligned} &=\sum\limits_{d=1}^{n}d^3\sum\limits_{k=1}^{\frac{n}{d}}\mu(k)k^2sum^2(\frac{n}{T}) \\ &=\sum\limits_{T=1}^{n}sum^2(\frac{n}{T})\sum\limits_{d|T}d^3\mu(\frac{T}{d})(\frac{T}{d})^2 \\ &=\sum\limits_{T=1}^{n}sum^2(\frac{n}{T})T^2\sum\limits_{d|T}d\mu(\frac{T}{d}) \\ \end{aligned}\)

发现后面的\(\sum\limits_{d|T}d\mu(\frac{T}{d})\)是\(id*\mu=\varphi\)

那么原式可以化为

\(\begin{aligned} &=\sum\limits_{T=1}^{n}sum^2(\frac{n}{T})T^2\varphi(T)\\ &=\sum\limits_{i=1}^{n}sum^2(\frac{n}{i})i^2\varphi(i)\\ \end{aligned}\)

后面的不会了,然后看题解发现要用杜教筛。。

所以学会杜教筛之后这题就成了简单的数学题

然后我们需要找到一个复杂度在线性以下的算法计算\(\sum\limits_{i=1}^{n}i^2\varphi(i)\)

不难考虑杜教筛,那么需要选择一个\(g\),发现可以选择\(g=id^2\),即\(g(n)=n^2\)

那么设\(f(n)=n^2\varphi(n)\),则

\((f*g)(n)=\sum\limits_{d|n}d^2\varphi(d)\times (\frac{n}{d})^2=n^2\sum\limits_{d|n}\varphi(d)\)

发现后面的\(\sum\limits_{d|n}\varphi(d)=(I*\varphi)(n)=id(n)=n\)

那么\((f*g)(n)=n^3\),然后杜教筛出前缀和再数论分块就可以解决了

code
#include<cmath>
#include<cstdio>
#include<cstring>
#include<iostream>
#include<algorithm>
#include<unordered_map>
#define int long long
using namespace std;
const int NN=5000005;
int n,mod,v6,v2,ans;
namespace AE86{
	auto read=[](){
		int x=0,f=1;char ch=getchar();
		while(ch<'0'||ch>'9'){if(ch=='-')f=-1;ch=getchar();}
		while(ch>='0'&&ch<='9'){x=(x<<1)+(x<<3)+(ch^48);ch=getchar();}
		return x*f;
	};
	auto write=[](int x,char opt='\n'){
		char ch[20];short len=0;if(x<0)x=~x+1,putchar('-');
		do{ch[len++]=x%10+(1<<5)+(1<<4);x/=10;}while(x);
		for(short i=len-1;i>=0;--i)putchar(ch[i]);putchar(opt);
	};
	auto qmo=[](int a,int b,int ans=1){
		int c=mod;for(;b;b>>=1,a=a*a%c)if(b&1)ans=ans*a%c;
		return ans;
	};
	inline int sig(int n){return (1+n)%mod*n%mod*v2%mod;}
	inline int pws(int n){return (2*n%mod+1)%mod*n%mod*(n+1)%mod*v6%mod;}
	inline int pps(int n){return sig(n)*sig(n)%mod;}
}using namespace AE86;
namespace Prime{
	signed len,prime[NN];
	int phi[NN],F[NN];
	bool vis[NN];
	auto getprime=[](){
		phi[1]=1; F[1]=1;
		for(signed i=2;i<NN;++i){
			if(!vis[i]) prime[++len]=i,phi[i]=i-1;
			for(signed j=1;j<=len&&prime[j]*i<NN;++j){
				vis[i*prime[j]]=1;
				if(!(i%prime[j])){phi[i*prime[j]]=phi[i]*prime[j];break;}
				phi[i*prime[j]]=phi[i]*(prime[j]-1);
			}
			F[i]=(phi[i]*i%mod*i%mod+F[i-1])%mod;
		}
	};
}using namespace Prime;
namespace dujiao_shai{
	unordered_map<int,int>sF;
	inline int getF(int n){
		if(n<NN) return F[n];
		if(sF[n]) return sF[n];
		int res=pps(n);
		for(int l=2,r;l<=n;l=r+1) r=n/(n/l),
			res=(res-(pws(r)-pws(l-1)+mod)%mod*getF(n/l)%mod+mod)%mod;
		res=(res%mod+mod)%mod;
		return sF[n]=res;
	}
}using namespace dujiao_shai;
namespace WSN{
	inline short main(){
		mod=read();n=read();getprime();
		v6=qmo(6,mod-2);v2=qmo(2,mod-2);
		for(int l=1,r;l<=n;l=r+1) r=n/(n/l),
			ans=(ans+(getF(r)-getF(l-1)+mod)%mod*pps(n/l)%mod)%mod;
		write(ans);
		return 0;
	}
}
signed main(){return WSN::main();}

标签:frac,浅谈,limits,int,sum,varphi,杜教,include
来源: https://www.cnblogs.com/hzoi-wsn/p/15491873.html

本站声明: 1. iCode9 技术分享网(下文简称本站)提供的所有内容,仅供技术学习、探讨和分享;
2. 关于本站的所有留言、评论、转载及引用,纯属内容发起人的个人观点,与本站观点和立场无关;
3. 关于本站的所有言论和文字,纯属内容发起人的个人观点,与本站观点和立场无关;
4. 本站文章均是网友提供,不完全保证技术分享内容的完整性、准确性、时效性、风险性和版权归属;如您发现该文章侵犯了您的权益,可联系我们第一时间进行删除;
5. 本站为非盈利性的个人网站,所有内容不会用来进行牟利,也不会利用任何形式的广告来间接获益,纯粹是为了广大技术爱好者提供技术内容和技术思想的分享性交流网站。

专注分享技术,共同学习,共同进步。侵权联系[81616952@qq.com]

Copyright (C)ICode9.com, All Rights Reserved.

ICode9版权所有