9

任何人都可以提供任何素数计数函数实现的计算上可行的伪代码吗？我最初尝试编写Hardy-Wright 算法，但它的阶乘开始产生严重的溢出，而且许多其他算法似乎必然会产生类似的问题。我在 Google 上搜索了实用的解决方案，但充其量只是发现了非常深奥的数学，我从未见过在传统程序中实现过这些数学。

4

2 回答 2

19

素数计数函数 pi(x) 计算不超过 x 的素数的数量，几个世纪以来一直让数学家着迷。18 世纪初，Adrien-Marie Legendre 给出了一个公式，它使用辅助函数 phi(x,a) 计算不大于 x 且未被第一个 a 素数筛选的数；例如，对于数字 1、7、11、13、17、19、23、29、31、37、41、43、47 和 49，phi(50,3) = 14。 phi 函数可以计算为 phi (x,a) = phi(x,a-1) - phi(x/P(a),a-1)，其中 phi(x,1) 是不超过 x 和 P(a) 的奇数个数是第 a 个素数（从 P(1)=2 开始计数）。

function phi(x, a)
  if (phi(x, a) is in cache)
    return phi(x, a) from cache
  if (a == 1)
    return (x + 1) // 2
  t := phi(x, a-1) - phi(x // p[a], a-1)
  insert phi(x, a) = t in cache
  return t

数组 p 存储小 a 的第 a 个素数，通过筛分计算得出。缓存很重要；没有它，运行时间将是指数级的。给定 phi，Legendre 的素数计数公式是 pi(x) = phi(x,a) + a - 1，其中 a = pi(floor(sqrt(x)))。勒让德使用他的公式计算 pi(10^6)，但他报告了 78526 而不是正确答案 78498，即使是错误的，对于复杂的手动计算来说也非常接近。

在 1950 年代，Derrick H. Lehmer 给出了一种改进的素数计数算法：

function pi(x)
  if (x < limit) return count(primes(x))
  a := pi(root(x, 4)) # fourth root of x
  b := pi(root(x, 2)) # square root of x
  c := pi(root(x, 3)) # cube root of x
  sum := phi(x,a) + (b+a-2) * (b-a+1) / 2
  for i from a+1 to b
    w := x / p[i]
    lim := pi(sqrt(w))
    sum := sum - pi(w)
    if (i <= c)
      for j from i to lim
        sum := sum - pi(w / p[j]) + j - 1
  return sum

例如，pi(10^12) = 37607912018。即使使用这些算法及其现代变体以及速度非常快的计算机，计算大的 pi 值仍然非常繁琐；在撰写本文时，已知的最大值为 pi(10^24) = 18435599767349200867866。

要使用此算法计算第 n 个素数，素数定理的推论将第 n 个素数 P(n) 限制在 n log n 和 n(log n + log log n) 之间，n > 5，因此计算pi 在边界处并使用二等分来确定第 n 个素数，当边界接近时切换到筛分。

我在博客的几个条目中讨论了质数。

于 2013-09-28T23:14:13.397 回答
3

维基百科也可能有所帮助。关于素数计数的文章包含一些提示。对于初学者，我会推荐 Meissel 在“评估 π(x) 的算法”一节中的算法，这是最简单的算法之一，不会生成所有素数。

我还发现 Pomerance 和 Crandall 的书“素数计算视角”很有帮助。这本书对素数计数方法进行了详细且易于理解的描述。但请记住，对于这里的大多数读者来说，该主题的性质有点太高级了。

于 2013-09-28T20:20:03.303 回答