是否有快速、功能强大的素数生成器？答案

【问题标题】：Is there a fast, functional prime generator?是否有快速、功能强大的素数生成器？
【发布时间】：2017-06-26 09:18:57
【问题描述】：

假设我有一个自然数 n，我想要一个列表（或其他），其中包含直到 n 的所有素数。

经典的素筛算法在O(n log n) 时间和O(n) 空间中运行——它适用于更多命令式语言，但从根本上需要对列表和随机访问进行就地修改。

有一个涉及优先级队列的功能版本，非常漂亮——您可以查看here。这在大约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) 纯算法，我仍然会很高兴。

【问题讨论】：

链接的文章提出了 Theta(n log n log log n)-time 算法，一个 log n 因子关闭 Eratosthenes 的标准筛，使用功能优先级队列将写入合并到数组.
虽然我相信这是这里的主题，但如果您在一段时间后没有得到任何答案，您可能还想查看 computerscience.stackexchange。
@DavidEisenstat 那么问题可能是我不能很好地计算复杂性。我确实很困惑解决这个问题。
不管怎样，this 声称是文章中算法的 Haskell 实现。
请参阅here 以获得非常好的埃拉托色尼 Clojure 筛子。

标签： algorithm functional-programming primes sieve-of-eratosthenes

【解决方案1】：

这是 Melissa O'Neill 算法的 Haskell 实现（来自链接的文章）。与 Gassa 所链接的实现不同，我对惰性的使用最少，因此性能分析很清晰 -- O(n log n log log n)，即 n log log 中的线性n，Eratosthenes 的命令式 Sieve 进行的写入次数。

堆实现只是一个锦标赛树。平衡逻辑在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)

【讨论】：

公平地说，性能分析并不完全是微不足道的。必须考虑有多少弹出/推送操作，以及操作发生时堆有多大。如果您忽略它（假设它始终是最大大小 ~n/ln(n)），那么您会大大高估时间复杂度（我在上面发布的错误估计）。
@RichardRast Pops 是每个 Eratosthenes 写入一个，即 Theta(n log log n)。在所检查的窗口中，推送是每个素数之一。两者都需要时间 log k，其中 k 是堆的当前大小；由于 log √n = (log n)/2 = Theta(log n)，绝大多数操作都花费 Theta(log n)，因为堆在考虑的前 √n log n 个数中增长到 √n 个素数。
事实证明，我真正的问题不是算法本身，而是计算其复杂性。但这很有帮助，我会将其标记为已接受，因为我认为它将继续对其他人有用。
@WillNess Haskell 素数生成器传统上提供所有素数的惰性列表。最小化堆的大小会增加复杂性。

【解决方案2】：

就是这样，如果（Haskell 的）纯数组算作纯数组（它们应该，IMO）。复杂性显然是 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;等等。 See also.

另外，Python generator-based sieve 也可能被认为是功能性的。 Python 的dicts 的性能非常好，empirically，虽然我不确定那里使用的用于避免重复合成的倍数过度生产方案的确切成本。

更新：testing it 确实产生了有利的结果，正如预期的那样：

{-     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
-}

empirical orders of growth 计算产生 n 个素数，其中 O(n log log n) 通常被视为 n^{1.05... 1.10} 和 O(n log n log log n) 为 n^1.20...1.25。 p>

"nested-feed" 变体实现了postponement 技术（也可以在上面链接的 Python answer 中看到）实现堆大小的二次减小，这显然对经验复杂性有显着影响，即使不是完全在这个答案的基于数组的代码中达到了更好的结果，在 ideone.com 上 is able to produce 10 秒内有 1000 万个素数（总体增长率仅为 n^1.09 在测试范围内）。

（"original heap"当然是the other answer这里的代码）。

【讨论】：

【解决方案3】：

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

这里是：

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

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

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

【讨论】：

It returns 0 for all composites and 1 for all primes (including 2) 听起来像 primality test 而不是 generating function。还是……
@greybeard 只需将整个函数乘以 x，然后它就是素数生成函数！ ? 我也想我可能搞砸了，它可能会返回 0 对于所有素数而不是 1。我不记得自从我研究这个方程以来已经有好几个月了。
@AlbertRenshaw 不，生成函数将为每个 n 生成 n-th prime，中间没有零。说到阶乘，一个简单的素数测试函数是基于Wilson's theorem：isPrime(n): (product([1..n-1]) % n) == (n-1)。
@WillNess 哦，我明白了。是的，除非你在无限和中工作并且有一些花哨的模技巧，否则我无法用我的函数来做到这一点。不过，这将是无限迭代，不会很快。
它是如何工作的？ (4Г(x)+4)/x 通常不是整数，(-1)^not-an-integer 定义不明确。