心得技巧
html5HTML/XhtmlCSSXML/XSLTDreamweaver教程Frontpage教程心得技巧

杜教筛（基础篇）

发布时间：2022-06-21 发布网站：脚本宝典
脚本宝典收集整理的这篇文章主要介绍了杜教筛（基础篇），脚本宝典觉得挺不错的，现在分享给大家，也给大家做个参考。

前置知识：莫比乌斯反演专题(基础篇)

杜教筛

设(f(n))为一积性函数 求(S(n)=sum_{i=1}^nf(i))，(nleq 10^{10}) 我们考虑给(f(n))卷上另一个积性函数(f(n))再做前缀和

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

[sum_{d=1}^ng(d)sum_{d|i}f(frac{i}{d}) ]

[sum_{d=1}^ng(d)sum_{d|i}f(frac{i}{d}) ]

[sum_{d=1}^ng(d)sum_{i=1}^{lfloor frac{n}{d}rfloor}f(i) ]

[sum_{d=1}^ng(d)S(lfloor frac{n}{d}rfloor) ]

这样就得到了一个子问题 接着我们特别看(d=1)这一项，也就是(g(1)times S(n))

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

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

这个式子中，(g(1))很容易求，(g)和(f)卷积的前缀和可以通过构造合适的(g)让他也很容易求，后边的可以通过整除分块+递归解决 如果我们事先预处理除一部分的(S)那么可以证明计算(S(n))的复杂度约为(O(n^{frac{2}{3}}))

[BZOJ3944]--Sum

最常见的

[sum_{i=1}^{n}mu(i) ]

[sum_{i=1}^{n}varphi(i) ]

先看(mu(n))，我们的任务是构造一个合适的(g) 我们知道(mu times I = e)（(times) 表示狄利克雷卷积） 因此我们构造(g(n)=I(n)=1) 那么(g)与(mu)的前缀和还是(e)，因此这一块的和为1 也就是说

[S(n)=frac{1-sum_{d=2}^nS(lfloor frac{n}{d}rfloor)}{1} ]

可以杜教筛解决 再看(varphi(n))，我们知道(varphi(n)times I=id) 那么我们还考虑卷上(g(n)=I(n)=1) 那么式子的值为

[S(n)=frac{sum_{i=1}^ni-sum_{d=2}^nS(lfloor frac{n}{d}rfloor)}{1} ]

左边是等差数列求和很容易算

#include<bits/stdc++.h>
using namespace std;
typedef long long LL;
const int N = 4641600;
LL s[N];
LL mu[N];
int v[N],prime[N],tot=0;
LL phi[N];
LL p[N];
void init(int n)
{
	phi[1]=1;
	mu[1]=1;
	for(int i=2;i<=n;i++)
	{
		if(!v[i])
		{
			prime[++tot]=i;
			v[i]=i;
			mu[i]=-1;
			phi[i]=i-1;
		}
		for(int j=1;j<=tot;j++)
		{
			if(prime[j]>v[i]||i*prime[j]>n) break;
			v[i*prime[j]]=prime[j];
			if(i%prime[j]==0)
			{
				mu[i*prime[j]]=0;
				phi[i*prime[j]]=phi[i]*prime[j];
				break;
			}
			else 
			{
				mu[i*prime[j]]=mu[i]*mu[prime[j]];
				phi[i*prime[j]]=1ll*phi[i]*(prime[j]-1);
			}
		}
	}
	for(int i=1;i<=n;i++)
	mu[i]=mu[i-1]+mu[i];
	for(int i=1;i<=n;i++)
	phi[i]=phi[i-1]+phi[i];
}
bool vis[N];
LL S(LL n,int k)
{
	if(n<=N-10) return mu[n];
	if(vis[k]) return s[k];
	LL res=1;
	LL l=2,r;
	for(;l<=n;l=r+1)
	{
		r=n/(n/l);
		res=res-1ll*(r-l+1)*S(n/l,1ll*k*l);
	}
	vis[k]=1;
	s[k]=res;
	return res;
}
LL P(LL n,int k)
{
	if(n<=N-10) return phi[n];
	if(vis[k]) return p[k];
	LL res=1ll*(1ll+n)*n/2;
	LL l=2,r;
	for(;l<=n;l=r+1)
	{
		r=n/(n/l);
		res=res-1ll*(r-l+1)*P(n/l,k*l);
	}
	vis[k]=1;
	p[k]=res;
	return  res;
}
int main()
{
	freopen("sum.in","r",stdin);
	freopen("sum.out","w",stdout);
	init(N-10);
	int T;
	cin>>T;
	while(T--)
	{
		LL n;
		cin>>n;
		cout<<P(n,1)<<' ';
		memset(vis,0,sizeof(vis));
		cout<<S(n,1)<<endl;
		memset(vis,0,sizeof(vis));		
	}

	return 0;
}

