杜教筛瞎推 学习笔记【杜教筛】【数学】
点击量：251
〇、前言
对于 bzoj3944 来说，和莫比乌斯反演等其他知识关系不大，但是 $\mu$ 函数在自变量较大情况下的前缀和在反演题中也是会被用到的。
接下来通过 bzoj3944 Sum 一题引入杜教筛的思想。
参考资料
一、例题
给定一个正整数 $N(0\le N\le 2^{31}-1)$，求 $\sum_{i=1}^n\varphi(i)$ 和 $\sum_{i=1}^n\mu(i)$。多组询问。
1.分析
对于多组询问求积性函数前缀和，我们可以 $O(n)$ 线筛预处理，$O(1)$ 询问。其中空间复杂度当然也是 $O(n)$ 的，所以我们无法直接线筛，但是之后的部分还是会用到的。
2.推导
其实地上本没有路，走的人多了，也便成了路。—— 鲁迅
对于一般的问题，我们要求 $S(n)=\sum_{i=1}^nf(i)$。
构造 $g(n)$，$(f*g)(n)=\sum_{d|n}f(d)g\left(\frac nd\right)$，我们要求的是 $S(n)$，式子中可以从右边分离出 $f(n)$，再进一步推导，求出方程左边的前缀和，有
$$
\begin{aligned}
\sum_{i=1}^n(f*g)(i)&=\sum_{i=1}^n\sum_{d|i}f(d)g\left(\frac id\right)\\
&=\sum_{i=1}^n\sum_{x=1}^i\sum_{xy=i}f(x)g(y)\\
\end{aligned}
$$
因为每对数 $(x,y)$ 最多只会做一次贡献，所以我们可以直接枚举 $x\cdot y$ 的结果来统计贡献。
$$
\begin{aligned}
\sum_{i=1}^n(f*g)(i)&=\sum_{y=1}^n\sum_{xy\le n}f(x)g(y)\\
&=\sum_{y=1}^ng(y)\sum_{xy\le n}f(x)\\
&=\sum_{y=1}^ng(y)\sum_{i=1}^{\left\lfloor\frac ny\right\rfloor}f(i)
\end{aligned}
$$
我们在上面定义了 $S(n)$，所以式子继续化成了
$$
\sum_{i=1}^n(f*g)(i)=\sum_{y=1}^ng(y)S\left(\left\lfloor\frac ny \right\rfloor\right)
$$
这时我们可以整出 $S(n)$ 来了，当右侧式子中 $y=1$ 时，变为 $g(1)S(n)$。
移项得到
$$
g(1)S(n)=\sum_{i=1}^n(f*g)(i)-\sum_{y=2}^ng(y)S\left(\left\lfloor\frac ny\right\rfloor\right)
$$
这时我们看到了可以数论分块的 $\left\lfloor\frac ny\right\rfloor$，看上去像个递归过程。假定 $g(n)$ 的前缀和可以在 $O(1)$ 的复杂度内求出，且 $(f*g)(i)$ 的前缀和——即 $\sum_{i=1}^n(f*g)(i)$ 可以在 $O(1)\sim O(\sqrt n)$ 的时间内求出，则可以保证复杂度。
因此构造 $g$ 是一个重点。
3.复杂度
在构造之前，先考虑一下复杂度。
假定在理想状态下，求出 $g(n)$ 的前缀和和 $(f*g)(i)$ 的前缀和都是 $O(1)$ 的。
复杂度分析为
$$
T(n)=\sum_{i=1}^\sqrt n \sqrt{\left\lfloor\frac ni\right\rfloor}=\int_1^\sqrt n\sqrt{\frac nx}\mathrm dx
$$
式子可以化为 $\sqrt n\times x^{-\frac 12} \mathrm dx$，找到函数 $\sqrt n\times 2x^{\frac 12}=2\sqrt {nx}$ 的导数是它，然后带入
$$
\begin{aligned}
T(n)&=\int_1^\sqrt n\sqrt n\times x^{-\frac 12}\mathrm dx\\
&=2\sqrt{n\sqrt n}-2\sqrt n\\
&=O\left(n^{\frac 34}\right)
\end{aligned}
$$
同时，我们可以线筛预处理出前面一段的答案。此外，我们从 $S(n)$ 进入整个递归过程，这之后也只会调用自变量为 $\left\lfloor\frac ny\right\rfloor$ 的 $S$ 函数。接着，由于 $\left\lfloor\frac{\left\lfloor\frac ny\right\rfloor}{x}\right\rfloor=\left\lfloor\frac{n}{xy}\right\rfloor$，自变量只可能是 $\left\lfloor\frac ni\right\rfloor$ 的 $O(\sqrt n)$ 种取值中的一种，因此可以记忆化。
假设我们预处理出了前 $n^c$ 的答案，那么后面一部分积分的上界就是 $n^{1-c}$ 了，因为只有当 $k<c$ 时需要计算 $S\left(\frac{n}{n^{k}}\right)=S\left(n^{1-k}\right)$。
则
$$
\begin{aligned}
T(n)&=O(n^c)+\int_1^{n^{1-c}}\sqrt n\times x^{-\frac 12}\mathrm dx\\
&=O(n^c)+2\sqrt{n\cdot n^{1-c}}-2\sqrt n\\
&=O(n^c)+O\left(n^{1-\frac c2}\right)
\end{aligned}
$$
由基本不等式得二者取等时 $T(n)$ 有最小值，此时 $c=1-\frac c2,\ c=\frac 23$，$T(n)_\min=O\left(n^\frac 23\right)$。
4.构造
对于一些常见的函数如 $id(n),1(n),\epsilon(n)$，$g(n)$ 的前缀和比较好求，因此在杜教筛时优先考虑这些函数。
对于 $\varphi(n)$，由于 $(\varphi*1)(n)=n$，$(\varphi*1)(n)$ 和 $1(n)$ 的前缀和都比较好求，所以可以利用它。
把它带入移项后的式子：
$$
\begin{aligned}
g(1)S(n)&=\sum_{i=1}^n(f*g)(i)-\sum_{y=2}^ng(y)S\left(\left\lfloor\frac ny\right\rfloor\right)\\
1\times S(n)&=\sum_{i=1}^n(\varphi*1)(i)-\sum_{y=2}^n1\times S\left(\left\lfloor\frac ny\right\rfloor\right)\\
S(n)&=\frac{n(n+1)}2-\sum_{y=2}^nS\left(\left\lfloor\frac ny\right\rfloor\right)
\end{aligned}
$$
然后递归就可以了。
对于 $S(x)$ 的计算（$n$ 是常数）
当 $x\le \left(2^{31}-1\right)^{\frac 23}$ 时，可以直接返回已经存下的值，其余情况根据关于 $n$ 的分母，即 $\left\lfloor\frac nx\right\rfloor$ 的值查询是否记忆过这个答案，此时只需要存 $O\left(n^\frac 13\right)$ 种状态。
5.代码
#include<cstdio>
#include<cstring>
#define ll long long
const int mxm=2000000;
int pri[1<<20],cnt=0,n;
ll phi[2000100],mu[2000100];
bool is[2000100];
ll f[2000100],g[2000100];
bool used[2000];
ll calc(int x)
{
if(x<=mxm)
return phi[x];
int t=n/x;
if(used[t])
return g[t];
used[t]=1;
ll ans=(ll)(1ll+x)*x/2;
for(int i=2;i<=x;++i)
{
int nxt=x/(x/i);
ans-=calc(x/i)*(nxt-i+1);
i=nxt;
if(i==2147483647)
break;
}
return g[t]=ans;
}
ll calc1(int x)
{
if(x<=mxm)
return mu[x];
int t=n/x;
if(used[t])
return g[t];
used[t]=1;
ll ans=1;
for(int i=2;i<=x;++i)
{
int nxt=x/(x/i);
ans-=calc1(x/i)*(nxt-i+1);
i=nxt;
if(i==2147483647)
break;
}
return g[t]=ans;
}
int main()
{
is[0]=is[1]=1;
phi[1]=1;
mu[1]=1;
for(int i=2;i<=mxm;++i)
{
if(!is[i])
{
pri[++cnt]=i;
mu[i]=-1;
phi[i]=i-1;
}
for(int j=1;j<=cnt&&i*pri[j]<=mxm;++j)
{
is[i*pri[j]]=1;
if(i%pri[j]==0)
{
mu[i*pri[j]]=0;
phi[i*pri[j]]=pri[j]*phi[i];
break;
}
else
{
mu[i*pri[j]]=-mu[i];
phi[i*pri[j]]=(pri[j]-1)*phi[i];
}
}
mu[i]+=mu[i-1];
phi[i]+=phi[i-1];
}
int T;
scanf("%d",&T);
while(T--)
{
scanf("%d",&n);
memset(used,0,sizeof(used));
printf("%lld ",calc(n));
memset(used,0,sizeof(used));
printf("%lld\n",calc1(n));
}
return 0;
}
在数论题中最好通过带入最大的数（在本题中是 $2147483647$）来检验程序的运行时间和数据规模。在本题中，当循环到 $i=2147483647$ 时，++i 就会爆 int，在这里需要特判。
数组不要开小了！！
二、练习
洛谷 P3768 简单的数学题
求 $\sum_{i=1}^n\sum_{j=1}^nij\gcd(i,j)\bmod p$，$n\le 10^{10},p$ 是质数。
首先正常地化到反演的部分。
$$
\begin{aligned}
&\sum_{i=1}^n\sum_{j=1}^nij\gcd(i,j)\\
=&\sum_{i=1}^n\sum_{j=1}^nij\sum_{k=1}^nk\cdot[\gcd(i,j)=k]\\
=&\sum_{k=1}^nk\sum_{i=1}^n\sum_{j=1}^nij[\gcd(i,j)=k]
\end{aligned}
$$
令 $f(k)=\sum_{i=1}^n\sum_{j=1}^nij[\gcd(i,j)=k]$，构造
$$
\begin{aligned}
g(k)&=\sum_{k|d}f(d)\\
&=\sum_{k|d}\sum_{i=1}^n\sum_{j=1}^nij[\gcd(i,j)=k,2k,3k,\dots]\\
&=\sum_{i=1}^n\sum_{j=1}^nij[k|\gcd(i,j)]\\
&=\sum_{i=1}^\left\lfloor\frac nk\right\rfloor\sum_{j=1}^\left\lfloor\frac nk\right\rfloor ijk^2\\
&=k^2\frac{\left\lfloor\frac nk\right\rfloor^2\left(\left\lfloor\frac nk\right\rfloor+1\right)^2}4\\\
f(k)&=\sum_{k|d}\mu\left(\frac dk\right)g(d)\\
&=\sum_{i=1}^\left\lfloor\frac nk\right\rfloor\mu(i)g(ki)\\
&=\sum_{i=1}^{\left\lfloor\frac nk\right\rfloor}\mu(i)k^2i^2\frac{\left\lfloor\frac n{ki}\right\rfloor^2\left(\left\lfloor\frac n{ki}\right\rfloor+1\right)^2}4
\end{aligned}
$$
回到我们要求的答案
$$
\begin{aligned}
&\sum_{k=1}^nkf(k)\\
=&\sum_{k=1}^nk\sum_{i=1}^{\left\lfloor\frac nk\right\rfloor}\mu(i)k^2i^2\frac{\left\lfloor\frac n{ki}\right\rfloor^2\left(\left\lfloor\frac n{ki}\right\rfloor+1\right)^2}4\\
=&\sum_{T=1}^n\sum_{d|T}\mu(d)\frac{T^3}d\frac{\left\lfloor\frac nT\right\rfloor^2\left(\left\lfloor\frac nT\right\rfloor+1\right)^2}4\\
=&\sum_{T=1}^n\left(\sum_{d|T}\mu(d)\frac Td\right)T^2\frac{\left\lfloor\frac nT\right\rfloor^2\left(\left\lfloor\frac nT\right\rfloor+1\right)^2}4
\end{aligned}
$$
可以发现，被括号框起来的部分，是一个卷积形式，它是 $(\mu*id)(T)$，可以通过推导发现，$\mu*id=\varphi$。
推导：
$(\mu*id)(n)=\sum_{d|n}\mu(d)\frac nd$
当且仅当一个数 $x$ 是无平方因子数时，$\mu(x)\ne 0$。此时 $x$ 对卷积做出贡献。
假设 $n=\prod_{i=1}^kp_i^{c_i}$，那么会有 $2^k$ 个因数做出贡献。根据莫比乌斯函数的计算式，当枚举出来的数 $x$ 有奇数个质因数时，贡献为 $-\frac nx$；有偶数个质因数时，贡献为 $\frac nx$。可以发现，$\frac nx$ 是 $[x,n]$ 中能被 $x$ 整除的数的个数。
然后根据容斥(不会的话意会一下)就可以看出是 $n-$ 因数个数的形式。
当然严格证明还需要进一步推导。
对上述式子进一步化简得
$$
\sum_{T=1}^n\varphi(T)T^2\frac{\left\lfloor\frac nT\right\rfloor^2\left(\left\lfloor\frac nT\right\rfloor+1\right)^2}4
$$
此时面对 $n\le 10^{10}$ 的数据规模，还需要对 $\varphi(T)T^2$ 的前缀和进行处理，后面一段可以数论分块。经过记忆化后希望做到 $O\left(n^\frac 23\right)$。
考虑杜教筛求前缀和。$f(n)=\varphi(T)T^2$，我们需要构造出一个 $g(n)$ 使得 $(f*g)(n),\ g(n)$ 可以在较低的时间复杂度内算出前缀和。
当 $F(n)$ 是积性函数时，有 $(F\cdot G)*(F\cdot H)=(F\cdot (G*H))$，下面是推导
$$
\begin{aligned}
&((F\cdot G)*(F\cdot H))(n)\\
=&\sum_{d|n}(F\cdot G)(d)\times(F\cdot H)\left(\frac nd\right)\\
=&\sum_{d|n}F(d)\times G(d)\times F\left(\frac nd\right)\times H\left(\frac nd\right)\\
=&F(n)\times \sum_{d|n}G(d)\times H\left(\frac nd\right)\\
=&F(n)\times(G*H)(n)\\
=&(F\cdot(H*H))(n)
\end{aligned}
$$
那么把 $id^2(n)$ 看成是 $F(n)$，$\varphi(n)$ 看成是 $G(n)$，那么我们可以构造 $H(n)=1$，那么 $f(n)=\varphi(n)n^2$，$g(n)=n^2$，$(f*g)(n)=n^3$，然后套杜教筛计算就可以了。
[…] 【杜教筛学习笔记】 […]