一种可行的计算质数数量函数的实现方法

9

有哪些计算上可行的伪代码可以实现素数计数函数

我最初尝试编写Hardy-Wright算法,但它的阶乘开始产生可怜的溢出,而其他许多算法似乎也会遇到类似的问题。我在谷歌上搜寻了实用的解决方案,但最多只找到了非常晦涩的数学方法,我从未见过这些方法在传统程序中得到实现。


1
“floor(x/j) * j == x - (x%j)” 这个式子是不是正确的呢?那么你提供的公式就变成了 “pi(x) = (-1) + SUM{j=3..n}( ((j-2)!) % j )” (?). 接下来使用模乘法(即 5! % 7 == (((((2*3)%7)*4)%7)*5)%7)。 - Will Ness
@SeverynKozak 他们说这不是一个编程问题的原因是因为你的问题中没有包含任何代码。 - brian beuning
除了伪代码之外,我在另一个SO线程中的回答提供了JavaScript中基本的Legarias-Miller-Odzlinko(LMO)1985算法/方法的可运行代码片段(可在浏览器中运行)。 LMO是所有现代实用通用素数计数函数的基础,而变体(Gourdon,2001)已被用于计算高达1e28的质数数量,尽管即使在现代计算机上使用多线程和大量内存,这也需要数月时间。 - GordonBGood
2个回答

22

素数计数函数pi(x)计算不超过x的质数数量,几个世纪以来一直吸引着数学家们。在18世纪初,阿德里安-玛丽·勒让德使用辅助函数phi(x,a)给出了一个公式,该函数计算不大于x且未被前a个质数筛去的数字数量;例如,phi(50,3) = 14表示数字1、7、11、13、17、19、23、29、31、37、41、43、47和49。可以将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的所有质数,通过筛法计算。缓存很重要,没有缓存,运行时间将呈指数级增长。根据欧拉函数phi和Legendre的计数公式,pi(x) = phi(x,a) + a - 1,其中a = pi(floor(sqrt(x)))。Legendre使用他的公式计算了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)限制在nlogn和n(logn + loglogn)之间,对于n > 5,因此在边界处计算pi并使用二分法确定第n个质数,在边界接近时切换到筛法。
我在我的博客的几篇文章中讨论了质数。

2
最近,已经计算出了pi(10^25)和pi(10^26)。请参见此处第40页 - qwr
1
你的帖子对我的实现非常有帮助。同时,记忆化pi(x)也是一个很大的加速。我在计算a、b、c时遇到了一些四舍五入问题,所以我使用了danaj的四舍五入方法:https://programmingpraxis.com/2011/07/22/counting-primes-using-legendres-formula/#comment-5958 - qwr

2
维基百科也可能有所帮助。素数计数的文章包含了一些指针。对于初学者,我建议阅读"用于评估π(x)的算法"部分中Meissel算法,这是一个最简单的不生成所有质数的算法。
我还发现Pomerance和Crandall的书"Prime numbers a computational perspective"很有帮助。这本书详细而易懂地描述了素数计数方法。但请记住,由于其性质,该主题对大多数读者来说有点太高级了。

网页内容由stack overflow 提供, 点击上面的
可以查看英文原文,
原文链接