[HDOJ5608]function

题目给出了一个函数(F(n)=n^2-3n+2) 并且他保证(F(n)=sum_{d|n}f(d)) 求(sum_{i=1}^nf(i)) 我们先求出(f(n))的递推式 莫比乌斯反演一下可得 (f(n)=sum_{d|n}F(d)mu(frac{n}{d})) 这显然是一个积性函数 但是(n)的范围很大，考虑杜教筛 发现题目给的这个(F(n))正好可以看成(f(n))和(I(n))的卷积 因此答案为

[S(n)=frac{sum_{i=1}^ni^2-3i+2-sum_{d=2}^nS(lfloor frac{n}{d}rfloor)}{1} ]

其中左半部分的三项可以分开算，其中 (1^2+2^2+3^2+……n^2=n(n+1)(2n+1)/6)

#include<bits/stdc++.h>
using namespace std;
typedef long long LL;
const int N = 1e6+7;
const LL mod = 1e9+7;
LL Pow(LL a,LL b)
{
	LL res=1;
	while(b)
	{
		if(b&1) res=1ll*res*a%mod;
		a=1ll*a*a%mod;
		b>>=1;	
	}
	return res;
}
LL inv(LL n)
{
	return Pow(n,mod-2);
}
LL inv6=Pow(3,mod-2),inv2=Pow(2,mod-2);
LL F(LL n){return (1ll*n*(n+1)/2%mod*(2ll*n+1)%mod*inv6%mod-3*(1+n)*n/2%mod+2*n+mod)%mod;}
int mu[N],prime[N],tot=0;
int v[N];
LL f[N];
void init(int n)
{
	mu[1]=1;
	f[1]=0;
	for(int i=2;i<=n;i++)
	{
		if(!v[i])
		{
			v[i]=i;
			mu[i]=-1;
			prime[++tot]=i;
		}
		for(int j=1;j<=tot;j++)
		{
			if(prime[j]>v[i]||i*prime[j]>n) break;
			v[i*prime[j]]=prime[j];
			if(i%prime[j]==0)
			{
				mu[i*prime[j]]=0;
				break;
			} 
			else mu[i*prime[j]]=mu[i]*mu[prime[j]]; 
		}
	}
	for(LL d=1;d<=n;d++)
	for(LL T=d;T<=n;T+=d)
	f[T]=(f[T]+(d*d-3*d+2)%mod*mu[T/d]+mod)%mod;
	for(int i=1;i<=n;i++)
	f[i]=(f[i-1]+f[i])%mod;
}
LL bin[N],top=0;
int tag[N];
int cnt=0;
unordered_map<LL,LL>s,vis;
LL S(LL n)
{
	if(n<=1e6) return f[n];
	if(vis[n]) return s[n];
	LL res=F(n);
	LL l=2,r;
	for(;l<=n;l=r+1)
	{
		r=n/(n/l);
		res=res-(1ll*(r-l+1)*S(n/l)%mod);
		res=(res%mod+mod)%mod;
	}
	vis[n]=1;
	s[n]=res;
	return res;
}
inline int read()
{
	int X=0; bool flag=1; char ch=getchar();
	while(ch<'0'||ch>'9') {if(ch=='-') flag=0; ch=getchar();}
	while(ch>='0'&&ch<='9') {X=(X<<1)+(X<<3)+ch-'0'; ch=getchar();}
	if(flag) return X;
	return ~(X-1);
}
void write(LL x)        
{        
    if(x < 0) {        
        putchar('-');        
        x = -x;        
    }        
    if(x > 9)        
        write(x/10);        
    putchar(x % 10 + '0');        
    return;        
}      
int main()
{
	freopen("a.in","r",stdin);
	freopen("a.out","w",stdout);
	init(1e6);
	LL T;
	cin>>T;
	LL n;
	while(T--)
	{
		n=read();
		if(n<=1e6)
		{
			write(f[n]);
			printf("n");
			continue;
		}
		write(S(n));
		printf("n");	
	}
	return 0;
}

