关于一种积性函数前缀和的通用筛法的时间复杂度证明

4月11日更新:修正了一些错误。同时,由于知乎排版功能比较欠缺,这里附上博客版本链接,该版本可能更有助于看懂本文的层次结构;同时该版本中加入了一些模板题。链接为:关于一种积性函数前缀和的通用筛法的时间复杂度证明 - zbh2047 - 博客园

写在前面


自从学了杜教筛,突然对这一类问题产生了浓厚的兴趣。于是花了好几天时间,期间也学习了数论相关的许多理论,总算是努力没有白费,把这一算法的时间复杂度分析出来了。由于网上并没有这一筛法的相关介绍(仅有的一篇文献已在附录列出),在此写一篇详尽的分析,也便于以后忘记时查阅。

为了叙述方便,引入以下记号:

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


通用筛法原理

算法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 {\frac{n}{x}} \right\rfloor } \right):1 \le x \le m\} \]

证明: \sum\limits_{F < p,ip \le n} {f(p)} =S'\left( {\left\lfloor {\frac{n}{i}} \right\rfloor } \right)-S'(F)

对于 S'\left( {\left\lfloor {\frac{n}{i}} \right\rfloor } \right) ,若 i\le m ,则属于 S_2 ,否则属于 S_1

对于 S'(F) ,由 F^2\le ip\le nF\le m ,属于 S_1

对于最后还需要加上的全体质数,显然 S'(n) 属于 S_2

证毕。

算法2

对于求 S'(i) ,采用欧拉筛法的思想,初始时将每个数都视为质数,将 f(j),2 \le j \le i 全部加到 S'(i) 中。然后从小到大遍历每一质数 p ,令 in 开始,从大到小执行以下操作:

