杜教筛学习笔记
杜教筛学习笔记
接着 Dirichlet 卷积，继续学习杜教筛。
本文中一些函数的定义见此博文
应用
通过杜教筛，我们可以快速求出某一数论函数 \(f\) 的前缀和，即，我们可以在低于线性的时间复杂度内求出 \(S(n)=\sum\limits_{i=1}^n f(i)\)
方法
杜教筛主要运用一个公式，通过这个公式我们建立了 \(S(n)\) 与 \(S(\frac{n}{i})\) 的关系式，已知两个数论函数 \(f,g\)，\(S(n)=\sum\limits_{i=1}^n f(i)\)，则有：
\]
我们给个证明：
\sum\limits_{i=1}^n (f\ast g)(i) \\=\sum\limits_{i=1}^n\sum\limits_{d\mid i} f(d)\cdot g(\frac{i}{d})\\=\sum\limits_{d=1}^n\sum\limits_{i=1}^{\lfloor \frac{n}{d}\rfloor} g(d)\cdot f(i)\\=\sum\limits_{d=1}^n g(d)\cdot\sum\limits_{i=1}^{\lfloor \frac{n}{d}\rfloor} f(i)\\=\sum\limits_{d=1}^n g(d)\cdot S(\lfloor \frac{n}{d}\rfloor)\\=\sum\limits_{i=1}^n g(i)\cdot S(\lfloor \frac{n}{i}\rfloor)\\&\square
\end{aligned}\]
证明的关键就在于中间对于 \(d,i\) 枚举顺序的变化，这个证明值得细品。
那么我们有了这个公式哟什么用呢？我们再来导一下：
\sum\limits_{i=1}^n (f\ast g)(i)=\sum\limits_{i=1}^n g(i)\cdot S(\lfloor\frac{n}{i}\rfloor)\\\sum\limits_{i=1}^n (f\ast g)(i)=g(1)\cdot S(n)+\sum\limits_{i=2}^n g(i)\cdot S(\lfloor\frac{n}{i}\rfloor)\\g(1)\cdot S(n)=\sum\limits_{i=1}^n (f\ast g)(i)-\sum\limits_{i=2}^n g(i)\cdot S(\lfloor\frac{n}{i}\rfloor)
\end{aligned}\]
我们可以轻松地求出 \(S(n)\) 的值。而 \(S(\lfloor\frac{n}{i}\rfloor)\) 可以通过数论分块来解决。由此看来，杜教筛的时间复杂度为低于线性复杂度。
例子
看了理论，可能还有点晕，我们来看几个例子。
求莫比乌斯函数前缀和
令 \(S(n)=\sum\limits_{i=1}^n \mu(i)\ \ \ \ (n<2^{31})\)
我们根据公式来求，我们令 \(f=\mu,g=1\)，由 Dirichlet 卷积知识，我们知道 \(\mu\ast 1=\varepsilon\)。那么：
g(1)\cdot S(n)=\sum\limits_{i=1}^n (f\ast g)(i)-\sum\limits_{i=2}^n g(i)\cdot S(\lfloor\frac{n}{i}\rfloor)\\1\cdot S(n)=\sum\limits_{i=1}^n \varepsilon(i)-\sum\limits_{i=2}^n 1\cdot S(\lfloor\frac{n}{i}\rfloor)\\S(n)=1-\sum\limits_{i=2}^n S(\lfloor\frac{n}{i}\rfloor)
\end{aligned}\]
根据数论分块的知识，\(\lfloor\frac{n}{i}\rfloor\) 的取值有 \(O(\sqrt{n})\) 个。那么时间复杂度为 \(O(\sum\limits_{i=1}2^{\frac{1}{2^i}}) \approx O(n^{\frac{3}{4}})\)
这明显会超时，我们需要进一步优化。其实可以这样想，一部分前缀和我们先预处理出来，另一部分再通过杜教筛来算。我们设预处理前 \(k\) 个，后 \(n-k\) 个用杜教筛做，那么此时时间复杂度是 \(O(k+(n-k)^{\frac{3}{4}})\)，由基本不等式得，\(k=(n-k)^{\frac{3}{4}}\) 原式取最小值。我们算一下发现 \(k=n^{\frac{2}{3}}\) 是不错的选择，此时复杂度大约为 \(O(n^{\frac{2}{3}})\)。
//Don't act like a loser.
//You can only use the code for studying or finding mistakes
//Or,you'll be punished by Sakyamuni!!!
#include<bits/stdc++.h>
#define int long long
using namespace std;
int read() {
char ch=getchar();
int f=1,x=0;
while(ch<'0'||ch>'9') {
if(ch=='-')
f=-1;
ch=getchar();
}
while(ch>='0'&&ch<='9') {
x=x*10+ch-'0';
ch=getchar();
}
return f*x;
}
const int maxn=4e7+10;
int n,mu[maxn],p[maxn/10],cnt,s[maxn];
bool is[maxn];
map<long,long> mp;//记忆化搜索
void sieve_mu(int n) {
fill(is,is+n+1,1);
mu[1]=1;
for(int i=2;i<=n;i++) {
if(is[i]) {
p[++cnt]=i;
mu[i]=-1;
}
for(int j=1;j<=cnt;j++) {
if(p[j]*i>n) {
break;
}
is[p[j]*i]=0;
if(i%p[j]==0) {
mu[i*p[j]]=0;
break;
}
mu[i*p[j]]=-mu[i];
}
}
for(int i=1;i<=n;i++){
s[i]=s[i-1]+mu[i];
}
}
int pre(int n) {
if(n<=4e7) {
return s[n];
}
if(mp[n]) {
return mp[n];
}
int j=2,ret=1;
for(int i=2;i<=n;i=j+1) {
j=(n/(n/i));//数论分块
ret-=(j-i+1)*pre(n/i);
}
return mp[n]=ret;
}
signed main() {
sieve_mu(4e7);
int t=read();
while(t--) {
printf("%lld\n",pre(read()));
}
return 0;
}
求欧拉函数前缀和
令 \(S(n)=\sum\limits_{i=1}^n \varphi(i)\ \ \ \ (n<2^{31})\)
我们根据公式来求，我们令 \(f=\varphi,g=1\)，由 Dirichlet 卷积知识，我们知道 \(\varphi\ast 1=id\)。那么：
g(1)\cdot S(n)=\sum\limits_{i=1}^n (f\ast g)(i)-\sum\limits_{i=2}^n g(i)\cdot S(\lfloor\frac{n}{i}\rfloor)\\1\cdot S(n)=\sum\limits_{i=1}^n id(i)-\sum\limits_{i=2}^n 1\cdot S(\lfloor\frac{n}{i}\rfloor)\\S(n)=\frac{n^2+n}{2}-\sum\limits_{i=2}^n S(\lfloor\frac{n}{i}\rfloor)
\end{aligned}\]
和之前一样，我们还是采取一部分直接算，一部分杜教筛的想法，时间复杂度为 \(O(n^{\frac{2}{3}})\)
//Don't act like a loser.
//You can only use the code for studying or finding mistakes
//Or,you'll be punished by Sakyamuni!!!>
#define int long long
using namespace std;
int read() {
char ch=getchar();
int f=1,x=0;
while(ch<'0'||ch>'9') {
if(ch=='-')
f=-1;
ch=getchar();
}
while(ch>='0'&&ch<='9') {
x=x*10+ch-'0';
ch=getchar();
}
return f*x;
}
const int maxn=1e7+10;
int n,phi[maxn],p[maxn/10],cnt,s[maxn];
bool is[maxn];
map<long,long> mp;
void sieve_euler(int n) {
fill(is,is+n+1,1);
phi[1]=1;
for(int i=2;i<=n;i++) {
if(is[i]) {
p[++cnt]=i;
phi[i]=i-1;
}
for(int j=1;j<=cnt;j++) {
if(p[j]*i>n) {
break;
}
is[p[j]*i]=0;
if(i%p[j]==0) {
phi[i*p[j]]=p[j]*phi[i];
break;
}
phi[i*p[j]]=(p[j]-1)*phi[i];
}
}
for(int i=1;i<=n;i++){
s[i]=s[i-1]+phi[i];
}
}
int pre(int n) {
if(n<=1e7) {
return s[n];
}
if(mp[n]) {
return mp[n];
}
int j=2,ret=(n*(n+1)/2);
for(int i=2;i<=n;i=j+1) {
j=(n/(n/i));
ret-=(j-i+1)*pre(n/i);
}
return mp[n]=ret;
}
signed main() {
sieve_euler(1e7);
int t=read();
while(t--) {
printf("%lld\n",pre(read()));
}
return 0;
}
这两个代码结合起来就可以通过洛谷上的模板题。
求莫比乌斯函数平方前缀和
令 \(S(n)=\sum\limits_{i=1}^n \mu^2(i)\ \ \ \ (n<2^{31})\)
我们继续根据公式来求，我们令 \(f=\mu^2,g=1\)，由 Dirichlet 卷积知识，我们知道 \(\mu^2*1=\mu\)。那么：
g(1)\cdot S(n)=\sum\limits_{i=1}^n (f\ast g)(i)-\sum\limits_{i=2}^n g(i)\cdot S(\lfloor\frac{n}{i}\rfloor)\\1\cdot S(n)=\sum\limits_{i=1}^n \mu(i)-\sum\limits_{i=2}^n 1\cdot S(\lfloor\frac{n}{i}\rfloor)\\S(n)=\sum\limits_{i=1}^n \mu(i)-\sum\limits_{i=2}^n S(\lfloor\frac{n}{i}\rfloor)
\end{aligned}\]
所以我们在求莫比乌斯函数平方和时，要求出莫比乌斯函数的前缀和，也是杜教筛，总时间复杂度 \(O(n^{\frac{2}{3}})\)。
//Don't act like a loser.
//You can only use the code for studying or finding mistakes
//Or,you'll be punished by Sakyamuni!!!
//#pragma GCC optimize("Ofast","-funroll-loops","-fdelete-null-pointer-checks")
//#pragma GCC target("ssse3","sse3","sse2","sse","avx2","avx")
#include<bits/stdc++.h>
#define int long long
using namespace std;
int read() {
char ch=getchar();
int f=1,x=0;
while(ch<'0'||ch>'9') {
if(ch=='-')
f=-1;
ch=getchar();
}
while(ch>='0'&&ch<='9') {
x=x*10+ch-'0';
ch=getchar();
}
return f*x;
}
const int maxn=4e6+10;
int n,mu[maxn],p[maxn],cnt,s[maxn],s2[maxn];
bool is[maxn];
map<long,long> mp,mp2;
void sieve_mu(int n) {
fill(is,is+n+1,1);
mu[1]=1;
for(int i=2;i<=n;i++) {
if(is[i]) {
cnt++;
p[cnt]=i;
mu[i]=-1;
}
for(int j=1;j<=cnt;j++) {
if(p[j]*i>n) {
break;
}
is[p[j]*i]=0;
if(i%p[j]==0) {
mu[i*p[j]]=0;
break;
}
mu[i*p[j]]=-mu[i];
}
}
s[0]=s2[0]=0;
for(int i=1;i<=n;i++){
s[i]=s[i-1]+mu[i];
s2[i]=s2[i-1]+mu[i]*mu[i];
}
}
int premu(int n) {
if(n<=maxn-10) {
return s[n];
}
if(mp[n]) {
return mp[n];
}
int j=2,ret=1;
for(int i=2;i<=n;i=j+1) {
j=(n/(n/i));
ret-=(j-i+1)*premu(n/i);
}
return mp[n]=ret;
}
int pre(int n) {
if(n<=maxn-10) {
return s2[n];
}
if(mp2[n]) {
return mp2[n];
}
int j=2,ret=premu(n);
for(int i=2;i<=n;i=j+1) {
j=(n/(n/i));
ret-=(j-i+1)*pre(n/i);
}
return mp2[n]=ret;
}
signed main() {
//freopen(".in","r",stdin);
//freopen(".out","w",stdout);
sieve_mu(maxn-10);
int t=read();
while(t--) {
printf("%lld\n",pre(read()));
}
//fclose(stdin);
//fclose(stdout);
return 0;
}
求莫比乌斯函数值平方前缀和
令 \(S(n)=\sum\limits_{i=1}^n (\mu(i))^2\ \ \ \ (n<2^{31})\)
这个题和之前的就有所不同了。我们令 \(f(i)=\mu(i)^2\)，但是要找一个方便计算的 \(g\) 才行。
我们选取函数 \(g(x)=[x=k^2\ \ \ k\in N^+]\)。我们来算一下 \(g\) 和 \(f\) 的狄利克雷卷积。我们会发现 \(f\ast g=1\)？！简要证明一下：
\((f\ast g)(n)=\sum\limits_{d\mid n} g(d)\cdot f(\frac{n}{d})\)。首先分类讨论，如果 \(d\) 是一个完全平方数，只有当 \(d\) 是 \(n\) 的所有约数中最大的完全平方数时，乘积为 \(1\)，否则 \(\mu(\frac{n}{d})=0\)。如果 \(d\) 不是完全平方数，那么 \(g(d)=0\)。得证。
那我们的计算就简化了不少。
g(1)\cdot S(n)=\sum\limits_{i=1}^n (f\ast g)(i)-\sum\limits_{i=2}^n g(i)\cdot S(\lfloor\frac{n}{i}\rfloor)\\1\cdot S(n)=\sum\limits_{i=1}^n 1-\sum\limits_{i=2}^{\sqrt{n}} 1\cdot S(\lfloor\frac{n}{i^2}\rfloor)\\S(n)=n-\sum\limits_{i=2}^{\sqrt{n}} S(\lfloor\frac{n}{i^2}\rfloor)
\end{aligned}\]
我们甚至不用数论分块就可以解决此题。采用同样的优化方法，时间复杂度还是约为 \(O(n^{\frac{2}{3}})\)
相关练习题：
