杜教筛学习笔记
杜教筛
参考资料：
https://blog.csdn.net/Ike940067893/article/details/84781307
https://www.luogu.com.cn/problemnew/solution/P4213
前置知识：
主要内容：
杜教筛
杜教筛能用低于线性的时间复杂度求积性函数前缀和。
比如要求积性函数 \(f\) 的前缀和 \(sum(n)=\sum\limits_{i=1}^nf_i\)。
需要找一个辅助积性函数 \(g\)，考虑 \((f*g)\)（\(*\) 是狄利克雷卷积，不懂自学）的前缀和 \(sum'\)：
\[\begin{split}
sum'(n)=&\sum\limits_{i=1}^n(f*g)(i)\\
=&\sum\limits_{i=1}^n\sum\limits_{d|i}f(\frac id)g(d)\\
=&\sum\limits_{d=1}^ng(d)\sum\limits_{i=1}^{\lfloor\frac nd\rfloor}f(i)\\
=&\sum\limits_{d=1}^ng(d)\color{#77cc33}{sum(\lfloor\frac nd\rfloor)}\\
\end{split}
\]
\[\begin{split}
\therefore g(1)sum(n)=&\sum\limits_{d=1}^ng(d)sum(\lfloor\frac nd\rfloor)-\sum\limits_{d=2}^ng(d)sum(\lfloor\frac nd\rfloor)\\
=&\sum\limits_{i=1}^n(f*g)(i)-\sum\limits_{d=2}^ng(d)sum(\lfloor\frac nd\rfloor)\\
\end{split}
\]
于是得到了最主要的式子，单独提出来：
\[g(1)sum(n)=\sum\limits_{i=1}^n(f*g)(i)-\sum\limits_{d=2}^ng(d)sum(\lfloor\frac nd\rfloor)
\]
如果 \(g\) 找得合适，大部分时候 \(g(1)=1\)，而 \(\sum\limits_{i=1}^n(f*g)(i)\) 和 \(\sum\limits_{i=1}^ng(i)\) 可以 \(\Theta(1)\) 算，而 \(\sum\limits_{d=2}^ng(d)sum(\lfloor\frac nd\rfloor)\) 可以整除分块算。
于是直接递归求 \(sum(n)\) 的时间复杂度为 \(\Theta(n^{\frac{3}{4}})\)。
证明
设求出 $sum(n)$ 的时间复杂度为 $\Theta(T(n))$：\[\begin{split}
T(n)=&\sqrt n+\sum\limits_{i=2}^{\sqrt n}\left(T(i)\cdot T(\lfloor\frac{n}{i}\rfloor)\right)\\
\ge&\sqrt n+\sum\limits_{i=2}^{\sqrt n}\left(\sqrt i\cdot \sqrt{\lfloor\frac{n}{i}\rfloor)}\right)\\
\ge&\sqrt n+\sum\limits_{i=2}^{\sqrt n}2\sqrt{\sqrt n}\\
=&n^{\frac{3}{4}}\\
\end{split}
\]
虽然等式步步是大于，但是如果记忆化 \(sum(n)\)，时间复杂度就是 \(\Theta(n^{\frac{3}{4}})\)。
优化：
线性筛求出前 \(n^{\frac{2}{3}}\) 个 \(sum(i)\)，然后再杜教筛剩下的，优化后的时间复杂度为 \(\Theta(n^{\frac{2}{3}})\)。
证明
设求出 $sum(n)$ 的时间复杂度为 $\Theta(T(n))$：\[\begin{split}
T(n)=&m+\sum\limits_{i=2}^{\lfloor\frac nm\rfloor}\sqrt{\lfloor\frac ni\rfloor}\\
=&m+\frac{n}{\sqrt m}\\
\ge&2n^{\frac{2}{3}}(=: \operatorname{while} m=n^{\frac{2}{3}})\\
\end{split}
\]
经典例题
【模板】杜教筛（Sum）
\(T\) 组测试数据，给定 \(n\)，求 \(\sum\limits_{i=1}^n\varphi(i)\) 和 \(\sum\limits_{i=1}^n\mu(i)\)。
数据范围：\(1\le T\le 10\)，\(1\le n\le 2^{31}-1\)。
求 \(\mu\) 的前缀和：
此时 \(f=\mu\)，如果取 \(g=1\)，那么 \(f*g=\mu*1=\epsilon\)。
同时满足：
- \(g(1)=1\)。
- \(\sum\limits_{i=1}^n(f*g)(i)=1\) 可以 \(\Theta(1)\) 算出。
- \(\sum\limits_{i=1}^ng(i)=n\) 可以 \(\Theta(1)\) 算出。
所以有：
\[sum(n)=1-\sum\limits_{d=2}^nsum(\lfloor\frac nd\rfloor)
\]
右边分块递归求解，\(sum\) 记忆化，时间复杂度 \(\Theta(n^{\frac{2}{3}})\)。
\(\texttt{Code}\)
hash<int,lng> Mu;
il lng DuMu(re int n){
if(n<=N) return mu[n];
if(Mu[n]) return Mu[n];
re lng res=1ll;
for(re int l=2,r;l<=n;l=r+1)
r=n/(n/l),res-=(r-l+1)*DuMu(n/l);
return Mu[n]=res;
}
求 \(\varphi\) 的前缀和：
此时 \(f=\varphi\)，如果取 \(g=1\)，那么 \(f*g=\varphi*1=ID\)。
同时满足：
- \(g(1)=1\)。
- \(\sum\limits_{i=1}^n(f*g)(i)=\frac{n(n+1)}{2}\) 可以 \(\Theta(1)\) 算出。
- \(\sum\limits_{i=1}^ng(i)=n\) 可以 \(\Theta(1)\) 算出。
所以有：
\[sum(n)=\frac{n(n+1)}{2}-\sum\limits_{d=2}^nsum(\lfloor\frac nd\rfloor)
\]
右边分块递归求解，\(sum\) 记忆化，时间复杂度 \(\Theta(n^{\frac{2}{3}})\)。
\(\texttt{Code}\)
hash<int,lng> Phi;
il lng DuPhi(re int n){
if(n<=N) return phi[n];
if(Phi[n]) return Phi[n];
re lng res=1ll*n*(n+1)/2;
for(re int l=2,r;l<=n;l=r+1)
r=n/(n/l),res-=(r-l+1)*DuPhi(n/l);
return Phi[n]=res;
}
$\texttt{Code}$
#include <bits/stdc++.h>
using namespace std;
//&Start
#define inf 0x3f3f3f3f
#define re register
#define il inline
#define hash unordered_map
typedef long long lng;
typedef unsigned long long ulng;
typedef vector<int> veci;
//&Data
#define N 5000000
//&Sieve
bitset<N+10> np;
int p[N+10];
lng mu[N+10],phi[N+10];
il void Sieve(){
np[1]=true;
mu[1]=1,phi[1]=1;
re int cnt=0;
for(re int i=2;i<=N;i++){
if(!np[i]) p[++cnt]=i,mu[i]=-1,phi[i]=i-1;
for(re int j=1;j<=cnt&&i*p[j]<=N;j++){
np[i*p[j]]=1;
if(i%p[j]==0){mu[i*p[j]]=0,phi[i*p[j]]=phi[i]*p[j];break;}
mu[i*p[j]]=-mu[i],phi[i*p[j]]=phi[i]*phi[p[j]];
}
}
for(re int i=2;i<=N;i++) mu[i]+=mu[i-1],phi[i]+=phi[i-1];
}
//&Du-Sieve
hash<int,lng> Mu,Phi; //记忆化sum
il lng DuPhi(re int n){
if(n<=N) return phi[n];
if(Phi[n]) return Phi[n];
re lng res=1ll*n*(n+1)/2;
for(re int l=2,r;l<=n;l=r+1)
r=n/(n/l),res-=(r-l+1)*DuPhi(n/l);
return Phi[n]=res;
}
il lng DuMu(re int n){
if(n<=N) return mu[n];
if(Mu[n]) return Mu[n];
re lng res=1ll;
for(re int l=2,r;l<=n;l=r+1)
r=n/(n/l),res-=(r-l+1)*DuMu(n/l);
return Mu[n]=res;
}
//&Main
int main(){
Sieve();
re int n,t;
scanf("%d",&t);
for(re int i=1;i<=t;i++){
scanf("%d",&n);
printf("%lld %lld\n",DuPhi(n),DuMu(n));
}
return 0;
}
毒瘤的数学题
给定 \(n\) 和 \(p\)，求
\[\left(\sum\limits_{i=1}^n\sum\limits_{j=1}^nij\gcd(i,j)\right)\bmod p \]
数据范围：\(1\le n\le \color{red}{10^{10}}\)，\(5\times10^{8}\le p\le 1.1\times 10^{9}\)。
稍微改了一下题目抬头，因为对这题的毒瘤深有体会。如果你看不到这道题的前半段推导，说明你没有学好莫比乌斯反演。
第一部分：普通の推导
\[\begin{split}
F(n)=&\sum\limits_{i=1}^n\sum\limits_{j=1}^nij\gcd(i,j)\\
=&\sum\limits_{d=1}^nd\sum\limits_{i=1}^ni\sum\limits_{j=1}^nj[\gcd(i,j)=d]\\
=&\sum\limits_{d=1}^nd\sum\limits_{i=1}^{\lfloor\frac{n}{d}\rfloor}id\sum\limits_{j=1}^{\lfloor\frac{n}{d}\rfloor}jd[\gcd(i,j)=1]\\
=&\sum\limits_{d=1}^nd\sum\limits_{i=1}^{\lfloor\frac{n}{d}\rfloor}id\sum\limits_{j=1}^{\lfloor\frac{n}{d}\rfloor}jd\sum\limits_{k|\gcd(i,j)}\mu(k)\\
=&\sum\limits_{d=1}^nd^3\sum\limits_{k=1}^{\lfloor\frac{n}{d}\rfloor}\mu(k)\sum\limits_{i=1}^{\lfloor\frac{n}{d}\rfloor}[k|i]i\sum\limits_{j=1}^{\lfloor\frac{n}{d}\rfloor}[k|j]j\\
=&\sum\limits_{d=1}^nd^3\sum\limits_{k=1}^{\lfloor\frac{n}{d}\rfloor}\mu(k)\sum\limits_{i=1}^{\lfloor\frac{n}{dk}\rfloor}ik\sum\limits_{j=1}^{\lfloor\frac{n}{dk}\rfloor}jk\\
=&\sum\limits_{d=1}^nd^3\sum\limits_{k=1}^{\lfloor\frac{n}{d}\rfloor}\mu(k)k^2\sum\limits_{i=1}^{\lfloor\frac{n}{dk}\rfloor}i\sum\limits_{j=1}^{\lfloor\frac{n}{dk}\rfloor}j\\
\end{split}
\]
设 \(S(n)=\sum\limits_{i=1}^ni=\frac{n(n+1)}{2}\)，所以
\[F(n)=\sum\limits_{d=1}^nd^3\sum\limits_{k=1}^{\lfloor\frac{n}{d}\rfloor}\mu(k)k^2S(\lfloor\frac{n}{dk}\rfloor)^2
\]
第二部分：代换
将 \(T=dk\) 带入：
\[\begin{split}
F(n)=&\sum\limits_{d=1}^nd^3\sum\limits_{k=1}^{\lfloor\frac{n}{d}\rfloor}\mu(k)k^2S(\lfloor\frac{n}{dk}\rfloor)^2\\
=&\sum\limits_{d=1}^nd^3\sum\limits_{k=1}^{\lfloor\frac{n}{d}\rfloor}\mu(k)k^2S(\lfloor\frac{n}{T}\rfloor)^2\\
=&\sum\limits_{T=1}^nS(\lfloor\frac{n}{T}\rfloor)^2\sum\limits_{d|T}d^3\mu(\frac Td)\left(\frac Td\right)^2\\
=&\sum\limits_{T=1}^nS(\lfloor\frac{n}{T}\rfloor)^2T^2\sum\limits_{d|T}d\mu(\frac Td)\\
=&\sum\limits_{T=1}^nS(\lfloor\frac{n}{T}\rfloor)^2T^2(\mu*ID)(T)\\
\end{split}
\]
将公式 \(\mu*ID=\varphi\) 带入得：
\[F(n)=\sum\limits_{T=1}^nS(\lfloor\frac{n}{T}\rfloor)^2T^2\varphi(T)
\]
第三部分：杜教卷积
这里给出其中一种方法，可能不是最优的：
将整个 \(T^2\varphi(T)\) 提出来杜教，即令 \(f(n)=n^2\varphi(n)\)。
为了防止遗忘，把最重要的套路式再拿出来遛一下：
\[g(1)sum(n)=\sum\limits_{i=1}^n(f*g)(i)-\sum\limits_{d=2}^ng(d)sum(\lfloor\frac nd\rfloor)
\]
为了抵消掉 \(f\) 中的 \(n^2\)，\(g\) 应该等于 \(n^2\)。
则有：
\[\begin{split}
(f*g)(n)=&\sum\limits_{d|n}f(d)g(\frac nd)\\
=&\sum\limits_{d|n}d^2\varphi(d)\cdot \left(\frac nd\right)^2\\
=&n^2\sum\limits_{d|n}\varphi(d)\\
\end{split}
\]
将公式 \(\varphi*1=ID\) 带入得：
\[(f*g)(n)=n^3
\]
于是便满足：
- \(g(1)=1\)。
- \(\sum\limits_{i=1}^n(f*g)(i)=S(n)^2\) 可以 \(\Theta(1)\) 算出。
- \(\sum\limits_{i=1}^ng(i)=\frac{n(n+1)(2n+1)}{6}\) 可以 \(\Theta(1)\) 算出。
所以可以杜教筛得：
\[sum(n)=\sum\limits_{i=1}^ni^2\varphi(i)=S(n)^2-\sum\limits_{d=2}^nd^2sum(\lfloor\frac{n}{d}\rfloor)
\]
然后再看原式：
\[F(n)=\sum\limits_{T=1}^nS(\lfloor\frac{n}{T}\rfloor)^2T^2\varphi(T)
\]
就可以分块加杜教解决了，复杂度 \(\Theta(n^{\frac {2}{3}})\)。
代码：要取模，要 \(\texttt{long long}\)。
$\texttt{Code}$
#include <bits/stdc++.h>
using namespace std;
//&Start
#define inf 0x3f3f3f3f
#define re register
#define il inline
#define hash unordered_map
typedef long long lng;
typedef unsigned long long ulng;
typedef vector<int> veci;
//&Data
#define N 5000000
int m;
//&Pow
il int Pow(re int a,re int x){ //用于求逆元
re int res=1;
for(;x;a=1ll*a*a%m,x>>=1)if(x&1) res=1ll*res*a%m;
return res;
}
int inv;
//&Sum
il int sum(re lng x){x%=m;return 1ll*x*(x+1)/2%m;} //即S
il int sum2(re lng x){x%=m;return 1ll*x*(x+1)%m*(2*x+1)%m*inv%m;} //即n(n+1)(2n+1)/6
//&Sieve
bitset<N+10> np;
int p[N+10],phi[N+10],f[N+10];
il void Sieve(){ //优化筛
np[1]=true,f[1]=phi[1]=1;
re int cnt=0;
for(re int i=2;i<=N;i++){
if(!np[i]) p[++cnt]=i,phi[i]=i-1;
f[i]=1ll*i*i%m*phi[i]%m;
for(re int j=1;j<=cnt&&i*p[j]<=N;j++){
np[i*p[j]]=1;
if(i%p[j]==0){phi[i*p[j]]=phi[i]*p[j];break;}
phi[i*p[j]]=phi[i]*phi[p[j]];
}
}
for(re int i=2;i<=N;i++) (f[i]+=f[i-1])%=m;
}
//&Du-Sieve0-------------->杜教筛
hash<lng,int> F;
il int DuF(re lng n){
if(n<=N) return f[n];
if(F[n]) return F[n];
re int res=sum(n); res=1ll*res*res%m;
for(re lng l=2,r;l<=n;l=r+1)
r=n/(n/l),(res-=1ll*(sum2(r)-sum2(l-1))%m*DuF(n/l)%m)%=m;
return F[n]=(res+m)%m;
}
//&Main
int main(){
re lng n;
scanf("%d%lld",&m,&n);
inv=Pow(6,m-2),Sieve();
re int ans=0,tmp;
for(re lng l=1,r;l<=n;l=r+1){
r=n/(n/l),tmp=sum(n/l);
(ans+=1ll*(DuF(r)-DuF(l-1))%m*tmp%m*tmp%m)%=m; //分块加杜教
}
printf("%d\n",(ans+m)%m);
return 0;
}
如果有感悟，会更新。
祝大家学习愉快！
