您也可以阅读体验更好、叙述更严谨的论文。
杜教筛是一种可以用低于线性的复杂度计算出积性函数前缀和的神奇算法。
在开始前，请您先阅读有关积性函数和莫比乌斯反演的文章，确保您已经理解了。
我们现在要求积性函\(f\)的前缀和，令我们要求的值为\(S(n)=\sum_{i=1}^{n}f(i)\)。我们找到另外两个积性函数\(g\)和\(h\)，使得\(h=f*g\)。
现在，我们写出\(h\)的前缀和的式子：
\[
\sum_{i=1}^{n}h(i)&=\sum_{i=1}^{n}(f*g)(i)\\
&=\sum_{i=1}^{n}\sum_{d|i}f(\frac{i}{d})\cdot g(d)\text{(展开狄利克雷卷积)}\\
&=\sum_{d=1}^{n}g(d)\cdot\sum_{i=1}^{\lfloor\frac{n}{d}\rfloor}f(i)\text{(交换求和号，用i替换\lfloor\frac{i}{d}\rfloor并直接枚举)}\\
&=\sum_{d=1}^{n}g(d)\cdot S(\lfloor\frac{n}{d}\rfloor)\\
&=g(1)\cdot S(n)+\sum_{d=2}^{n}g(d)\cdot S(\lfloor\frac{n}{d}\rfloor)\text{(把第一项提出来)}
\]
然后把\(g(1)\cdot S(n)\)扔到等号一边，得到：
\[g(1)\cdot S(n)=\sum_{i=1}^{n}h(i)-\sum_{d=2}^{n}g(d)\cdot S(\lfloor\frac{n}{d}\rfloor)\]
其中\(g(1)\)是\(1\)，\(\sum_{i=1}^{n}h(i)\)一般用一个\(O(1)\)的式子计算，至于右边的求和可以整除分块一下，其中的\(S\)要递归计算，同时记忆化一下。记忆化时要哈希或是用C++11的unordered_map，尽量别用带一个\(\log\)的map。
这样做的复杂度，积分一下可以得到是\(O(n^{\frac{3}{4}})\)的，如果用线筛筛到\(n^{\frac{2}{3}}\)之后再开始杜教筛，复杂度是\(O(n^{\frac{2}{3}})\)。（实践证明，杜教筛的常数略大，所以一般会多筛一些）
以上就是杜教筛的通式，对于一个要求的函数我们只要找到一个前缀和好求的\(h\)和对应的适合整除分块的\(g\)，然后套进去就好了。对于莫比乌斯函数\(\mu\)，我们很容易想到一个\(\mu*I=\epsilon\)；欧拉函数\(\varphi\)也有\(\varphi*I=id\)。于是我们就可以得到：
\[S_{\mu}(n)=1-\sum_{i=2}^{n}S_{\mu}(\lfloor\frac{n}{i}\rfloor)\]\[S_{\varphi}(n)=\frac{n(n+1)}{2}-\sum_{i=2}^{n}S_{\varphi}(\lfloor\frac{n}{i}\rfloor)\]
现在考虑一个不那么简单的例子：\(f(n)=n\times\varphi(n)\)，我们考虑把它和一个假想的\(g\)卷上并展开：\[h(n)=(f*g)(n)=\sum_{d|n}d\cdot\varphi(d)\cdot g(\frac{n}{d})\]
我们考虑配一个\(g=id\)把多出来的这个\(d\)消掉：\[h(n)=\sum_{d|n}n\cdot\varphi(d)=n\cdot\sum_{d|n}\varphi(d)=n^2\]
我们就有答案了：\[S_{f}(n)=\frac{n(n+1)(2n+1)}{6}-\sum_{i=2}^{n}i\cdot S_{f}(\lfloor\frac{n}{i}\rfloor)\]
对于后面的整除分块，我们套用等差数列求和公式即可。
以下是洛谷模板P4213 【模板】杜教筛的代码，不得不说洛谷的数据真是弱，连2147483647都没有
代码
#include <cstdlib>
#include <iostream>
#include <cstdio>
#include <algorithm>
#include <unordered_map>
#include <climits>
using namespace std;
template<typename IntType>
void read(IntType& val) {
val = 0;
int c;
bool inv = false;
while (!isdigit(c = getchar()))
if (c == '-')
inv = true;
do {
val = (val << 1) + (val << 3) + c - '0';
} while (isdigit(c = getchar()));
if (inv)
val = -val;
}
typedef long long ll;
typedef unsigned long long ull;
const int MaxN=3200000+10;
//const int MaxN=1664510+10;
int linearx=3200000;
//int linearx=1664510;
bool flag[MaxN];
int primes[MaxN],pcnt;
ll phi[MaxN],mu[MaxN];
ll sp[MaxN],sm[MaxN];
void linearsieve(){
for(int i=2;i<=linearx;i++)
flag[i]=true;
flag[1]=false;
phi[1]=mu[1]=1;
for(int i=2;i<=linearx;i++){
if(flag[i]){
primes[++pcnt]=i;
phi[i]=i-1;
mu[i]=-1;
}
for(int j=1;j<=pcnt;j++){
int p=primes[j];
if(p*i>linearx)
break;
flag[p*i]=false;
if(i%p==0){
phi[p*i]=p*phi[i];
mu[p*i]=0;
break;
}else{
phi[p*i]=phi[p]*phi[i]; // =phi[i]*(p-1)
mu[p*i]=-mu[i];
}
}
}
for(int i=1;i<=linearx;i++){
sp[i]=sp[i-1]+phi[i];
sm[i]=sm[i-1]+mu[i];
}
}
unordered_map<int,ll> phids,muds;
ll dusievephi(int x){
if(x<=linearx)
return sp[x];
auto it=phids.find(x);
if(it!=phids.end())
return it->second;
ull ans=(ull)x*((ull)x+1)/2ll;
int l,r;
for(int l=2;l<=x;l=r+1){
r=x/(x/l);
ans-=(r-l+1)*(dusievephi(x/l));
}
return phids[x]=ans;
}
ll dusievemu(int x){
if(x<=linearx)
return sm[x];
auto it=muds.find(x);
if(it!=muds.end())
return it->second;
ll ans=1;
int l,r;
for(int l=2;l<=x;l=r+1){
r=x/(x/l);
ans-=(r-l+1)*(dusievemu(x/l));
}
return muds[x]=ans;
}
int main(int argc, char* argv[]) {
// linearx=10;
linearsieve();
int t;
read(t);
int a;
for(int i=1;i<=t;i++){
read(a);
printf("%lld %lld\n",dusievephi(a),dusievemu(a));
}
return 0;
}
LaTeX 好像炸了欸