关于一种积性函数前缀和的通用筛法的时间复杂度证明
注：本篇博客是从我知乎搬过来的，一方面是blog的排版不知道比知乎高到哪里去了，另外感觉知乎也不太适合发这种较为理论的内容，遂转战博客啦。
这是我的第一篇博客（顺便学习了各种格式和排版技巧），大家多多包涵！
最后，转载请注明出处，谢谢。
1.写在前面
自从学了杜教筛，突然对这一类问题产生了浓厚的兴趣。于是花了好几天时间，期间也学习了数论相关的许多理论，总算是努力没有白费，把这一算法的时间复杂度分析出来了。由于网上并没有这一筛法的相关介绍（仅有的一篇文献已在附录列出），在此写一篇详尽的分析，也便于以后忘记时查阅。
为了叙述方便，引入以下记号：
1. \(\left\lfloor x \right\rfloor \) 和 \(\left\lceil x \right\rceil \) 分别表示 \(x \) 的下取整和上取整，\([condition]\)表示条件 \(condition\)是否满足，若满足则值为1，不满足为0；
2.设需要考虑的任意积性函数为 \( f(x) \) ，前缀和记为 \( S(n)=\sum\limits_{i = 1}^n {f(i)} \) ，并令 \(m = \left\lceil {\sqrt n } \right\rceil \)；
3.以下用 \( p \) 表示质数，并设 \( maxfactor_i \) 表示 \( i \) 的最大质因子，特别地，\( maxfactor_p=p \)；
4.用 \(i = \prod\limits_{j = 1}^s {p_j^{{k_j}}} \)表示正整数 \(i\) 的质因数分解，其中 \(p_j\) 按升序排列， \(k_j\ge 1\) 。
2.通用筛法原理
2.1 总体思想
我们按如下的算法（以后简称算法1）分类进行统计。
在该算法中，我们从2开始从小到大遍历每一整数 \(i\) ，设 \(F=maxfactor_i\) ，我们进行如下2类操作：
操作一：对每一大于 \(F\) 的质数 \(p\) ，计算 \(f(ip)=f(i)f(p)\) ，并累加进和中，直到 \(ip>n\) 为止；
操作二：如果 \(i\) 含有至少2个质因子 \(F\) ，那么将 \(f(i)\) 累加进和中。
那么全部整数 \(i\) 遍历完后，所得的总和再加上 \(\sum\limits_{p \le n} {f(p)}\) ，以及 \(f(1)=1\) ，便是 \(S(n)\) 。
命题1：上述算法求得的答案正好是 \(S(n)\) 。
证明：只需证明对于每个整数 \(i\) ，都不重不漏的正好统计了一次即可。
（1）对于 \(f(1)\) ，显然只被统计了一次；
（2）所有 \(f(p)\) 均被统计了恰好一次，因为操作一以及操作二均不会统计 \(f(p)\) ；
（3）其它情况，设 \(i = \prod\limits_{j = 1}^s {p_j^{{k_j}}}\) ，其中 \(s \ge 2\) ，那么分两种情况：
情况1：若 \(k_s=1\) ，那么操作一中 \(\prod\limits_{j = 1}^{s-1} {p_j^{{k_j}}}\) 必然会转移到 \(i\) 而将 \(f(i)\) 统计，且其它数不会转移到 \(i\) ；显然操作二也不会统计 \(f(i)\) ；因而被恰好统计了一次；
情况2：若 \(k_s>1\) ，那么操作一不会统计 \(f(i)\) ，但操作二会统计，因而被恰好统计了一次。
证毕。
接来下来考虑一个关键问题：如何统计操作一中的 \(\sum\limits_{F < p,ip \le n} {f(p)}\) 。为此我们定义 \(S'(x) = \sum\limits_{p \le x} {f(p)}\) ，那么只要转化为求 \(S(x)\) 即可。
命题2：上述算法只需要用到 \(S_1=\{ S'(x):1 \le x \le m\} \) 以及 \(S_2=\{ S'\left( {\left\lfloor {\dfrac{n}{x}} \right\rfloor } \right):1 \le x \le m\} \) 。
证明： \(\sum\limits_{F < p,ip \le n} {f(p)} =S'\left( {\left\lfloor {\dfrac{n}{i}} \right\rfloor } \right)-S'(F) \)。
对于 \(S'\left( {\left\lfloor {\dfrac{n}{i}} \right\rfloor } \right)\) ，若 \(i\le m\) ，则属于 \(S_2\) ，否则属于 \(S_1\) ；
对于 \(S'(F)\) ，由 \(F^2\le ip\le n\) 得 \(F\le m\) ，属于 \(S_1\) 。
对于最后还需要加上的全体质数，显然 \(S'(n)\) 属于 \(S_2\) 。
证毕。
2.2 积性函数素数前缀和算法
现在问题已经归结为求1到\(t\)内所有素数\(p\)的前缀和了。我们采用如下算法（以下简称算法2）进行求解。
对于求 \(S'(i)\) ，采用欧拉筛法的思想，初始时将每个数都视为质数，将 \(f(j),2 \le j \le i\) 全部加到 \(S'(i)\) 中。然后从小到大遍历每一质数 \(p\) ，令 \(i\) 从 \(n\) 开始，从大到小执行以下操作： \[S'(i) - = {\rm{ }}\left( {S'\left( {\left\lfloor {\dfrac{i}{p}} \right\rfloor } \right){\rm{ }} - {\rm{ }}S'(p - 1)} \right)f\left( p \right)\] 直到 \(i < p^2\) 为止。
定理3：当 \(f(i)\) 是完全积性函数（即 \(f(ij) = f(i)f(j),\forall i,j\) ）时，上述算法求出的每一个 \(S'(i)\) 是正确的。
证明：我们证明一个循环不变式，即算法执行完当前质数 \(p\) 后，对于每一个 \(i\) ， \(S'(i)\) 均扣除且仅扣除了所有的 \(f(x)\) ，其中 \(x\) 是合数且 \(maxfactor_x \le p\) 。
初始时，显然是正确的。考虑当前的质数 \(p\) ： \[{\rm{ }}\left( {S'\left( {\left\lfloor {\dfrac{i}{p}} \right\rfloor } \right){\rm{ }} - {\rm{ }}S'(p - 1)} \right)f\left( p \right)\\=f(p)\sum\limits_{j = p}^{\left\lfloor {i/p} \right\rfloor } {[maxfacto{r_j} \ge p]f(j)} \\=\sum\limits_{j = p}^{\left\lfloor {i/p} \right\rfloor } {[maxfacto{r_j} = p]f(jp)} \] 因而质数 \(p\) 对应的一轮循环执行完后，所有 \(x\) 是合数且\(maxfactor_x = p\) 的数会被筛去，且仅筛去了这些数。证毕。
然而，上述证明具有极为苛刻的条件， \(f(i)\) 是完全积性函数。如果这一条件不满足，结论是不成立的。例如， \(f(p)=-1\) （即莫比乌斯函数 \(\mu\) ）就不满足，可以令 \(n=4\) 执行上述流程进行验证。
然而好消息是，函数 \(f(p)=p^a\) 满足这一要求，于是我们可以采取将一般积性函数表示为 \(p^a\) 的线性组合的形式，对绝大多数积性函数均能这么做。例如： \[\mu(p)=-1\] \[\phi(p)=p-1\] \[\sigma_0(p)=2\] \[\sigma_1(p)=p+1\] 等等。我们对每一项 \(p^a\) 执行一次算法，并将结果保存起来，最后将所有项结果线性叠加即可。注意： \(S'(i)\) 仅统计质数的 \(f(i)\) 的和，与 \(f(p^k),k>1\) 完全无关。在这一步中完全不必关心它， \(f(p^k)\) 在算法1中才有用。
3.通用筛法实现
根据上述原理，这里仍分两部分考虑实现。
3.1 算法2的实现
首先考虑计算 \(S'(i)\) 的部分。根据结论2，我们只需要两个大小均为 \(m\) 的数组，其中 \(a[i]=S'(i)\) ， \(b[i]=S'\left( {\left\lfloor {\dfrac{n}{i}} \right\rfloor } \right)\) ，这只需要 \(\Theta (m)\) 的空间。
对于筛法的公式，不难发现以下三种情况：
（1）若 \(\dfrac{i}{p} > \left\lfloor {\sqrt n } \right\rfloor\) ，那么 \(S'(i)\) 便对应了 \(b[n/i]\)，而\(S'\left( {\left\lfloor {\dfrac{i}{p}} \right\rfloor } \right)\) 对应了 \(b[n/i*p]\) ；
（2）若 \(i > \left\lfloor {\sqrt n } \right\rfloor \) 且 \(\dfrac{i}{p} \le \left\lfloor {\sqrt n } \right\rfloor \) ，那么 \(S'(i)\) 便对应了 \(b[n/i]\)，而 \(S'\left( {\left\lfloor {\dfrac{i}{p}} \right\rfloor } \right)\) 便对应了 \(a[i/p]\) ；
（3）若 \(i \le \left\lfloor {\sqrt n } \right\rfloor \) ，则 \(S'(i)\) 便对应了 \(a[i]\)，而 \(S'\left( {\left\lfloor {\dfrac{i}{p}} \right\rfloor } \right)\) 便对应了 \(a[i/p]\) 。
由此可见，更新 \(a,b\) 数组不会用到 \(a,b\) 数组之外的元素 \(S'(i)\) 。
3.2 算法1的实现
计算完 \(S'(i)\) 后，根据之前的算法，我们要枚举每一整数 \(i\) ，计算 \(f(i)\left(S'\left( {\left\lfloor {\dfrac{n}{i}} \right\rfloor } \right)-S'(maxfactor_i)\right)\) 。为了快速枚全部 \(i\) ，可以考虑递归实现。
设函数形式为 \(fun(x,p)\) ，表示 \(x\) 的最大质因子是 \(p\) ，那么我们先计算上面的式子，然后转移到更大的 \(x\) ，即递归的调用 \(fun(y,q)\) ，其中 \(y=xq^k\) ， \(q>p\) 是更大质数。若 \(k\ge 2\) ，那么我们顺便计算 \(f(y)\) 。这样我们就能在 \(O(1)\) 的时间复杂度处理一个 \(i\) 了。
最后，采用一个重要的剪枝：若 \(i \times maxfactor_i>n\) ，显然这样的 \(i\) 是无需考虑的，因为不会对答案有贡献。于是在递归调用时加上这一条件即可。显然不加这一优化，我们需要遍历1到 \(n\) 中的所有 \(i\) ，时间复杂度将退化为线性。
4.通用筛法的时间复杂度证明
终于到了最具有理论价值（研究了最久）的部分了。同样地，我们将分别考虑两个算法各自的时间复杂度。为此，需要引入许多数论中的定理，尤其是素数计数相关的定理。
4.1 素数幂和相关引理
引理4：1到 \(n\) 中质数的个数具有估计式 \(\pi(n)=\Theta \left( {\dfrac{n}{{\ln n}}} \right)\) ，同时隐含常数为1。
引理的证明从略。
引理5：\(\sum\limits_{p \le n} {{p^2}} = \Theta \left( {\dfrac{{{n^3}}}{{\ln n}}} \right)\) ，同时隐含常数为\(\dfrac 1 3\)。
引理的证明从略，以上两引理可在最后参考文献中查阅。
引理6：对任意整数 \(k\) ，当 \(n\rightarrow\infty\)时，\(\displaystyle \sum\limits_{p \le n} {\dfrac{1}{{{p^k}}}} = \int_1^n {\dfrac{{\mathrm{d}\pi (x)}}{{{x^k}}}} \) 。
该公式称为Stieltjes积分公式，相关证明可查阅文献。
推论7： \(\displaystyle \sum\limits_{ p > n} {\dfrac{1}{{{p^2}}}} =\Theta \left( {\dfrac{{1 }}{{n\ln n}}} \right)\)，同时隐含常数为1 。
由于并没有查到相关的公式，因此以下对该结论证明。
证明：利用引理6得 \[\sum\limits_{p > n} {\dfrac{1}{{{p^2}}}} = \int_n^\infty {\dfrac{{\mathrm{d}\pi (x)}}{{{x^2}}}} = \dfrac{{\pi (x)}}{{{x^2}}}\bigg|_n^\infty + 2\int_n^\infty {\dfrac{{\mathrm{d}x}}{{{x^2}\ln x}}}\] 这一步是在 \(n\rightarrow\infty\)的前提下成立的，此时分部积分中 \(\pi(x)\) 可替换为 \(\dfrac{x}{\ln x}\) 。继续计算可以得到 \[\sum\limits_{p > n} {\dfrac{1}{{{p^2}}}} =- \dfrac{1}{{n\ln n}} - 2\int_0^{\dfrac{1}{n}} {\dfrac{{\mathrm{d}t}}{{\ln t}}}=-\dfrac{1}{{n\ln n}} - 2Li\left(\dfrac 1 n\right)\] 其中右边积分用 \(t=\dfrac 1 x\) 替换得到， \(\displaystyle Li(x)=\int_0^{x} {\dfrac{{\mathrm{d}x}}{{\ln x}}}\) 称为对数积分。根据 \(Li(x)\) 的泰勒展开式， \[Li(x) = \dfrac{x}{{\ln x}}\sum\limits_{k = 0}^\infty {\dfrac{{k!}}{{{{(\ln x)}^k}}}}\] 因而 \[Li\left(\dfrac 1 n\right)\sim-\dfrac 1 {n\ln n}\] 故 \[\sum\limits_{ p > n} {\dfrac{1}{{{p^2}}}} =\Theta \left( {\dfrac{{1 }}{{n\ln n}}} \right)\] 证毕。
4.2 算法2的时间复杂度证明
前面列举了一堆引理，现在终于到核心部分了。下面定理给出了算法2中的时间复杂度的精确估计。
定理8：计算 \(S'(i)\) 的时间复杂度为 \(\Theta \left( {\dfrac{{{n^{3/4}}}}{{\ln n}}} \right)\) ，精确到常数，其基本语句执行次数为 \( \dfrac{{{32n^{3/4}}}}{{3\ln n}}\) 。
证明：我们将1到 \(\sqrt n\) 中的质数分为两部分，对应算法2的流程，分别考察其对时间复杂度的贡献。
第一部分： \(2\le p \le \sqrt m\)
此时对于每一个 \(1\le i\le \sqrt m\) ， \(b[i]\) 均需要花 \(\Theta (1)\) 的时间更新。
而对于 \(p^2\le i\le \sqrt m\) ， \(a[i]\) 均需要花 \(\Theta (1)\) 的时间更新。
因而对于这个 \(p\) ，总时间复杂度为 \(\Theta (2m-p^2)\) 。
第二部分： \(\sqrt m < p \le m\)
此时 \(a[i]\) 不需要继续更新；
而对于每一个 \(\dfrac n i\ge p^2\) ，即\(1\le i\le \dfrac n {p^2}\) ， \(b[i]\) 均需要花 \(\Theta (1)\) 的时间更新。
因而对于这个 \(p\) ，总时间复杂度为 \(\Theta \left( {\dfrac{n}{{{p^2}}}} \right)\) 。
将两部分对 \(p\) 求和，得总时间复杂度为 \[T(n)=\sum\limits_{p \le \sqrt m } {(2m - {p^2})} + \sum\limits_{\sqrt m < p \le m} {\left( {\dfrac{n}{{{p^2}}}} \right)} \\= \left(2m\sum\limits_{p \le \sqrt m } 1\right) - \left(\sum\limits_{p \le \sqrt m } {{p^2}}\right) + \left(n\sum\limits_{\sqrt m < p \le m} {\dfrac{1}{{{p^2}}}} \right)\] 根据引理4、引理5和推论7，得 \[T(n)=\Theta\left(2m\dfrac{{\sqrt m }}{{\ln \sqrt m }} - \dfrac{{{{\left( {\sqrt m } \right)}^3}}}{{3\ln \sqrt m }} + n\dfrac{1}{{\sqrt m \ln \sqrt m }}\right) \\=\Theta\left(\dfrac {8m\sqrt m} {3\ln {\sqrt m}} \right)=\Theta\left(\dfrac {32n^{3/4}} {3\ln n} \right)\] 由于上述证明使用的均为等价形式，因而时间复杂度估计使用的是记号 \(\Theta\) 。证毕。
以下列出了一些典型 \(n\) 时，理论分析和实际测试的计算次数。
当 \(n=10^{11}\) 时，计算结果7489万，测试结果8029万；
当 \(n=10^{12}\) 时，计算结果3.86亿，测试结果4.13亿。
可以看出，两者是十分接近的。差距主要在于上述引理是在 \(n\) 趋于无穷大时才成立。同时，由于已经考虑常数，故 \(n=10^{11}\) 时算法执行时间仅为300到500毫秒左右。
4.3 Dickman函数
接下来考虑算法1。为此，仍需要引入一些引理。
引理9：对于给定的实常数 \(0 < a < 1\) ，令 \(Q(n) = \{ i \le n:maxfactor_i \le {i^a}\}\) ，那么 \(\left| Q(n) \right| \sim n\rho \left(\dfrac 1 a\right)\) ，其中 \(\rho (x)\) 称为迪克曼函数(Dickman function)。
性质10：对于任意 \(x>1\) ， \(\rho (x)>0\) 且随 \(x\) 增大迅速衰减，具有近似式 \(\rho (x) \approx {x^{ - x}}\) 。
以上两条引理证明从略，可在最后参考文献中查阅。
上述引理表明，当 \(a\) 较小时， \(Q\) 集合中的元素个数急剧衰减；即对于一个大整数 \(i\)，它的最大质因子通常都比较大。
4.4 关于算法1的时间复杂度
现在回到算法1，之前已说明算法需要枚举每一个 \(i\) ，而通过一个剪枝 \(i\times maxfactor_i\le n\) 可以大大减少枚举的 \(i\) 的数量，这一原理便可由上面指出。
然而通过一番分析，仍然不能给出较好的时间复杂度证明。最终不幸的事情还是发生了，我们有如下定理：
定理11：令 \(M(n) = \{ i:i \times maxfacto{r_i} \le n\}\) ，对于任意的 \(\alpha<1\) ，有 \(\left| M(n) \right| = \Omega ({n^\alpha })\) 。
证明：根据引理9， 令 \(Q(n) = \{ i \le n:maxfactor_i \le {i^{1-\alpha}}\}\) ，那么 \(\left| Q(n^\alpha) \right| \sim n^\alpha\rho \left(\dfrac 1 {1-\alpha}\right)\) 。对于每一个数 \(x\in Q(n^\alpha)\) ，显然 \(x \times maxfacto{r_x} \le n\) ，这就意味着1到 \(n\) 中至少有 \(n^\alpha\rho \left(\dfrac 1 {1-\alpha}\right)\) 个数在 \(M(n)\) 中，即 \(\left| M(n) \right| = \Omega ({n^\alpha })\) 。证毕。
尽管理论十分悲观，但是不妨代入一个 \(\alpha=3/4\) 算一下。可以查阅到 \(\rho (4) \approx 5 \times 10^{-3}\) ，当 \(n=10^{12}\) 时， \(n^\alpha\rho \left(\dfrac 1 {1-\alpha}\right) \approx 5 \times 10^6\) ，仍是十分小的，这正是得益于\(\rho (x)\) 衰减速度非常快，使得在平常使用时该算法仍具有价值。
注意到定理11给出的明显是一个下界，实际的值要大不少。下面列出了一些典型 \(n\) 时，实际测试结果。
当 \(n=10^{11}\) 时，测试结果3400万；
当 \(n=10^{12}\) 时，测试结果1.88亿。
可以看出，这一值甚至小于算法2的计算量。但是由于递归等原因，算法2往往跑的更快。
命题12：如果对每一质数 \(p\leq \sqrt m\) ， \(i \le n\) 中满足 \(maxfactor_i=p\) 且 \(ip \le n\) 的 \(i\) 平均小于 \(\sqrt n\) 个，那么算法1就具有 \(O\left(\dfrac {n^{3/4}} {\log n} \right)\) 的时间复杂度。
证明：设 \(maxfactor_i=p\) ，将\(p\) 分为两段分别考虑。
当 \(\sqrt m < p\leq m\) 时，算法1的时间复杂度必然小于 \(\displaystyle \sum\limits_{\sqrt m < p \le m} {\dfrac{n}{{{p^2}}}}\) ，因为对于每个 \(p\) 只有不超过 \(\left\lfloor {\dfrac{n}{{{p^2}}}} \right\rfloor\) 个 \(i\) 满足 \(maxfactor_i=p\) 。根据推论7，这一求和式不超过 \(O\left(\dfrac {n^{3/4}} {\log n} \right)\) 。
当 \(p\le \sqrt m\) 时，根据引理5， \(p\) 的个数有 \(\Theta \left( {\dfrac{{{n^{1/4}}}}{{\ln n}}} \right)\) ，再根据命题假设，计算量不超过 \(O\left(\dfrac {n^{3/4}} {\log n} \right)\) 。
综上可知，整个算法1的时间复杂度不超过 \(O\left(\dfrac {n^{3/4}} {\log n} \right)\) 。证毕。
通过打表知，在\(n \le 10^{12}\) 下命题2的假设满足，因而可以“认为”整个求积性函数前缀和算法的时间复杂度为 \(\Theta \left( {\dfrac{{{n^{3/4}}}}{{\ln n}}} \right)\) 。
5.一些模板题目
HDU5091 求1到\(n\)中素数个数
LOJ6202 求素数的和
51nod1239 求欧拉函数前缀和
51nod1244 求莫比乌斯函数前缀和
以上是4个模板题。其中求欧拉函数的题使用本文的方法时间上完全可以和杜教筛持平。
SPOJ 求\(\sigma_0(i^3)\)之和
这道题采用本文方法便瞬间退化为一道毫无思维难度的题了。
ZOJ5340 求一般意义下的积性函数前缀和
值得注意的是，此题使用本文方法虽然直接退化为了一道模板题，但却并不能过（若有人使用此方法能过麻烦联系我~），其原因在于多组数据。而杜教筛在面对多组数据时是很灵活的，只需改变分块大小（预处理更多的数）即可。因此面对实际问题要选用合适的方法。
6.后记
尽管最后的结果表明算法时间复杂度似乎是线性的（即对于任意\(\alpha < 1\)，该算法时间复杂度均高于\(n^{\alpha}\)），这已经将算法的理论价值大打折扣，但是：一方面，该算法仍具有巨大的应用价值。例如，用该算法求莫比乌斯函数前缀和，在oj上的提交结果为 \(n=10^{10}\) 下运行350ms，已接近杜教筛的运行时间；另一方面，算法2具有严格的理论时间复杂度证明（精确到常数），从而具有很强的理论价值。例如，可用来求1到 \(n\) 中的质数个数，或1到 \(n\) 中的全部质数之和等。
个人感觉，算法2是一个非常优美的算法，因为其时间复杂度的证明过程中，恰好最后的三项均为 \(\Theta \left( {\dfrac{{{n^{3/4}}}}{{\ln n}}} \right)\) ，即各个组成部分时间复杂度平衡。这就类似于杜教筛的 \(n^{2/3}\) 预处理来配平时间复杂度一样，不同之处在于，这里没有任何人为因素，完全是天然平衡的，而且是三项求和平衡。
总之，希望该算法的价值能够在日后不断体现出来，也更希望能有巨佬将算法1改进，使得理论时间复杂度真正达到 \(\Theta \left( {\dfrac{{{n^{3/4}}}}{{\ln n}}} \right)\) !
最后，由于本人才疏学浅，文章中如有任何不当之处，欢迎批评指正，谢谢！
