「算法笔记」杜教筛
upd on 2022.1.13：小修一波，并加了复杂度证明。
一、前置概念
具体在 「算法笔记」莫比乌斯反演 写过，所以「前置概念」就简单写写。积性函数和完全积性函数就不写了。
狄利克雷卷积：对于两个数论函数 \(f,g\)，定义它们的狄利克雷卷积 \(h=f*g\) 为
狄利克雷卷积满足交换律、结合律、对加法的分配律，有单位元 \(\varepsilon\)（其中 \(\varepsilon\) 为单位函数 \(\varepsilon(x)=[x=1]\)）。若 \(f,g\) 是积性函数，则 \(f*g\) 也是积性函数。
在狄利克雷卷积的意义下，\(\mu*1=\varepsilon\)，即 \(\mu\) 和 \(1\) 互为逆元。
常用卷积：
- \(\varepsilon=\mu*1\)，即 \(\varepsilon(n)=\sum_{d\mid n}\mu(d)\)。
- \(\varphi*1=\text{Id}\)，即 \(\sum_{d\mid n}\varphi(d)=n\)。两边同时卷一个 \(\mu\) 可得 \(\varphi=\text{Id}*\mu\)，即 \(\varphi(n)=\sum_{d\mid n}d\cdot \mu(\frac n d)\)。
- \(d=1*1 \Leftrightarrow d(n)=\sum_{d\mid n}1\)，\(\sigma_k=\text{Id}_k*1 \Leftrightarrow \sigma_k(n)=\sum_{d|n} d^k\)
提共因式：对于任意数论函数 \(f,g\) 和 完全积性函数 \(h\) 均有 \((f*g)\cdot h=(f\cdot h)*(g\cdot h)\)（注意中间是点乘，\((f\cdot g)(i)=f(i)g(i)\)），展开即可证。比如 \((\varphi\cdot \text{Id}_k)* (1\cdot \text{Id}_k)=(\varphi*1)\cdot \text{Id}_k=\text{Id}\cdot \text{Id}_k=\text{Id}_{k+1}\)。
二、杜教筛
1. 基本思想
对于数论函数 \(f\)，求 \(S(n)=\sum_{i=1}^nf(i)\)。
假如我们能快速计算 \(g\) 和 \(f*g\) 的前缀和，考虑怎么求 \(S\)。
则可以得到 \(S(n)\) 关于 \(S(\lfloor \frac{n}{i}\rfloor)\) 的递推式：
如果我们可以快速求出 \(\sum_{i=1}^n(f*g)(i)\)，再用整除分块来求 \(\sum_{i=2}^ng(i)S\left(\left\lfloor\frac{n}{i}\right\rfloor\right)\)，就能求得 \(g(1)S(n)\)。也就是说，我们要找到一个合适的 \(g\)，使得 \(f*g\) 和 \(g\) 都能快速求前缀和。
\(S(\lfloor\frac{n}{i}\rfloor)\) 是一个子问题，对 \(S(n)\) 递归求解并记忆化即可。
线性筛预处理前 \(n^{\frac{2}{3}}\) 个 \(f(x)\) 的值，则时间复杂度为 \(\mathcal O(n^{\frac{2}{3}})\)。
（当然如果我们能快速求出 \(f,g\) 的前缀和，也可以求 \(f*g\) 的前缀和，根据最初那个式子整除分块右边推左边就好了）
2. 复杂度证明
一种不用积分的证法（不保证我没写错 qwq）。
先不考虑预处理。首先我们每次记忆化搜索用到的 \(S\) 可以写成 \(S(\lfloor\large\frac n k\normalsize\rfloor)\) 的形式，再加上每次整除分块用的 \(\sqrt x\)，总复杂度就是 \(\sum_{\large x,\exists k, \lfloor\frac n k\rfloor=x}\sqrt x\)。
-
Lemma：\(\lfloor\sqrt x\rfloor=p\)，那么 \(x\in[p^2,(p+1)^2)\)，即 \(x\) 有 \((p+1)^2-p^2=\mathcal O(p)\) 个不同的取值（下面的根号为向下取整）。
-
若 \(x=\lfloor\large\frac n k\normalsize\rfloor\leq \sqrt n\)，对和式的贡献为 \(\sum\limits_{x=1}^{\sqrt x}\sqrt x\)，枚举 \(i=\sqrt x\)，根据 Lemma 有 \(\sum\limits_{i=1}^{n^{\frac 1 4}}i^2\)（有 \(\mathcal O(i)\) 个 \(\sqrt x=i\) 的 \(x\)，每个贡献 \(i\)），平方求和式个三次方的东西所以有 \(\mathcal O((n^\frac 1 4)^3)=\mathcal O(n^{\frac 3 4})\)。
-
若 \(x=\lfloor\large\frac n k\normalsize\rfloor>\sqrt n\)，那么对应的 \(k\leq \sqrt n\)，对和式的贡献为 \(\sum\limits_{k=1}^{\sqrt n}\sqrt{\large\frac n k\normalsize }=\sqrt n\sum\limits_{k=1}^{\sqrt n}\frac{1}{\sqrt k}\)，枚举 \(i=\sqrt k\)，根据 Lemma 有 \(\sqrt n\sum\limits_{i=1}^{n^{\frac 1 4}}\frac{i}{i}=n^{\frac 3 4}\)。
因此总复杂度也是 \(\mathcal O(n^{\frac 3 4})\) 级别的。还可以进一步优化杜教筛，即先线性筛出前 \(m\) 个答案，之后再用杜教筛。假设 \(m>\sqrt n\)，优化后复杂度变为 \(\mathcal O(m+\sum\limits_{k=1}^{\large\lfloor\frac n m\rfloor}\sqrt{\large\frac n k})=\mathcal O(m+\sqrt n\cdot \sqrt{\large\frac n m})\)，根据均值不等式当 \(m=\sqrt n\cdot \sqrt{\large\frac n m}\) 时复杂度最优，此时 \(m=n^{\frac 2 3}\)，时间复杂度达到最小值 \(\mathcal O(n^{\frac 2 3})\)。
实测 1s 杜教筛能跑 \(n=10^{10}\)。
upd on 2022.1.14：
pwj 指正：并不能按 \(m=\sqrt n\cdot \sqrt{\large\frac n m}=\large\frac{n}{\sqrt m}\) 算，因为 \(m+\large\frac{n}{\sqrt m}\normalsize\geq 2\sqrt{m\cdot \large\frac{n}{\sqrt m}}\) 会发现右边的 \(m\) 消不掉，而 \(m\) 是我们自己定的。但是数量级没错。
这里有个拆项凑系数用来证不等式的 trick：
\[m+\frac{n}{\sqrt m}=m+\frac{\frac n 2}{\sqrt m}+\frac{\frac n 2}{\sqrt m}\geq 3\sqrt[3]{m\cdot \frac{\frac n 2}{\sqrt m}\cdot \frac{\frac n 2}{\sqrt m}}=3\sqrt[3]{\frac{n^2}{4}}=\mathcal O(n^{\frac 2 3}) \]当 \(m=\large\frac{\frac n 2}{\sqrt m}\) 时复杂度最优。
3. 例子
求 \(S(n)=\sum_{i=1}^n\varphi(i)\)。
\(\varphi*1=\text{Id}\)。显然 \(1\)、\(\text{Id}(x)=x\) 都可以快速求前缀和。取 \(g(x)=1\) 即可。
\(S(n)=\sum_{i=1}^n i-\sum_{i=2}^nS\left(\left\lfloor\frac{n}{i}\right\rfloor\right)=\frac{n(n+1)}{2}-\sum_{i=2}^nS\left(\left\lfloor\frac{n}{i}\right\rfloor\right)\)。
求 \(S(n)=\sum_{i=1}^n\mu(i)\)。
\(\mu*1=\varepsilon\)。显然 \(1\) 和 \(\varepsilon(x)=[x=1]\) 可以快速求前缀和。取 \(g(x)=1\) 即可。
\(S(n)=1-\sum_{i=2}^nS\left(\left\lfloor\frac{n}{i}\right\rfloor\right)\)。
//Luogu P4213 #include<bits/stdc++.h> #define int long long using namespace std; const int N=2e6+5; int t,n=2e6,cnt,p[N],phi[N],mu[N]; bool vis[N]; map<int,int>Smu,Sphi; int S_mu(int n){ if(n<=2e6) return mu[n]; if(Smu[n]) return Smu[n]; int ans=0; for(int l=2,r=0;l<=n;l=r+1) //注意从 2 开始 r=n/(n/l),ans+=(r-l+1)*S_mu(n/l); return Smu[n]=1-ans; } int S_phi(int n){ if(n<=2e6) return phi[n]; if(Sphi[n]) return Sphi[n]; int ans=0; for(int l=2,r=0;l<=n;l=r+1) r=n/(n/l),ans+=(r-l+1)*S_phi(n/l); return Sphi[n]=n*(n+1)/2-ans; } signed main(){ scanf("%lld",&t); vis[0]=vis[1]=1,phi[1]=mu[1]=1; for(int i=2;i<=n;i++){ if(!vis[i]) p[++cnt]=i,phi[i]=i-1,mu[i]=-1; for(int j=1;j<=cnt&&i*p[j]<=n;j++){ vis[i*p[j]]=1; if(i%p[j]==0){phi[i*p[j]]=phi[i]*p[j],mu[i*p[j]]=0;break;} phi[i*p[j]]=phi[i]*phi[p[j]],mu[i*p[j]]=-mu[i]; } } for(int i=1;i<=n;i++) mu[i]+=mu[i-1],phi[i]+=phi[i-1]; while(t--){ scanf("%lld",&n); printf("%lld %lld\n",S_phi(n),S_mu(n)); } return 0; }
再比如说求 \(S(n)=\sum_{i=1}^n i\cdot \varphi(i)\)。
很明显 \(f(x)=x\cdot\varphi(x)\) 是积性函数。我们尝试构造出合适的 \(g(x)\)。
\[(f*g)(x)=\sum_{d\mid x}d\cdot \varphi(d)\cdot g(\frac x d) \]考虑将 \(d\) 消掉，这样可以利用 \(\sum_{d\mid x}\varphi(d)=x\) 推出一些东西。尝试 \(g=\text{Id}\)。
\[(f*g)(x)=\sum_{d\mid x}d\cdot \varphi(d)\cdot \frac{x}{d}=x\sum_{d\mid x}\varphi(d)=x^2 \]取 \(g=\text{Id}\) 即可。
求 \(S(n)=\sum_{i=1}^n\varphi(i)\cdot i^k\)。
按照提公因式卷上一个 \(1\cdot \text{Id}_k\)，那么 \((\varphi\cdot \text{Id}_k)*(1\cdot\text{Id}_k)=(\varphi*1)\cdot \text{Id}_k=\text{Id}\cdot\text{Id}_k=\text{Id}_{k+1}\)。
取 \(g=\text{Id}_k\) 即可。
求 \(S(n)=\sum_{i=1}^n \mu(i)\cdot i^k\)。
按照提公因式卷上一个 \(1\cdot \text{Id}_k\)，那么 \((\mu\cdot \text{Id}_k)*(1\cdot \text{Id}_k)=(\mu*1)\cdot \text{Id}_k=\varepsilon\cdot \text{Id}_k=\varepsilon\)。取 \(g=\text{Id}_k\) 即可。
三、狄利克雷生成函数
有些时候不容易看出 \(g(x)\) 取什么，比如：求 \(S(n)=\sum_{i=1}^n i\cdot \varphi(i)\)。凑是一种方法；而使用 狄利克雷生成函数（DGF） 可以从另一个角度求出 \(g(x)\)。
对于无穷序列 \(f_1,f_2,\cdots\)，定义其狄利克雷生成函数为
与普通生成函数以及指数生成函数的对比：
普通生成函数：\(F(x)=\sum_{i\geq 0}f_ix^i\)。
指数生成函数：\(\hat F(x)=\sum_{i\geq 0}f_i\frac{x^i}{i!}\)。
对于两个序列 \(f,g\)，其 DGF 的乘积对应 \(f,g\) 的狄利克雷卷积序列：
如果序列 \(f\) 满足积性（积性函数），其 DGF 可以由质数幂处的取值表示：
可以考虑 \(i=\prod p_j^{e_j}\)，\(f(i)=\prod f(p_j^{e_j})\)。那么有：
\[\frac{f_i}{i^x}=\prod \frac{f(p_j^{e_j})}{(p_j^{e_j})^x} \]
咕咕咕……
四、例题
Problem：给定 \(n,p\)，求：
\[\left(\sum_{i=1}^n\sum_{j=1}^n i\cdot j\cdot \gcd(i,j)\right)\bmod p \]\(n\leq 10^{10},5\times 10^8\leq p\leq 1.1\times 10^9\)，\(p\) 是质数。
Solution：
可以用 \(\varphi*1=\text{Id}\) 推式子，即 \(\sum_{d\mid x}\varphi(d)=x\)。
设 \(s(n)=(\sum_{i=1}^n i)^2=(\frac{(1+n)n}{2})^2\)，则 \(ans=\sum_{d=1}^n\varphi(d)d^2s(\lfloor\frac{n}{d}\rfloor)\)，如果能快速求出 \(f(x)=\varphi(x)x^2\) 的前缀和，那么整除分块就做完了。
考虑用杜教筛来筛 \(f(x)=\varphi(x)x^2\)，要找到一个合适的 \(g\)，使得 \(f*g\) 和 \(g\) 都可以快速求前缀和。
考虑把这里的 \(d^2\) 消掉，我们尝试 \(g=\text{Id}_2\)（幂函数 \(\text{Id}_k(n)=n^k\)），那么：
取 \(g=\text{Id}_2\) 即可，\(f*g=\text{Id}_3\)。当然也可以用 狄利克雷生成函数（DGF） 推出。
#include<bits/stdc++.h> #define int long long using namespace std; const int N=2e6+5; int n,m=2e6,mod,cnt,p[N],phi[N],inv2,inv6,ans; bool vis[N]; map<int,int>s; int mul(int x,int n,int mod){ int ans=mod!=1; for(x%=mod;n;n>>=1,x=x*x%mod) if(n&1) ans=ans*x%mod; return ans; } int s2(int n){ //平方和 n%=mod; return n*(n+1)%mod*(2*n+1)%mod*inv6%mod; } int s3(int n){ //立方和 n%=mod; int x=n%mod*(n+1)%mod*inv2%mod; return x*x%mod; } int S(int n){ if(n<=m) return phi[n]; if(s[n]) return s[n]; int ans=0; for(int l=2,r=0;l<=n;l=r+1) r=n/(n/l),ans=(ans+(s2(r)-s2(l-1))%mod*S(n/l)%mod)%mod; return s[n]=((s3(n)-ans)%mod+mod)%mod; } signed main(){ scanf("%lld%lld",&mod,&n); vis[0]=vis[1]=1,phi[1]=1; for(int i=2;i<=m;i++){ if(!vis[i]) p[++cnt]=i,phi[i]=i-1; for(int j=1;j<=cnt&&i*p[j]<=m;j++){ vis[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(int i=1;i<=m;i++) phi[i]=(phi[i-1]+i*i%mod*phi[i]%mod)%mod; inv2=mul(2,mod-2,mod),inv6=mul(6,mod-2,mod); for(int l=1,r=0;l<=n;l=r+1) r=n/(n/l),ans=(ans+s3(n/l)*(S(r)-S(l-1))%mod)%mod; printf("%lld\n",(ans+mod)%mod); return 0; }
