有快速、实用的素数生成器吗？

Question

假设我有一个自然数

n

并且我想要一个包含直到

n

的所有素数的列表（或其他）。

经典的素数筛算法在

O(n log n)

时间和

O(n)

空间中运行——对于更多命令式语言来说很好，但需要从根本上对列表和随机访问进行就地修改。

有一个涉及优先级队列的功能版本，它非常灵活——您可以在这里查看它。这具有更好的空间复杂度，约为

O(n / log(n))

（渐近更好，但在实际规模上存在争议）。不幸的是，时间分析很糟糕，但非常接近

O(n^2)

（实际上，我认为这是关于

O(n log(n) Li(n))

，但

log(n) Li(n)

大约是

n

）。

渐近地讲，实际上最好在生成每个数字时使用连续的试除法来检查每个数字的素数，因为这只需要

O(1)

空间和

O(n^{3/2})

时间。有更好的办法吗？

编辑：事实证明我的计算完全不正确。文章中的算法是

O(n (log n) (log log n))

，文章对此进行了解释和证明（并参见下面的答案），而不是我上面提出的复杂混乱。如果有的话，我仍然很高兴看到一个真正的

O(n log log n)

纯算法。

Answer 1

这是 Melissa O'Neill 算法的 Haskell 实现（来自链接的文章）。与 Gassa 链接的实现不同，我最大限度地减少了惰性，因此性能分析很清晰 - O(n log n log log n)，即 n log log n 中的线性对数，由埃拉托斯特尼的命令式筛法所写。

堆实现只是一个锦标赛树。平衡逻辑在

push

；通过每次交换子树，我们确保对于每个分支，左子树的大小相同或比右子树大一，这确保了深度 O(log n)。

module Sieve where
 
type Nat = Int
 
data Heap = Leaf !Nat !Nat
          | Branch !Nat !Heap !Heap
          deriving Show
 
top :: Heap -> Nat
top (Leaf n _) = n
top (Branch n _ _) = n
 
leaf :: Nat -> Heap
leaf p = Leaf (3 * p) p
 
branch :: Heap -> Heap -> Heap
branch h1 h2 = Branch (min (top h1) (top h2)) h1 h2
 
pop :: Heap -> Heap
pop (Leaf n p) = Leaf (n + 2 * p) p
pop (Branch _ h1 h2)
  = case compare (top h1) (top h2) of
        LT -> branch (pop h1) h2
        EQ -> branch (pop h1) (pop h2)
        GT -> branch h1 (pop h2)
 
push :: Nat -> Heap -> Heap
push p h@(Leaf _ _) = branch (leaf p) h
push p (Branch _ h1 h2) = branch (push p h2) h1
 
primes :: [Nat]
primes
  = let helper n h
          = case compare n (top h) of
                LT -> n : helper (n + 2) (push n h)
                EQ -> helper (n + 2) (pop h)
                GT -> helper n (pop h)
      in 2 : 3 : helper 5 (leaf 3)

Answer 2

这就是，如果（Haskell 的）纯数组算作纯数组（在我看来，它们应该算作纯数组）。复杂度显然是 O(n log (log n))，前提是

accumArray

确实为给出的每个条目花费 O(1) 时间，因为它应该这样做：

import Data.Array.Unboxed 
import Data.List (tails, inits)

primes = 2 : [ n |
   (r:q:_, px) <- zip (tails (2 : [p^2 | p <- primes]))
                      (inits primes),
   (n,True)    <- assocs ( accumArray (\_ _ -> False) True
                             (r+1,q-1)
                             [ (m,()) | p <- px
                                      , let s = div (r+p) p * p
                                      , m <- [s,s+p..q-1]]
                             :: UArray Int Bool ) ]

通过素数的连续平方之间的段计算素数，通过枚举素数列表相应前缀的倍数（使用

inits

）来生成复合数，就像埃拉托色尼的任何适当的筛法一样，通过重复相加。

因此，素数 {2,3} 用于筛选从 10 到 24 的线段； {2,3,5}从26到48；等等。 另请参阅。

此外，Python 基于生成器的筛子也可能被认为是功能性的。 Python 的

dict

的性能非常好，根据经验，尽管我不确定用于避免重复复合的倍数过度生产方案的确切成本。

---

更新：测试它确实产生了良好的结果，正如预期的那样：

{-     original heap       tweaked           nested-feed         array-based
          (3*p,p)         (p*p,2*p)            JBwoVL              abPSOx
          6Uv0cL          2x speed-up     another 3x+ speed-up
                n^                n^                  n^                  n^
100K:  0.78s             0.38s               0.13s              0.065s    
200K:  2.02s   1.37      0.97s   1.35        0.29s   1.16       0.13s    1.00
400K:  5.05s   1.32      2.40s   1.31        0.70s   1.27       0.29s    1.16
800K: 12.37s   1.29                     1M:  2.10s   1.20       0.82s    1.13
 2M:                                                            1.71s    1.06
 4M:                                                            3.72s    1.12
10M:                                                            9.84s    1.06
    overall in the tested range:
               1.33                                  1.21                1.09
-}

使用经验增长阶计算产生n素数，其中O(n log log n)通常被视为n^1.05...1.10和O(n log n log log n) ) 为 n^1.20...1.25.

“nested-feed”变体实现了postponement技术（也可以在上面链接的Pythonanswer中看到），该技术实现了堆大小的二次减少，这显然对经验复杂性有显着的影响，即使对于这个答案的基于数组的代码来说，还没有完全达到更好的结果，它能够在 ideone.com 上在 10 秒内产生 1000 万个素数（总体增长率仅为 n1.09）测试范围）。

（

“原始堆”当然是来自另一个答案的代码）。

Answer 3

不久前我导出了一个素数生成函数（按顺序生成所有素数）。还为其创建了 6 页的校样。我认为它实际上是历史上第一个素数生成函数（至少我找不到任何其他例子）。

这里是：

(-1)^((4*gamma(x)+4)/x)-1

不确定计算速度有多快。它对所有素数返回 0（或者可能是 1，记不清了）。伽玛函数本质上是阶乘，因此可以在早期快速运行。不过，将负 1 提高到小数指数完全是另一回事，我相信它可能使用 base_e 中的积分，或者可能使用一些三角函数；不记得了。

我不懂 LaTeX，所以如果有人想编辑我的帖子并包含 LaTeX 版本，那就太棒了！