【学习笔记】杜教筛

补充几个前置知识

莫比乌斯反演
$$\phi \ast 1=id$$

杜教筛

杜教筛能够在 $\Theta (n^{\frac{2}{3}})$求出一个积性函数的前缀和。

给出一个积性函数$f(i)$，求$$\sum _{i=1}^{n}f(i)$$

构造一个前缀和容易求的积性函数$g$，令$h=f\ast g$（其中$\ast$是狄利克雷卷积），使得$h$的前缀和也容易求，那么

$$\begin{gathered} \sum _{i=1}^{n}h(i) \\ =\sum _{i=1}^n\sum _{d|i}g(d)f(\frac{i}{d}) \\ =\sum _{d=1}^{n}g(d)\sum _{i=1}^{\lfloor \frac{n}{d}\rfloor }f(i) \\ =\sum _{d=1}^{n}g(d)S(\lfloor \frac{n}{d}\rfloor ) \\ =g(1)S(n)+\sum _{d=2}^{n}g(d)S(\lfloor \frac{n}{d}\rfloor ) \end{gathered}$$

所以

$$g(1)S(n)=\sum _{i=1}^{n}h(i)-\sum _{d=2}^{n}g(d)S(\lfloor \frac{n}{d}\rfloor )$$
不难发现左边部分就是$h(i)$的前缀和，右边可以通过整除分块做。
实现的时候需要先预处理前$\Theta (n^{\frac{2}{3}})$项，然后用map存下求出的已知项。时间复杂度$\Theta (n^{\frac{2}{3}})$（不会证qwq）

int S(int n)
{
	if(n<=5000000) return sum[n];
	if(mp[n]!=0) return mp[n];
	int ans=g(n);
	for(int l=2,r; l<=n; l=r+1)
	{
		r=n/(n/l);
		(ans-=S(n/l)*(t(r)-t(l-1))%mod)%=mod;
	}
	return mp[n]=ans;
}
1
2
3
4
5
6
7
8
9
10
11
12

如何构造g？

我们的目的是让$h=f\ast g$的前缀和也容易求，所以可以先把$(f\ast g)(n)$写出来，看看$h(n)$长啥样，努力凑出一个可以求前缀和的形式。

例题

【模板】杜教筛

题目描述

给定一个正整数，求

$$an{s}_1=\sum _{i=1}^n\phi (i)$$

$$an{s}_2=\sum _{i=1}^n\mu (i)$$

输入格式

本题单测试点内有多组数据。

输入的第一行为一个整数，表示数据组数 $T$。

接下来 $T$ 行，每行一个整数 $n$，表示一组询问。

输出格式

对于每组询问，输出一行两个整数，分别代表 $an{s}_1$ 和 $an{s}_2$。

样例 #1

样例输入 #1

6
1
2
8
13
30
2333
1
2
3
4
5
6
7

样例输出 #1

1 1
2 0
22 -2
58 -3
278 -3
1655470 2
1
2
3
4
5
6

提示

数据规模与约定

对于全部的测试点，保证 $1\le T\le 10$，$1\le n<2^{31}$。

题解

莫比乌斯函数前缀和

$$\mu \ast 1=\epsilon$$
不难想到令$g=1$，那么$h=\epsilon$
所以式子就变成了
$$1\times S(n)=1-\sum _{d=2}^{n}S(\lfloor \frac{n}{d}\rfloor )$$

欧拉函数前缀和

方法一：根据定义用莫比乌斯反演做。
$$\sum _{i=1}^n\phi (i)=\sum _{i=1}^n\sum _{j=i+1}^n[gcd(i,j)=1]$$
这不就是莫比乌斯反演板子题吗

方法二：杜教筛
$$\phi \ast 1=id$$
令$g=1$得
$$S(n)=\sum _{i=1}^{n}i-\sum _{i=2}^{n}S(\lfloor \frac{n}{d}\rfloor )$$

代码

#include
#define int long long
using namespace std;
const int N=1e6+7;
int mu[N],phi[N],s1[N],s2[N];
bool bz[N];
vector<int> p;
map<int,int> mp1,mp2;
void init(int n)
{
	mu[1]=1; phi[1]=1;
	for(int i=2; i<=n; i++)
	{
		if(!bz[i])
		{
			p.push_back(i);
			mu[i]=-1; phi[i]=i-1;
		}
		for(int j:p)
		{
			if(i*j>n) break;
			bz[i*j]=1;
			if(i%j==0)
			{
				phi[i*j]=phi[i]*j;
				break;
			}
			mu[i*j]=-mu[i];
			phi[i*j]=phi[i]*(j-1);
		}
	}
	for(int i=1; i<=n; i++)
	{
		s1[i]=mu[i]+s1[i-1];
		s2[i]=phi[i]+s2[i-1];
	}
}
int S1(int n)
{
	if(n<2) return n;
	if(n<=N-7) return s1[n];
	if(mp1[n]) return mp1[n];
	int ans=1;
	for(int l=2,r; l<=n; l=r+1)
	{
		r=n/(n/l);
		ans-=S1(n/l)*(r-l+1);
	}
	return mp1[n]=ans;
}
int S2(int n)
{
	if(n<2) return n;
	if(n<=N-7) return s2[n];
	if(mp2[n]) return mp2[n];
	int ans=n*(n+1)/2;
	for(int l=2,r; l<=n; l=r+1)
	{
		r=n/(n/l);
		ans-=S2(n/l)*(r-l+1);
	}
	return mp2[n]=ans;
}
void O_o()
{
	int n;
	cin>>n;
	cout<<S2(n)<<" "<<S1(n)<<"\n";
}
signed main()
{
	ios::sync_with_stdio(false);
	cin.tie(0); cout.tie(0);	
	int T=1;
	init(N-7);
	cin>>T;
	while(T--)
	{
		O_o();
	}
}
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81

