其实求\(x\)以内素数个数,有时间复杂度为\(O(x^{2/3}loglogx)\) ,空间复杂度为\(O(x^{1/3})\)优秀算法

\(prime(i)\)表示第\(i\)个素数。

\(\pi(x)\)表示小于等于\(x\)的素数个数。

\(\phi(x, a)\)表示小于等于\(x\)的数中,不被前\(a\)个素数整除的数的个数,特别地,\(\phi(x, 0) = x\)

\(\phi(x, a)\)分类,其中\(P_k(x, a)\)表示小于等于\(x\)的数中,不被前\(a\)个素数整除且有且仅有\(k\)个素数因子的数的个数,特别地,\(P_0(x, a) = 1\)。那么有,\(\phi(x,a) = P_0(x, a) + P_1(x, a) + \dots + P_k(x, a) + \dots\)

考虑到,

\[\pi(x) = P_1(x, a) + a = \phi(x, a) - P_0(x, a) - \sum_{k = 2}^{\infty}P_k(x, a) + a = \phi(x, a) -1 - \sum_{k = 2}^{\infty}P_k(x, a) + a\]

所以,对于特定的\(a\),我们只需要知道\(\phi(x, a)\)\(\sum_{k = 2}^{\infty}P_k(x, a)\),就可以计算出\(\pi(x)\)

\(a = \pi(x ^ {1 / 3})\),那么\(P_3(x, a) = P_4(x, a) = \dots = 0\),也就是说我们只需要计算\(\phi(x, a)\)\(P_2(x, a)\)

先考虑如何计算\(\phi(x, a)\)

注意到\(\phi(x, a - 1)\)\(\phi(x, a)\)多了那些不被前\(a - 1\)个素数整除但可以被\(prime(a)\)整除的数,即多了\(\phi(\frac{x}{prime(a)}, a - 1)\)。因此,\(\phi(x, a) = \phi(x, a - 1) - \phi(\frac{x}{prime(a)}, a - 1)\)

我们解这个递归式就可以计算出\(\phi(x, a)\)

假设我们能够解这个递归式,那么我们如何计算时间复杂度呢?

因为这个递归树是一棵满二叉树,所以叶子结点个数\(n_0\)与非叶结点个数\(n_2\)满足\(n_0 = n_2 + 1\),也就是说时间复杂度为\(O(n_0 + n_2) = O(2n_0 - 1) = O(n_0)\),因此我们只需要考虑递归树的叶子结点就好了。

因为递归树上的任一结点都可以被唯一地表示为\(\phi(\frac{x}{n}, b)\),那么有一个显然的递归边界如下:

  1. \(b = 0\)\(\phi(\frac{x}{n}, b) = \frac{x}{n}\)
  2. \(\frac{x}{n} \leq prime(b)\)\(\phi(\frac{x}{n}, b) = 1\)

不过可以证明,这样的递归边界会使时间复杂度至少为\(\Omega(\frac{x}{log^3x})\)。考虑\(n\)\(3\)个小于等于\(x^{1/3}\)的素数的乘积,由于\(prime(b) \leq x^{1/3}\),因此必定存在叶子结点\(\phi(\frac{x}{n}, b)\)。根据素数定理,小于等于\(x^{1/3}\)的素数大概有\(\frac{x^{1/3}}{logx^{1/3}}\)个,故递归树至少有\(\Omega(\frac{x}{log^3x})\)个叶子结点,达不到之前所说\(O(x^{2/3}loglogx)\)这样的时间复杂度。

现在考虑修改递归边界如下:

  1. \(b = 0\)并且\(n \leq x^{1/3}\)
  2. \(n > x^{1/3}\)

再考虑如何计算\(P_2(x, a)\)