【51nod 2026】Gcd and Lcm

题目要求

[sum_{i=1}^nsum_{j=1}^nf(gcd(i,j))f(lcm(i,j)) ]

其中(f(n)=sum_{d|n}mu(d)d) 一开始想把这个式子莫反一下，但是发现做不下去了…… 考虑分析一席(f(n))的性质 这个函数显然是积性函数，因此我们对每个质因子分别考虑 (f(p^k)=1-p)，也就是说对于每个质因子而言，他的指数是多少根本没有关系，因此一个数的(f)只跟他质因子的种类数有关 接着我们分析(gcd(i,j))和(lcm(i,j)) 1：对于(i,j)共有的因子，他们被乘了两次,设他们是A 2：对于(i,j)各自有的因子，他们都只乘了一次，假设分别是B,C 也就是 A(BCA) 我们交换一下顺序可以得到(AB)(AC) 也就是(f(gcd(i,j))times f(lcm(i,j))=f(i)times f(j)) 式子变为

[sum_{i=1}^nsum_{j=1}^nf(i)f(j) ]

[sum_{i=1}^nf(i)sum_{j=1}^nf(j) ]

[(sum_{i=1}^nf(i))^2 ]

问题变成了求(f)的前缀和 也就是

[sum_{i=1}^nsum_{d|i}mu(d)d ]

[sum_{d=1}^nmu(d)dsum_{d|i}^n1 ]

[sum_{d=1}^nmu(d)dlfloor frac{n}{d}rfloor ]

后半部分可以整除分块 问题进一步转化为求(mu(d)d)的前缀和 考虑杜教筛 我们考虑构造一个(g(n)=Id(n)=n) 也就是那么$$gtimes (mu(d)d)$$

[=sum_{d|n};mu(d)d;g(frac{n}{d}) ]

[=sum_{d|n};mu(d)d;frac{n}{d} ]

[=nsum_{d|n};mu(d) ]

[=e ]

因此前缀和就是

[S(n)=frac{1-sum_{d=2}^niS(lfloor frac{n}{d}rfloor)}{1} ]

带进去算就行了

#include<bits/stdc++.h>
using namespace std;
typedef long long LL;
const int mod = 1e9+7;
const int N =1e6+7;
int mu[N];
LL F[N];
int v[N],prime[N],tot=0;
void init(int n)
{
	mu[1]=1;
	for(int i=2;i<=n;i++)
	{
		if(!v[i])
		{
			v[i]=i;
			prime[++tot]=i;
			mu[i]=-1;
		}
		for(int j=1;j<=tot;j++)
		{
			if(v[i]<prime[j]||i*prime[j]>n) break;
			v[i*prime[j]]=prime[j];
			if(i%prime[j]==0)
			{
				mu[i*prime[j]]=0;
				break;
			}
			else mu[i*prime[j]]=mu[i]*mu[prime[j]];
		}
	}
	for(LL i=1;i<=n;i++)
	F[i]=(F[i-1]+mu[i]*i%mod)%mod;
}
unordered_map<LL,LL> s,vis;
LL S(LL n)
{
	if(n<=1e6) return F[n];
	if(vis[n]) return s[n];
	LL res=1;
	LL l=2,r;
	for(;l<=n;l=r+1)
	{
		r=n/(n/l);
		res=(res-(l+r)*(r-l+1)/2%mod*S(n/l)%mod+mod)%mod;
	}
	vis[n]=1;
	s[n]=res;
	return res;
}
LL f(LL n)
{
	LL l=1,r;
	LL ans=0;
	for(;l<=n;l=r+1)
	{
		r=(n/(n/l));
		ans=(ans+(S(r)-S(l-1)+mod)%mod*(n/l)%mod)%mod;
	}
	return ans;
}
int main()
{
	freopen("c.in","r",stdin);
	freopen("c.out","w",stdout);
	init(1e6);
	LL n;
	cin>>n;
	cout<<f(n)*f(n)%mod;
	return 0;
}

