杜教筛原理，实现与时间复杂度分析

引例

洛谷 P4213 【模板】杜教筛

题目描述

给定一个正整数，求

$$an{s}_1=\sum _{i=1}^n\phi (i)$$

$$an{s}_2=\sum _{i=1}^n\mu (i)$$

输入格式

本题单测试点内有多组数据。

输入的第一行为一个整数，表示数据组数 $T$。

接下来 $T$ 行，每行一个整数 $n$，表示一组询问。

输出格式

对于每组询问，输出一行两个整数，分别代表 $an{s}_1$ 和 $an{s}_2$。

输入输出样例 #1

输入 #1

6
1
2
8
13
30
2333

输出 #1

1 1
2 0
22 -2
58 -3
278 -3
1655470 2

说明/提示

数据规模与约定

对于全部的测试点，保证 $1\le T\le 10$，$1\le n<2^{31}$。对于全部的测试点，保证 $1\le T\le 10$，$1\le n<2^{31}$。

题中要求以低于线性的时间复杂度对数论函数求和。

原理

将例写为更一般的形式：

求 $$S(n)=\sum _{i=1}^{n}f(i)$$ ，其中，$f$ 是数论函数。

构造函数 $g$ 并令 $h=f\ast g=g\ast f$ ，
故，
$$\begin{array}{rl}\sum _{i=1}^{n}h(i)& =\sum _{i=1}^n\sum _{d|i}g(d)f(\frac{i}{d})\\ & =\sum _{d=1}^n\sum _{d|i}^{n}g(d)f(\frac{i}{d})\\ & =\sum _{d=1}^{n}g(d)\sum _{d|i}^{n}f(\frac{i}{d})\\ & =\sum _{d=1}^{n}g(d)\sum _{i=1}^{\lfloor n/d\rfloor }f(i)\\ & =\sum _{d=1}^{n}g(d)S(\lfloor n/d\rfloor )\\ & =g(1)S(n)+\sum _{i=2}^{n}g(i)S(\lfloor n/i\rfloor )\end{array}$$
因此，
$$g(1)S(n)=\sum _{i=1}^{n}h(i)-\sum _{i=2}^{n}g(i)S(\lfloor n/i\rfloor )$$

采用数论分块加速第二个求和运算，这要求：

• 可以快速计算 ${\sum }_{i=1}^{n}h(i)$ ；
• 可以快速计算 ${\sum }_{i=l}^{r}g(i)$ ；

实现

实现时构造出来的 $h,g$ 已经满足了上述要求，采用记忆化搜索计算 $S(n)$ ，数论分块加速求和。
众所周知，$\mu \ast I=ϵ$，$\phi \ast I=id$ 。
–>To learn more about it
于是，带入到杜教筛公式中，
$$I(1)S(n)=\sum _{i=1}^nϵ(i)-\sum _{i=2}^{n}I(i)S(\lfloor n/i\rfloor )$$
$$I(1)S(n)=\sum _{i=1}^{n}id(i)-\sum _{i=2}^{n}I(i)S(\lfloor n/i\rfloor )$$
即，
$$S(n)=1-\sum _{i=2}^{n}S(\lfloor n/i\rfloor )$$
$$S(n)=\frac{n(n+1)}{2}-\sum _{i=2}^{n}S(\lfloor n/i\rfloor )$$
考虑进一步优化，对于前 $k$ 项使用线性筛直接计算，可以证明 $k=O(n^{2/3})$ 最优（证明参见最后一部分）。
在实现中，一般使用 unordered_map/gp_hash_table 进行记忆化。由于常数较大，一般将预处理部分的 $k$ 在不炸空间的情况下开得尽量大一些。
以下为C++参考实现，尽管全开了long long并且未经过卡常，常数依然中规中矩，完全可以接受。record link 。