\[S'(i) - = {\rm{ }}\left( {S'\left( {\left\lfloor {\frac{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 后,对于每一个 iS'(i) 均扣除且仅扣除了所有的 f(x) ,其中 x 是合数且maxfactor_x \le p

初始时,显然是正确的。考虑当前的质数 p

\[{\rm{ }}\left( {S'\left( {\left\lfloor {\frac{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中才有用。


通用筛法实现

根据上述原理,这里仍分两部分考虑实现。

首先考虑计算 S'(i) 的部分。根据结论2,我们只需要两个大小均为 m 的数组,其中 a[i]=S'(i)b[i]=S'\left( {\left\lfloor {\frac{n}{i}} \right\rfloor } \right) ,这只需要 \[\Theta (m)\] 的空间。

对于筛法的公式,不难发现以下三种情况:

(1)若 \[\frac{i}{p} > \left\lfloor {\sqrt n } \right\rfloor \] ,那么 S'(i) 便对应了 b[n/i],而S'\left( {\left\lfloor {\frac{i}{p}} \right\rfloor } \right) 对应了 b[n/i*p]

(2)若 \[i > \left\lfloor {\sqrt n } \right\rfloor \]\[\frac{i}{p} \le \left\lfloor {\sqrt n } \right\rfloor \] ,那么 S'(i) 便对应了 b[n/i],而 S'\left( {\left\lfloor {\frac{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 {\frac{i}{p}} \right\rfloor } \right) 便对应了 a[i/p]

由此可见,更新 a,b 数组不会用到 a,b 数组之外的元素 S'(i)

计算完 S'(i) 后,根据之前的算法,我们要枚举每一整数 i ,计算 f(i)\left(S'\left( {\left\lfloor {\frac{n}{i}} \right\rfloor } \right)-S'(maxfactor_i)\right) 。为了快速枚全部 i ,可以考虑递归实现。

设函数形式为 fun(x,p) ,表示 x 的最大质因子是 p ,那么我们先计算上面的式子,然后转移到更大的 x ,即递归的调用 fun(y,q) ,其中 y=xq^kq>p 是更大质数。若 k\ge 2 ,那么我们顺便计算 f(y) 。这样我们就能在 O(1) 的时间复杂度处理一个 i 了。

最后,采用一个重要的剪枝:若 i \times maxfactor_i>n ,显然这样的 i 是无需考虑的,因为不会对答案有贡献。于是在递归调用时加上这一条件即可。显然不加这一优化,我们需要遍历1到 n 中的所有 i ,时间复杂度将退化为线性。


通用筛法的时间复杂度证明

终于到了最具有理论价值(研究了最久)的部分了。同样地,我们将分别考虑两个算法各自的时间复杂度。为此,需要引入许多数论中的定理,尤其是素数计数相关的定理。

引理4:1到 n 中质数的个数具有估计式 \pi(n)=\Theta \left( {\frac{n}{{\ln n}}} \right) ,同时隐含常数为1。

引理的证明从略。

引理5:\[\sum\limits_{p \le n} {{p^2}} = \Theta \left( {\frac{{{n^3}}}{{3\ln n}}} \right)\]

引理的证明从略,以上两引理可在最后参考文献中查阅。

引理6:对任意整数 k ,当 n\rightarrow\infty时,\[\sum\limits_{p \le n} {\frac{1}{{{p^k}}}} = \int_1^n {\frac{{\mathrm{d}\pi (x)}}{{{x^k}}}} \]

该公式称为Stieltjes积分公式,相关证明可查阅文献。

推论7: \[\sum\limits_{ p > n} {\frac{1}{{{p^2}}}} =\Theta \left( {\frac{{1 }}{{n\ln n}}} \right)\]

由于并没有查到相关的公式,因此以下对该结论证明。

证明: \sum\limits_{p > n} {\frac{1}{{{p^2}}}} = \int_n^\infty {\frac{{\mathrm{d}\pi (x)}}{{{x^2}}}} = \frac{{\pi (x)}}{{{x^2}}}\bigg|_n^\infty + 2\int_n^\infty {\frac{{\mathrm{d}x}}{{{x^2}\ln x}}}

这一步是在 n\rightarrow\infty的前提下成立的,此时分部积分中 \pi(x) 可替换为 \frac{x}{\ln x}

继续计算可以得到 \sum\limits_{p > n} {\frac{1}{{{p^2}}}} =- \frac{1}{{n\ln n}} - 2\int_0^{\frac{1}{n}} {\frac{{\mathrm{d}t}}{{\ln t}}}=-\frac{1}{{n\ln n}} - 2Li\left(\frac 1 n\right)

其中右边积分用 t=\frac 1 x 替换得到, Li(x)=\int_0^{x} {\frac{{\mathrm{d}x}}{{\ln x}}} 称为对数积分。

根据 Li(x) 的泰勒展开式, Li(x) = \frac{x}{{\ln x}}\sum\limits_{k = 0}^\infty {\frac{{k!}}{{{{(\ln x)}^k}}}}

因而 Li\left(\frac 1 n\right)\sim-\frac 1 {n\ln n}

\[\sum\limits_{ p > n} {\frac{1}{{{p^2}}}} =\Theta \left( {\frac{{1 }}{{n\ln n}}} \right)\] ,证毕。


前面列举了一堆引理,现在终于到核心部分了。下面定理给出了算法2中的时间复杂度的精确估计。

定理8:计算 S'(i) 的时间复杂度为 \Theta \left( {\frac{{{n^{3/4}}}}{{\ln n}}} \right) ,精确到常数,其基本语句执行次数为 \Theta \left( {\frac{{{32n^{3/4}}}}{{3\ln n}}} \right)

证明:我们将1到 \sqrt n 中的质数分为两部分,对应算法2的流程,分别考察其对时间复杂度的贡献。

第一部分: 2\le p \le \sqrt m

此时对于每一个 1\le i\le \sqrt mb[i] 均需要花 \Theta (1) 的时间更新。

而对于 p^2\le i\le \sqrt ma[i] 均需要花 \Theta (1) 的时间更新。

因而对于这个 p ,总时间复杂度为 \Theta (2m-p^2)

第二部分: \sqrt m < p \le m

此时 a[i] 不需要继续更新;

而对于每一个 \frac n i\ge p^2 ,即1\le i\le \frac n {p^2}b[i] 均需要花 \Theta (1) 的时间更新。

因而对于这个 p ,总时间复杂度为 \Theta \left( {\frac{n}{{{p^2}}}} \right)

将两部分对 p 求和,得总时间复杂度为

T(n)=\sum\limits_{p \le \sqrt m } {(2m - {p^2})} + \sum\limits_{\sqrt m < p \le m} {\left( {\frac{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} {\frac{1}{{{p^2}}}} \right)

根据引理4、引理5和推论7,得

T(n)=\Theta\left(2m\frac{{\sqrt m }}{{\ln \sqrt m }} - \frac{{{{\left( {\sqrt m } \right)}^3}}}{{3\ln \sqrt m }} + n\frac{1}{{\sqrt m \ln \sqrt m }}\right) \\=\Theta\left(\frac {8m\sqrt m} {3\ln {\sqrt m}} \right)=\Theta\left(\frac {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毫秒左右。


接下来考虑算法1。为此,仍需要引入一些引理。

引理9:对于给定的实常数 0<a<1 ,令 Q(n) = \{ i \le n:maxfactor_i \le {i^a}\} ,那么 \left| Q(n) \right| \sim n\rho \left(\frac 1 a\right) ,其中 \rho (x) 称为迪克曼函数(Dickman function)。

性质10:对于任意 x>1\rho (x)>0 且随 x 增大迅速衰减,具有近似式 \rho (x) \approx {x^{ - x}}

以上两条引理证明从略,可在最后参考文献中查阅。

上述引理表明,当 a 较小时, Q 集合中的元素个数急剧衰减;即对于一个大整数 i,它的最大质因子通常都比较大。

现在回到算法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(\frac 1 {1-\alpha}\right) 。对于每一个数 x\in Q(n^\alpha) ,显然 x \times maxfacto{r_x} \le n ,这就意味着1到 n 中至少有 n^\alpha\rho \left(\frac 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(\frac 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 mi \le n 中满足 maxfactor_i=p ip \le n i 小于 \sqrt n 个,那么算法1就具有 O\left(\frac {n^{3/4}} {\log n} \right) 的时间复杂度。

证明:设 maxfactor_i=p ,将p 分为两段分别考虑。

\sqrt m <p\leq m 时,算法1的时间复杂度必然小于 \sum\limits_{\sqrt m < p \le m} {\frac{n}{{{p^2}}}} ,因为对于每个 p 只有不超过 \left\lfloor {\frac{n}{{{p^2}}}} \right\rfloor i 满足 maxfactor_i=p 。根据推论7,这一求和式不超过 O\left(\frac {n^{3/4}} {\log n} \right)

p\le \sqrt m 时,根据引理5, p 的个数有 \Theta \left( {\frac{{{n^{1/4}}}}{{\ln n}}} \right) ,再根据命题假设,计算量不超过 O\left(\frac {n^{3/4}} {\log n} \right)

综上可知,整个算法1的时间复杂度不超过 O\left(\frac {n^{3/4}} {\log n} \right) 。证毕。

通过打表知,在n \le 10^{12} 下命题2的假设满足,因而可以“认为”整个求积性函数前缀和算法的时间复杂度为 \Theta \left( {\frac{{{n^{3/4}}}}{{\ln n}}} \right)


后记

尽管最后的结果表明算法时间复杂度似乎是线性的,(即对于任意 \alpha<1 ,该算法时间复杂度均高于 n^\alpha ),这已经将算法的理论价值大打折扣,但是:一方面,该算法仍具有巨大的应用价值。例如,用该算法求莫比乌斯函数前缀和,在oj上的提交结果为 n=10^{10} 下运行350ms,已接近杜教筛的运行时间;另一方面,算法2具有严格的理论时间复杂度证明(精确到常数),从而具有很强的理论价值。例如,可用来求1到 n 中的质数个数,或1到 n 中的全部质数之和等。

个人感觉,算法2是一个非常优美的算法,因为其时间复杂度的证明过程中,恰好最后的三项均为 \Theta \left( {\frac{{{n^{3/4}}}}{{\ln n}}} \right) ,即各个组成部分时间复杂度平衡。这就类似于杜教筛的 n^{2/3} 预处理来配平时间复杂度一样,不同之处在于,这里没有任何人为因素,完全是天然平衡的,而且是三项求和平衡。

总之,希望该算法的价值能够在日后不断体现出来,也更希望能有巨佬将算法1改进,使得理论时间复杂度真正达到 \Theta \left( {\frac{{{n^{3/4}}}}{{\ln n}}} \right) !

最后,由于本人才疏学浅,文章中如有任何不当之处,欢迎批评指正,谢谢!


参考文献

1.Prime Sums -- from Wolfram MathWorld

2.Dickman Function

3.SPOJ.com - Problem TEES

编辑于 2018-04-11