【51Nod 1238最小公倍数之和 V3

毒瘤 求

[sum_{i=1}^nsum_{j=1}^nlcm(i,j) ]

莫反的部分详见莫比乌斯反演专题(基础篇)中的【BZOJ2693】jzptab 我么这里直接用最后的的式子

[= sum_{T=1}^nS(lfloor frac{n}{T} rfloor)S(lfloor frac{m}{T} rfloor)Tsum_{k|T}^Tmu(k)k ]

其中(S(n)=sum_{i=1}^ni) 所以关键的问题还是筛出来(f(n)=nsum_{d|n}mu(d)d)的前缀和

[sum_{i=1}^nsum_{d|i}imu(d)d ]

[sum_{d=1}^nmu(d)dsum_{d|i}i ]

[sum_{d=1}^nmu(d)dsum_{i=1}^{lfloor frac{n}{d}rfloor}id ]

[sum_{d=1}^nmu(d)d^2sum_{i=1}^{lfloor frac{n}{d}rfloor}i ]

前边的部分显然我们卷上一个(g(n)=n^2)就可以杜教筛了，整个式子整除分块一下就行了 但是，别忘了外边还有一个整除分块，这样就成(O(n))的了 因此我们只能考虑别的方法 通过观察打表 等方式我们发现(ftimes (n^2) = (n)) 至于具体证明抱歉我不懂，因此我们卷上一个(g(n)=n^2) 前缀和变为

[S(n)=frac{sum_{i=1}^ni-sum_{d=2}^ni^2S(lfloor frac{n}{d}rfloor)}{1} ]

就可以解决了

#include<bits/stdc++.h>
using namespace std;
const int N = 4641600;
typedef unsigned long long LL;
int mu[N],prime[N],tot=0;
int v[N];
LL F[N];
const LL mod = 1e9+7;
void init(LL n)
{
	F[1]=1;
	for(int i=2;i<=n;i++)
	{
		if(!v[i])
		{
			v[i]=i;
			F[i]=1LL*i*((1-i+mod)%mod)%mod%mod; 
			prime[++tot]=i;
		}
		for(int j=1;j<=tot;j++)
		{
			if(prime[j]>v[i]||i*prime[j]>n) break;
			v[i*prime[j]]=prime[j];
			if(i%prime[j]==0)
			{
				F[i*prime[j]]=1LL*F[i]*prime[j]%mod;
				break;
			}
			else F[i*prime[j]]=1LL*F[i]*F[prime[j]]%mod;
		}
	}
	for(int i=1;i<=n;i++)
	F[i]=(F[i-1]+F[i])%mod;
}
unordered_map<LL,LL>s,vis;
LL Pow(LL a,LL b)
{
	LL res=1;
	while(b)
	{
		if(b&1) res=1LL*res*a%mod;
		a=1LL*a*a%mod;
		b>>=1LL;	
	}
	return res;
}
LL inv6=Pow(6,mod-2),inv2=Pow(2,mod-2);
LL sqr(LL n){n%=mod;return 1LL*n*(n+1)%mod*(2LL*n+1)%mod*inv6%mod;}
LL sum(LL n){n%=mod;return 1LL*(1+n)%mod*n%mod*inv2%mod;}
bool flag=1;
LL Sum(LL n)
{
	if(n<=N-10) return F[n];
	if(vis[n]) return s[n];
	LL res=sum(n);
	LL l=2,r;
	for(;l<=n;l=r+1)
	{
		r=n/(n/l);
		res=(res-1LL*((sqr(r)-(sqr(l-1)))+mod)%mod*Sum(n/l)%mod+mod)%mod;
	}
	vis[n]=1;
	s[n]=res;
	return res;
}
int main()
{
	freopen("d.in","r",stdin);
	freopen("d.out","w",stdout);
	init(N-10);
	LL n;
	cin>>n;
	LL ans=0;
	LL l=1,r;
	for(;l<=n;l=r+1)
	{
		r=(n/(n/l));
		ans=(ans+1LL*sum(n/l)*sum(n/l)%mod*((Sum(r)-Sum(l-1)+mod)%mod)%mod)%mod;
	}
	cout<<ans;
	return 0;
}