#include <bits/stdc++.h>
#define int long long
using namespace std;const int K=5e6; 
int t,n,vh1[K],vh2[K];
bool notprime[K];
vector<int>prime;
unordered_map<int,int>s1,s2;inline int read(){int s=0;bool f=0;char c=getchar();while(!isdigit(c)){if(c=='-')f^=1;c=getchar();}while(isdigit(c)){s=(s<<1)+(s<<3)+(c^48);c=getchar();}return f==0?s:-s;
}void pre(){vh1[1]=1;vh2[1]=1;for(int i=2;i<K;++i){if(!notprime[i]){prime.push_back(i);vh1[i]=i-1;vh2[i]=-1;}for(int p:prime){if(i*p>=K)break;notprime[i*p]=1;vh1[i*p]=vh1[i]*vh1[p];vh2[i*p]=-vh2[i];if(i%p==0){vh1[i*p]=vh1[i]*p;vh2[i*p]=0; break;}}}for(int i=1;i<K;++i)vh1[i]+=vh1[i-1],vh2[i]+=vh2[i-1];
}int S1(int n){if(n<K)return vh1[n];if(s1[n]!=0)return s1[n];int res=(n+1)*n/2,l=2,r;while(l<=n){r=n/(n/l);res-=S1(n/l)*(r-l+1);l=r+1;}return s1[n]=res;
}int S2(int n){if(n<K)return vh2[n];if(s2[n]!=0)return s2[n];int res=1,l=2,r;while(l<=n){r=n/(n/l);res-=S2(n/l)*(r-l+1);l=r+1;}return s2[n]=res;
}signed main(){pre();t=read();while(t--){n=read();cout<<S1(n)<<" "<<S2(n)<<endl;}return 0;
}

时间复杂度分析

此处伪证极多，点名《算法竞赛》中的“高阶小量”伪证。

先考虑有哪些 $S(i)$ 会被计算。
设记搜时 $S(i)$ 会使用 R ( i ) = { j ∣ j = ⌊ i k ⌋ , k ∈ N + , k > 1 } R(i)=\{j|j=\lfloor \frac ik \rfloor,k∈\mathbb{N^+},k>1\} R(i)={j∣j=⌊ki​⌋,k∈N+,k>1} 中的元素的 $S$ 。

下证若 $m\in R(n)$ ，则 $R(m)\subseteq R(n)$ 。
对于任意 $i\in R(m)$ , $i=\lfloor \frac{m}{x}\rfloor$ ，
又 $m=\lfloor \frac{n}{y}\rfloor$ ，
故 $i=\lfloor \frac{\lfloor \frac{n}{y}\rfloor }{x}\rfloor =\lfloor \frac{n}{xy}\rfloor \in R(n)$ 。

由上知只有 $i\in R(n)$ 才会被计算。

下令 $T(n)$ 为计算 $S(n)$ 的时间复杂度，并忽略计算 $h,g$ 相关的消耗。则由数论分块相关可知 $T(n)=\sqrt{n}$ 。

类似数论分块，若 $i=\lfloor \frac{n}{x}\rfloor \in R(n)$，
则 $i\le \sqrt{n}$ 或 $ i=\lfloor \frac nj \rfloor,j<\sqrt n$ 。

对于朴素的杜教筛（不使用线性筛优化），其时间复杂度为，
$$\begin{array}{rl}\sum _{i\in R(n)}T(i)& =\sum _{i\in R(n)}\sqrt{i}\\ & =\sum _{i=1}^{\lfloor \sqrt{n}\rfloor }\sqrt{i}+\sum _{j=2}^{\lfloor \sqrt{n}\rfloor }\sqrt{\frac{n}{j}}\\ & ＜\sum _{i=1}^{\lfloor \sqrt{n}\rfloor }\sqrt{i}+\sqrt{\frac{n}{i}}\\ & \approx {\int }_0^{\sqrt{n}}(\sqrt{x}+\sqrt{\frac{n}{x}})dx\\ & =n^{3/4}\end{array}$$
即为 $O(n^{3/4})$ 。

若使用线性筛预处理前 $k$ 项，预处理部分时间复杂度为 $O(k)$ 。当 $k>\sqrt{n}$ 时时间复杂度变为，
$$\begin{array}{rl}k+\sum _{i\in R(n),i>k}T(i)& =k+\sum _{i\in R(n),i>k}\sqrt{i}\\ & =k+\sum _{j=2}^{\lfloor n/k\rfloor }\sqrt{\frac{n}{j}}\\ & ＜k+\sum _{i=1}^{\lfloor n/k\rfloor }\sqrt{\frac{n}{i}}\\ & \approx k+{\int }_0^{n/k}\sqrt{\frac{n}{x}}dx\\ & =k+\frac{n}{\sqrt{k}}\end{array}$$
为求上式极值，以 $k$ 为自变量求导，令 $\frac{dO}{dk}=1-\frac{1}{2}n{k}^{-\frac{3}{2}}=0$ ,
得到 $k=O(n^{2/3})$ ，此时总时间复杂度为 $O(n^{2/3})$ 。