注：本篇博客是从我知乎搬过来的，一方面是blog的排版不知道比知乎高到哪里去了，另外感觉知乎也不太适合发这种较为理论的内容，遂转战博客啦。
这是我的第一篇博客（顺便学习了各种格式和排版技巧），大家多多包涵！
最后，转载请注明出处，谢谢。
1.写在前面
自从学了杜教筛，突然对这一类问题产生了浓厚的兴趣。于是花了好几天时间，期间也学习了数论相关的许多理论，总算是努力没有白费，把这一算法的时间复杂度分析出来了。由于网上并没有这一筛法的相关介绍（仅有的一篇文献已在附录列出），在此写一篇详尽的分析，也便于以后忘记时查阅。
为了叙述方便，引入以下记号：
1. (leftlfloor x ight floor ) 和 (leftlceil x ight ceil ) 分别表示 (x ) 的下取整和上取整，([condition])表示条件 (condition)是否满足，若满足则值为1，不满足为0；
2.设需要考虑的任意积性函数为 ( f(x) ) ，前缀和记为 ( S(n)=sumlimits_{i = 1}^n {f(i)} ) ，并令 (m = leftlceil {sqrt n } ight ceil )；
3.以下用 ( p ) 表示质数，并设 ( maxfactor_i ) 表示 ( i ) 的最大质因子，特别地，( maxfactor_p=p )；
4.用 (i = prodlimits_{j = 1}^s {p_j^{{k_j}}} )表示正整数 (i) 的质因数分解，其中 (p_j) 按升序排列， (k_jge 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) 遍历完后，所得的总和再加上 (sumlimits_{p le n} {f(p)}) ，以及 (f(1)=1) ，便是 (S(n)) 。
命题1：上述算法求得的答案正好是 (S(n)) 。
证明：只需证明对于每个整数 (i) ，都不重不漏的正好统计了一次即可。
（1）对于 (f(1)) ，显然只被统计了一次；
（2）所有 (f(p)) 均被统计了恰好一次，因为操作一以及操作二均不会统计 (f(p)) ；
（3）其它情况，设 (i = prodlimits_{j = 1}^s {p_j^{{k_j}}}) ，其中 (s ge 2) ，那么分两种情况：
情况1：若 (k_s=1) ，那么操作一中 (prodlimits_{j = 1}^{s-1} {p_j^{{k_j}}}) 必然会转移到 (i) 而将 (f(i)) 统计，且其它数不会转移到 (i) ；显然操作二也不会统计 (f(i)) ；因而被恰好统计了一次；
情况2：若 (k_s>1) ，那么操作一不会统计 (f(i)) ，但操作二会统计，因而被恰好统计了一次。
证毕。
接来下来考虑一个关键问题：如何统计操作一中的 (sumlimits_{F < p,ip le n} {f(p)}) 。为此我们定义 (S'(x) = sumlimits_{p le x} {f(p)}) ，那么只要转化为求 (S(x)) 即可。
命题2：上述算法只需要用到 (S_1={ S'(x):1 le x le m} ) 以及 (S_2={ S'left( {leftlfloor {dfrac{n}{x}} ight floor } ight):1 le x le m} ) 。
证明： (sumlimits_{F < p,ip le n} {f(p)} =S'left( {leftlfloor {dfrac{n}{i}} ight floor } ight)-S'(F) )。
对于 (S'left( {leftlfloor {dfrac{n}{i}} ight floor } ight)) ，若 (ile m) ，则属于 (S_2) ，否则属于 (S_1) ；
对于 (S'(F)) ，由 (F^2le iple n) 得 (Fle 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) - = { m{ }}left( {S'left( {leftlfloor {dfrac{i}{p}} ight floor } ight){ m{ }} - { m{ }}S'(p - 1)} ight)fleft( p ight)] 直到 (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) ： [{ m{ }}left( {S'left( {leftlfloor {dfrac{i}{p}} ight floor } ight){ m{ }} - { m{ }}S'(p - 1)} ight)fleft( p ight)\=f(p)sumlimits_{j = p}^{leftlfloor {i/p} ight floor } {[maxfacto{r_j} ge p]f(j)} \=sumlimits_{j = p}^{leftlfloor {i/p} ight floor } {[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( {leftlfloor {dfrac{n}{i}} ight floor } ight)) ，这只需要 (Theta (m)) 的空间。
对于筛法的公式，不难发现以下三种情况：
（1）若 (dfrac{i}{p} > leftlfloor {sqrt n } ight floor) ，那么 (S'(i)) 便对应了 (b[n/i])，而(S'left( {leftlfloor {dfrac{i}{p}} ight floor } ight)) 对应了 (b[n/i*p]) ；
（2）若 (i > leftlfloor {sqrt n } ight floor ) 且 (dfrac{i}{p} le leftlfloor {sqrt n } ight floor ) ，那么 (S'(i)) 便对应了 (b[n/i])，而 (S'left( {leftlfloor {dfrac{i}{p}} ight floor } ight)) 便对应了 (a[i/p]) ；
（3）若 (i le leftlfloor {sqrt n } ight floor ) ，则 (S'(i)) 便对应了 (a[i])，而 (S'left( {leftlfloor {dfrac{i}{p}} ight floor } ight)) 便对应了 (a[i/p]) 。
由此可见，更新 (a,b) 数组不会用到 (a,b) 数组之外的元素 (S'(i)) 。
3.2 算法1的实现
计算完 (S'(i)) 后，根据之前的算法，我们要枚举每一整数 (i) ，计算 (f(i)left(S'left( {leftlfloor {dfrac{n}{i}} ight floor } ight)-S'(maxfactor_i) ight)) 。为了快速枚全部 (i) ，可以考虑递归实现。
设函数形式为 (fun(x,p)) ，表示 (x) 的最大质因子是 (p) ，那么我们先计算上面的式子，然后转移到更大的 (x) ，即递归的调用 (fun(y,q)) ，其中 (y=xq^k) ， (q>p) 是更大质数。若 (kge 2) ，那么我们顺便计算 (f(y)) 。这样我们就能在 (O(1)) 的时间复杂度处理一个 (i) 了。
最后，采用一个重要的剪枝：若 (i imes maxfactor_i>n) ，显然这样的 (i) 是无需考虑的，因为不会对答案有贡献。于是在递归调用时加上这一条件即可。显然不加这一优化，我们需要遍历1到 (n) 中的所有 (i) ，时间复杂度将退化为线性。
4.通用筛法的时间复杂度证明
终于到了最具有理论价值（研究了最久）的部分了。同样地，我们将分别考虑两个算法各自的时间复杂度。为此，需要引入许多数论中的定理，尤其是素数计数相关的定理。
4.1 素数幂和相关引理
引理4：1到 (n) 中质数的个数具有估计式 (pi(n)=Theta left( {dfrac{n}{{ln n}}} ight)) ，同时隐含常数为1。
引理的证明从略。
引理5：(sumlimits_{p le n} {{p^2}} = Theta left( {dfrac{{{n^3}}}{{ln n}}} ight)) ，同时隐含常数为(dfrac 1 3)。
引理的证明从略，以上两引理可在最后参考文献中查阅。
引理6：对任意整数 (k) ，当 (n ightarrowinfty)时，(displaystyle sumlimits_{p le n} {dfrac{1}{{{p^k}}}} = int_1^n {dfrac{{mathrm{d}pi (x)}}{{{x^k}}}} ) 。
该公式称为Stieltjes积分公式，相关证明可查阅文献。
推论7： (displaystyle sumlimits_{ p > n} {dfrac{1}{{{p^2}}}} =Theta left( {dfrac{{1 }}{{nln n}}} ight))，同时隐含常数为1 。
由于并没有查到相关的公式，因此以下对该结论证明。
证明：利用引理6得 [sumlimits_{p > n} {dfrac{1}{{{p^2}}}} = int_n^infty {dfrac{{mathrm{d}pi (x)}}{{{x^2}}}} = dfrac{{pi (x)}}{{{x^2}}}igg|_n^infty + 2int_n^infty {dfrac{{mathrm{d}x}}{{{x^2}ln x}}}] 这一步是在 (n ightarrowinfty)的前提下成立的，此时分部积分中 (pi(x)) 可替换为 (dfrac{x}{ln x}) 。继续计算可以得到 [sumlimits_{p > n} {dfrac{1}{{{p^2}}}} =- dfrac{1}{{nln n}} - 2int_0^{dfrac{1}{n}} {dfrac{{mathrm{d}t}}{{ln t}}}=-dfrac{1}{{nln n}} - 2Lileft(dfrac 1 n ight)] 其中右边积分用 (t=dfrac 1 x) 替换得到， (displaystyle Li(x)=int_0^{x} {dfrac{{mathrm{d}x}}{{ln x}}}) 称为对数积分。根据 (Li(x)) 的泰勒展开式， [Li(x) = dfrac{x}{{ln x}}sumlimits_{k = 0}^infty {dfrac{{k!}}{{{{(ln x)}^k}}}}] 因而 [Lileft(dfrac 1 n ight)sim-dfrac 1 {nln n}] 故 [sumlimits_{ p > n} {dfrac{1}{{{p^2}}}} =Theta left( {dfrac{{1 }}{{nln n}}} ight)] 证毕。
4.2 算法2的时间复杂度证明
前面列举了一堆引理，现在终于到核心部分了。下面定理给出了算法2中的时间复杂度的精确估计。
定理8：计算 (S'(i)) 的时间复杂度为 (Theta left( {dfrac{{{n^{3/4}}}}{{ln n}}} ight)) ，精确到常数，其基本语句执行次数为 ( dfrac{{{32n^{3/4}}}}{{3ln n}}) 。
证明：我们将1到 (sqrt n) 中的质数分为两部分，对应算法2的流程，分别考察其对时间复杂度的贡献。
第一部分： (2le p le sqrt m)
此时对于每一个 (1le ile sqrt m) ， (b[i]) 均需要花 (Theta (1)) 的时间更新。
而对于 (p^2le ile sqrt m) ， (a[i]) 均需要花 (Theta (1)) 的时间更新。
因而对于这个 (p) ，总时间复杂度为 (Theta (2m-p^2)) 。
第二部分： (sqrt m < p le m)
此时 (a[i]) 不需要继续更新；
而对于每一个 (dfrac n ige p^2) ，即(1le ile dfrac n {p^2}) ， (b[i]) 均需要花 (Theta (1)) 的时间更新。
因而对于这个 (p) ，总时间复杂度为 (Theta left( {dfrac{n}{{{p^2}}}} ight)) 。
将两部分对 (p) 求和，得总时间复杂度为 [T(n)=sumlimits_{p le sqrt m } {(2m - {p^2})} + sumlimits_{sqrt m < p le m} {left( {dfrac{n}{{{p^2}}}} ight)} \= left(2msumlimits_{p le sqrt m } 1 ight) - left(sumlimits_{p le sqrt m } {{p^2}} ight) + left(nsumlimits_{sqrt m < p le m} {dfrac{1}{{{p^2}}}} ight)] 根据引理4、引理5和推论7，得 [T(n)=Thetaleft(2mdfrac{{sqrt m }}{{ln sqrt m }} - dfrac{{{{left( {sqrt m } ight)}^3}}}{{3ln sqrt m }} + ndfrac{1}{{sqrt m ln sqrt m }} ight) \=Thetaleft(dfrac {8msqrt m} {3ln {sqrt m}} ight)=Thetaleft(dfrac {32n^{3/4}} {3ln n} ight)] 由于上述证明使用的均为等价形式，因而时间复杂度估计使用的是记号 (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) ight| sim n ho left(dfrac 1 a ight)) ，其中 ( ho (x)) 称为迪克曼函数(Dickman function)。
性质10：对于任意 (x>1) ， ( ho (x)>0) 且随 (x) 增大迅速衰减，具有近似式 ( ho (x) approx {x^{ - x}}) 。
以上两条引理证明从略，可在最后参考文献中查阅。
上述引理表明，当 (a) 较小时， (Q) 集合中的元素个数急剧衰减；即对于一个大整数 (i)，它的最大质因子通常都比较大。
4.4 关于算法1的时间复杂度
现在回到算法1，之前已说明算法需要枚举每一个 (i) ，而通过一个剪枝 (i imes maxfactor_ile n) 可以大大减少枚举的 (i) 的数量，这一原理便可由上面指出。
然而通过一番分析，仍然不能给出较好的时间复杂度证明。最终不幸的事情还是发生了，我们有如下定理：
定理11：令 (M(n) = { i:i imes maxfacto{r_i} le n}) ，对于任意的 (alpha<1) ，有 (left| M(n) ight| = Omega ({n^alpha })) 。
证明：根据引理9， 令 (Q(n) = { i le n:maxfactor_i le {i^{1-alpha}}}) ，那么 (left| Q(n^alpha) ight| sim n^alpha ho left(dfrac 1 {1-alpha} ight)) 。对于每一个数 (xin Q(n^alpha)) ，显然 (x imes maxfacto{r_x} le n) ，这就意味着1到 (n) 中至少有 (n^alpha ho left(dfrac 1 {1-alpha} ight)) 个数在 (M(n)) 中，即 (left| M(n) ight| = Omega ({n^alpha })) 。证毕。
尽管理论十分悲观，但是不妨代入一个 (alpha=3/4) 算一下。可以查阅到 ( ho (4) approx 5 imes 10^{-3}) ，当 (n=10^{12}) 时， (n^alpha ho left(dfrac 1 {1-alpha} ight) approx 5 imes 10^6) ，仍是十分小的，这正是得益于( ho (x)) 衰减速度非常快，使得在平常使用时该算法仍具有价值。
注意到定理11给出的明显是一个下界，实际的值要大不少。下面列出了一些典型 (n) 时，实际测试结果。
当 (n=10^{11}) 时，测试结果3400万；
当 (n=10^{12}) 时，测试结果1.88亿。
可以看出，这一值甚至小于算法2的计算量。但是由于递归等原因，算法2往往跑的更快。
命题12：如果对每一质数 (pleq sqrt m) ， (i le n) 中满足 (maxfactor_i=p) 且 (ip le n) 的 (i) 平均小于 (sqrt n) 个，那么算法1就具有 (Oleft(dfrac {n^{3/4}} {log n} ight)) 的时间复杂度。
证明：设 (maxfactor_i=p) ，将(p) 分为两段分别考虑。
当 (sqrt m < pleq m) 时，算法1的时间复杂度必然小于 (displaystyle sumlimits_{sqrt m < p le m} {dfrac{n}{{{p^2}}}}) ，因为对于每个 (p) 只有不超过 (leftlfloor {dfrac{n}{{{p^2}}}} ight floor) 个 (i) 满足 (maxfactor_i=p) 。根据推论7，这一求和式不超过 (Oleft(dfrac {n^{3/4}} {log n} ight)) 。
当 (ple sqrt m) 时，根据引理5， (p) 的个数有 (Theta left( {dfrac{{{n^{1/4}}}}{{ln n}}} ight)) ，再根据命题假设，计算量不超过 (Oleft(dfrac {n^{3/4}} {log n} ight)) 。
综上可知，整个算法1的时间复杂度不超过 (Oleft(dfrac {n^{3/4}} {log n} ight)) 。证毕。
通过打表知，在(n le 10^{12}) 下命题2的假设满足，因而可以“认为”整个求积性函数前缀和算法的时间复杂度为 (Theta left( {dfrac{{{n^{3/4}}}}{{ln n}}} ight)) 。
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}}} ight)) ，即各个组成部分时间复杂度平衡。这就类似于杜教筛的 (n^{2/3}) 预处理来配平时间复杂度一样，不同之处在于，这里没有任何人为因素，完全是天然平衡的，而且是三项求和平衡。
总之，希望该算法的价值能够在日后不断体现出来，也更希望能有巨佬将算法1改进，使得理论时间复杂度真正达到 (Theta left( {dfrac{{{n^{3/4}}}}{{ln n}}} ight)) !
最后，由于本人才疏学浅，文章中如有任何不当之处，欢迎批评指正，谢谢！