51【Nod 1237】最大公约数之和 V3

同上莫反部分详见莫比乌斯反演专题(基础篇)中的 [NOI2010day1]T1-能量采集 最后的式子为

[=sum_{d=1}^nvarphi(d)lfloor frac{n}{d} rfloorlfloor frac{m}{d} rfloor ]

直接杜教筛处理(varphi(n))的前缀和就行了

#include<bits/stdc++.h>
using namespace std;
typedef long long LL;
const LL mod = 1e9+7;
const LL N = 4641600;
LL phi[N];
LL v[N],prime[N],tot=0;
void init(LL n)
{
	phi[1]=1;
	for(LL i=2;i<=n;i++)
	{
		if(!v[i])
		{
			v[i]=i;
			phi[i]=i-1;
			prime[++tot]=i;
		}
		for(LL j=1;j<=tot;j++)
		{
			if(prime[j]>v[i]||i*prime[j]>n) break;
			v[i*prime[j]]=prime[j];
			if(i%prime[j]==0)
			{
				phi[i*prime[j]]=1LL*phi[i]*prime[j]%mod;
				break;
			} 
			else phi[i*prime[j]]=1LL*phi[i]*(prime[j]-1)%mod; 
		}
	}
	for(LL i=1;i<=n;i++)
	phi[i]=(phi[i-1]+phi[i])%mod;
}
LL Pow(LL a,LL b)
{
	LL res=1;
	while(b)
	{
		if(b&1) res=1LL*res*a%mod;
		a=1LL*a*a%mod;
		b>>=1; 
	}
	return res;
}
LL inv2=Pow(2ll,mod-2);
unordered_map<LL,LL>vis,s;
LL mul(LL a,LL b)
{
	LL res=0;
	while(b)
	{
		if(b&1) res=(res+a)%mod;
		a=(a+a)%mod;
		b>>=1;
	}
	return res;
}
LL Sum(LL n)
{
	if(n<=N-10) return phi[n];
	if(vis[n]) return s[n];
	LL res=mul(mul(n,n+1),inv2);
	LL l=2,r;
	for(;l<=n;l=r+1)
	{
		r=(n/(n/l));
		res=(res-1LL*(r-l+1)%mod*Sum(n/l)%mod)%mod;
	}
	res=(res+mod)%mod;
	vis[n]=1;
	s[n]=res;
	return res;
}
int main()
{
	freopen("d.in","r",stdin);
	freopen("d.out","w",stdout);
	init(N-10);
	LL n;
	cin>>n;
	LL res=0;
	LL l=1,r;
	for(;l<=n;l=r+1)
	{
		r=n/(n/l);
		res=(res+mul(Sum(r)-Sum(l-1),mul(n/l,n/l))+mod)%mod;
	}
	cout<<res;
	return 0;
}

[bzoj 4176] Lucas的数论

同上莫反部分详见莫比乌斯反演专题(基础篇)中的 [SDOI2015] 约数个数和

[ans=sum_{d=1}^nmu(d)sum_{x=1}^{lfloor frac{n}{d} rfloor}lfloor frac{n}{dx} rfloor sum_{y=1}^{lfloor frac{m}{d} rfloor}lfloor frac{m}{dy} rfloor ]

