其实求 \(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