根据定义,\(P_2(x, a) = \sum_{b = a + 1}^{\pi(x^{1/2})}(\pi(\frac{x}{prime(b)}) - b + 1)\)。也就是说,我们只需要求出\(\pi(x^{1/3})\)\(\pi(x^{1/2})\)\(\sum_{b = a + 1}^{\pi(x^{1/2})}\pi({\frac{x}{prime(b)}})\)就可以求出\(P_2(x, a)\)

一个简单的想法是这样的,因为\(\frac{x}{prime(b)} \in [x^{1/2}, x^{2/3}]\),所以我们可以通过欧拉筛\(O(x^{2/3})\)的时间复杂度和空间复杂度下,预处理出\(\{prime(i) | 1 \leq i \leq x^{2/3} \}\)\(\{\pi(i) | 1 \leq i \leq x^{2/3} \}\),从而计算出\(P_2(x, a)\)

实际上还存在着时间复杂度\(O(x^{2/3}loglogx)\),空间复杂度\(O(x^{1/3})\)的做法。

主要的思想是将\([1, x^{2/3}]\)区间分块,每块长度\(len = x^{1/3}\),第\(j\)块表示区间\(B_j = [(j - 1)len + 1, jlen]\),那么最多有\(x^{1/3} + 1\)块,注意最后一块的长度可能小于\(x^{1/3}\)

在第\(1\)轮,我们通过欧拉筛预处理出\(P = \{prime(i)|prime(i) \leq x^{1/3} \}\),这样同时我们也知道了\(\pi(x^{1/3})\)

在第\(j\)轮,我们进行两个操作。

操作一。我们首先记录在第\(j - 1\)轮中所求得的\(\pi((j - 1)len)\),然后运用区间筛处理出第\(j\)块的素数,同时再根据\(\pi((j - 1)len)\)计算出\(\{\pi(y)|(j - 1)len \leq y \leq jlen \}\)。注意,当\((j - 1)len + 1 \leq x^{1/2} < jlen + 1\)时,我们可以确定\(\pi(x^{1/2})\)

操作二。考虑若素数\(p(x^{1/3} < p \leq x^{1/2})\)使得\((j - 1)len + 1 \leq \frac{x}{p} < jlen + 1\),那么设\(I_j = (\frac{x}{jlen + 1}, \frac{x}{(j - 1)len + 1}] \cap (x^{1/3}, x^{1/2}]\),我们再次运用区间筛处理出\(I_j\)中的素数,然后统计\(\pi(\frac{x}{p})(p \in I_j)\)的和即可。注意,要满足\(I_j \neq \emptyset\)\(\frac{x}{jlen + 1} < x^{1/2}\),即\(jlen + 1 > x ^{1/2}\)

每轮进行操作一的时间复杂度为\(O(x^{1/3}loglogx)​\),最多进行\(x^{1/3}​\)轮,故操作一的总时间复杂度为\(O(x^{2/3}loglogx)​\)。每轮进行操作一的空间复杂度为\(O(x^{1/3})​\),由于每一轮只保存当前的\(B_j​\),故操作一的总空间复杂度为\(O(x^{1/3})​\)。由于\(\frac{x}{(j - 1)len + 1} - \frac{x}{jlen + 1} + 1 \leq \frac{xlen}{((j - 1)len + 1)^2} \leq 2x^{1/3}​\),所以每轮进行操作二的时间复杂度为\(O(x^{1/3}loglogx)​\),最多进行\(x^{1/3}​\)轮,故操作二的总时间复杂度为\(O(x^{2/3}loglogx)​\)。每轮进行操作二的空间复杂度为\(O(x^{1/3})​\),由于每一轮只保存当前的\(I_j​\),故操作二的总空间复杂度为\(O(x^{1/3})​\)

综上所述,我们得到了一个求\(x\)以内素数个数,时间复杂度为\(O(x^{2/3}loglogx)\) ,空间复杂度为\(O(x^{1/3})\)的算法。

代码大概是这个样子的:

练习题:Count primes

其实还有更优秀的算法(逃

comments powered by Disqus