(mu)的部分可以杜教筛处理，后边的部分可以整除分块 话说这样不是变成(O(n))了吗

#include<bits/stdc++.h>
using namespace std;
typedef long long LL;
const int N = 2500005;
const int mod = 1e9+7;
int mu[N],prime[N],tot=0;
LL v[N];
void init(int n)
{
	mu[1]=1;
	for(int i=2;i<=n;i++)
	{
		if(!v[i])
		{
			v[i]=i;
			mu[i]=-1;
			prime[++tot]=i;
		}
		for(int j=1;j<=tot;j++)
		{
			if(prime[j]>v[i]||i*prime[j]>n) break;
			v[i*prime[j]]=prime[j];
			if(i%prime[j]==0)
			{
				mu[i*prime[j]]=0;
				break;
			}
			else mu[i*prime[j]]=mu[i]*mu[prime[j]];
		}
	}
	for(int i=1;i<=n;i++)
	mu[i]=(mu[i]+mu[i-1]);
}
unordered_map<LL,LL> s;
LL Sum(LL n)
{
	if(n<=N-10) return mu[n];
	if(s[n]) return s[n];
	LL res=1;
	LL l=2,r;
	for(;l<=n;l=r+1)
	{
		r=n/(n/l);
		res=(res-1ll*(r-l+1)*Sum(n/l)%mod)%mod;
	}
	res=(res+mod)%mod;
	s[n]=res;
	return res;
}
LL S(LL n)
{
	int l=1,r;
	LL res=0;
	for(;l<=n;l=r+1)
	{
		r=n/(n/l);
		res=(res+(LL)(n/l)*(r-l+1)%mod)%mod;
	}
	return res;
}
int main()
{
	freopen("math.in","r",stdin);
	freopen("math.out","w",stdout);
	init(N-10);
	LL n;
	cin>>n;
	LL res=0;
	int l=1,r;
	for(;l<=n;l=r+1)
	{
		r=n/(n/l);
		LL val=S(n/l);
		res=(res+1ll*val*val%mod*(Sum(r)-Sum(l-1)+mod)%mod)%mod;
	}
	cout<<res;
	return 0;
}

脚本宝典总结

以上是脚本宝典为你收集整理的杜教筛（基础篇）全部内容，希望文章能够帮你解决杜教筛（基础篇）所遇到的问题。

如果觉得脚本宝典网站内容还不错，欢迎将脚本宝典推荐好友。

本图文内容来源于网友网络收集整理提供，作为学习参考使用，版权属于原作者。
如您有任何意见或建议可联系处理。小编QQ：384754419，请注明来意。
标签：
猜你在找的心得技巧相关文章
其他相关热搜词更多
微信小程序PHPVueVue.jsPythonPython3SQLJSJavaScriptCSSCSS3Web前端C++C语言HTMLJavaTypeScriptes6鸿蒙OS
全站导航更多
html5HTML/XhtmlCSSXML/XSLTDreamweaver教程Frontpage教程心得技巧JavaScriptASP.NETPHP编程正则表达式AJAX相关ASP编程JSP编程编程10000问CSS/HTMLFlexvbsDOS/BAThtahtcpythonperl游戏相关VBA远程脚本ColdFusionMsSqlMysqlmariadboracleDB2mssql2008mssql2005SQLitePostgreSQLMongoDB星外虚拟主机华众虚拟主机Linuxwin服务器FTP服务器DNS服务器Tomcatnginxzabbix云和虚拟化bios系统安装系统系统进程Windows系列LINUXRedHat/CentosUbuntu/DebianFedoraSolaris红旗Linux建站经验微信营销网站优化网站策划网络赚钱网络创业站长故事alexa域名photoshop教程摄影教程Fireworks教程CorelDraw教程Illustrator教程Painter教程Freehand教程IndesignSketch笔记本主板内存CPU存储显卡显示器光存储鼠标键盘平板电脑安全教程杀毒防毒安全设置病毒查杀脚本攻防入侵防御工具使用业界动态Exploit漏洞分析
最新心得技巧教程
热门心得技巧教程