• 欢迎光临~
> > 杜教筛（基础篇）

杜教筛（基础篇）

开发技术 开发技术 2022-01-28 139次浏览

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

杜教筛

设(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;
}
程序员灯塔
转载请注明原文链接：杜教筛（基础篇）
喜欢 (0)
违法和不良信息举报电话：022-22558618 举报邮箱：dljd@tidljd.com