简单的数学题

题目描述

由于出题人懒得写背景了，题目还是简单一点好。

输入一个整数 $n$ 和一个整数 $p$，你需要求出：

$$\left(\sum _{i=1}^n\sum _{j=1}^{n}ij\gcd (i,j)\right)\mathrm{mod}p$$

其中 $\gcd (a,b)$ 表示 $a$ 与 $b$ 的最大公约数。

输入格式

一行两个整数 $p,n$。

输出格式

一行一个整数表示答案。

样例 #1

样例输入 #1

998244353 2000
1

样例输出 #1

883968974
1

提示

对于20%的数据，$n\le 1000$。

对于30%的数据，$n\le 5000$。

对于60%的数据，$n\le 10^6$，时限1s。

对于另外20%的数据，$n\le 10^9$，时限3s。

对于最后20%的数据，$n\le 10^{10}$，时限4s。

对于100%的数据，$5\times 10^8\le p\le 1.1\times 10^9$ 且 $p$ 为质数。

题解

引入欧拉反演（用莫比乌斯反演能做，推出来是一样的，但是很麻烦）

$$\begin{gathered} \sum _{i=1}^n\sum _{j=1}^{n}gcd(i,j)ij \\ =\sum _{i=1}^n\sum _{j=1}^{n}id(gcd(i,j))ij \\ =\sum _{i=1}^n\sum _{j=1}^n\sum _{d|gcd(i,j)}\phi (d)ij \\ =\sum _{d=1}^n\phi (d)\ast d^2\sum _{i=1}^{\lfloor \frac{n}{d}\rfloor }\sum _{j=1}^{\lfloor \frac{n}{d}\rfloor }ij \end{gathered}$$

令$sum(n)={\sum }_{i=1}^n{\sum }_{j=1}^{n}ij$
显然$$sum(n)=\frac{n^2(n+1{)}^2}{4}$$

令$f(i)=\phi (i)\ast i^2$，现在我们如果能求出它的前缀和，我们就可以整除分块解决问题。这个时候就要杜教筛了：
$$\begin{gathered} (f\ast g)(n) \\ =\sum _{d|n}\phi (d)\ast d^2\ast g(\frac{n}{d}) \end{gathered}$$
$d^2$特别烦人怎么办？消掉！
令$g(i)=i^2$
$$\begin{gathered} (f\ast g)(n) \\ =\sum _{d|n}\phi (d)\ast d^2\ast g(\frac{n}{d}) \\ =n^2\sum _{d|n}\phi (d) \\ =n^3 \end{gathered}$$

用我们小学二年级就学过的立方和公式得
$$\sum _{i=1}^n{i}^3=\frac{n^2(n+1{)}^2}{4}=sum(n)$$
于是乎，整除分块即可

代码

#include
#define int long long
using namespace std;
const int N=5e6+7;
int mod,inv6,phi[N],sum[N];
bool bz[N];
vector<int> p;
map<int,int> mp;
int power(int x,int t)
{
	int b=1;
	while(t)
	{
		if(t&1) b=b*x%mod;
		x=x*x%mod; t>>=1;
	}
	return b;
}
void init(int n)
{
	phi[1]=1;
	for(int i=2; i<=n; i++)
	{
		if(!bz[i])
		{
			p.push_back(i); phi[i]=i-1;
		}
		for(int j:p)
		{
			if(i*j>n) break;
			bz[i*j]=1;
			if(i%j==0)
			{
				phi[i*j]=phi[i]*j%mod;
				break;
			}
			phi[i*j]=phi[i]*phi[j]%mod;
		}
	}
	for(int i=1; i<=n; i++) sum[i]=(sum[i-1]+phi[i]*i%mod*i%mod)%mod;
}
int t(int n)
{
	n%=mod;
	return n*(n+1)%mod*(2*n%mod+1)%mod*inv6%mod;
}
int g(int n)
{
	n%=mod;
	return (n*(n+1)/2%mod)*(n*(n+1)/2%mod)%mod;
}
int S(int n)
{
	if(n<=5000000) return sum[n];
	if(mp[n]!=0) return mp[n];
	int ans=g(n);
	for(int l=2,r; l<=n; l=r+1)
	{
		r=n/(n/l);
		(ans-=S(n/l)*(t(r)-t(l-1))%mod)%=mod;
	}
	return mp[n]=ans;
}
void O_o()
{
	int n,ans=0;
	cin>>mod>>n;
	inv6=power(6,mod-2);
	init(5000000);
	for(int l=1,r; l<=n; l=r+1)
	{
		r=n/(n/l);
		(ans+=(S(r)-S(l-1))*g(n/l)%mod)%=mod;
	}
	cout<<(ans+mod)%mod;
}
signed main()
{
	ios::sync_with_stdio(false); cin.tie(0); cout.tie(0);
	cout<<fixed<<setprecision(2);
	int T=1;
	
//	cin>>T;
	while(T--)
	{
		O_o();
	}
}
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88