在给定所有先前的情况下找到下一个素数

Question

我正在编写一个递归无限素数生成器，我几乎可以肯定我可以更好地优化它。

现在，除了前十几个素数的查找表之外，对递归函数的每次调用都会收到所有先前素数的列表。

因为它是一个惰性生成器，所以现在我只是过滤掉之前任何素数的模 0 的任何数字，并获取第一个未过滤的结果。 （我正在使用短路的检查，因此前一个素数第一次将当前数字平均除以该信息会中止。）

现在，在搜索第 400 个素数 (37,813) 时，我的性能会下降。我正在寻找方法来使用我拥有所有先验素数列表的独特事实，并且只搜索下一个，以改进我的过滤算法。 （我能找到的大多数信息都提供了非惰性筛子来查找限制下的素数，或者在给定 p_n-1 的情况下找到第 p_n 个素数的方法，而不是优化给定 2...p_n-1 个素数，求 p_n。）

例如，我知道第 p_n 个素数必须位于 (p_n-1 + 1)...(p_{n -1}+p_n-2)。现在我开始在 p_n-1 + 2 处过滤整数（因为 p_n-1 + 1 只能是 p_{n-1 = 2，这是预先计算的）。但由于这是一个惰性生成器，因此知道范围的终端边界 (pn-1+pn-2) 并不能帮助我过滤任何内容。}

在给定所有先前的素数的情况下，我可以做些什么来更有效地过滤？

代码示例

  @doc """
  Creates an infinite stream of prime numbers.

    iex> Enum.take(primes, 5)
    [2, 3, 5, 7, 11]

    iex> Enum.take_while(primes, fn(n) -> n < 25 end)
    [2, 3, 5, 7, 11, 13, 17, 19, 23]

  """
  @spec primes :: Stream.t
  def primes do
    Stream.unfold( [], fn primes ->
      next = next_prime(primes)
      { next, [next | primes] }
    end )
  end

  defp next_prime([]),      do: 2
  defp next_prime([2 | _]), do: 3
  defp next_prime([3 | _]), do: 5
  defp next_prime([5 | _]), do: 7
  # ... etc

  defp next_prime(primes) do
    start = Enum.first(primes) + 2
    Enum.first(
      Stream.drop_while(
        Integer.stream(from: start, step: 2),
        fn number ->
          Enum.any?(primes, fn prime ->
            rem(number, prime) == 0
          end )
        end
      )
    )
  end

primes 函数从一个空数组开始，获取它的下一个素数（最初是2），然后 1) 从 Stream 中发出它，2) 将它添加到用于下一个电话。 （我确信这个堆栈是一些减速的根源。）

next_primes 函数接收该堆栈。从最后一个已知的素数 +2 开始，它创建一个无限的整数流，并丢弃每个整数除以列表的任何已知素数，然后返回第一个出现的整数。

我想，这类似于惰性增量 Eratosthenes 的筛子。

您可以看到一些基本的优化尝试：我从 p_n-1+2 开始检查，然后跨过偶数。

我尝试了一个更逐字的 Eratosthenes 筛子，方法是将 Integer.stream 传递给每个计算，并在找到一个素数后，将 Integer.stream 包装在一个新的 Stream.drop_while 中，该 Stream.drop_while 只过滤掉该素数的倍数。但由于 Streams 是作为匿名函数实现的，因此破坏了调用堆栈。

值得注意的是，我并不是假设您需要所有先前的素数来生成下一个素数。多亏了我的实现，我碰巧有它们。

Answer 1

当通过 Eratosthenes 算法从素数生成复合时，您不需要所有先前的素数，只需低于当前生产点的平方根即可。

这大大降低了内存需求。那么素数就是那些不属于复合数的奇数。

每个素数 p 产生一个它的倍数链，从它的平方开始，以 2p 的步长枚举（因为我们只处理奇数）。这些倍数，每个都有其步长值，存储在字典中，从而形成优先级队列。此优先级队列中仅存在直到当前候选的平方根的素数（与 E 的分段筛的内存需求相同）。

象征性地，埃拉托色尼的筛子是

     P = {3,5,7,9, ...} \ ⋃ {{p², p²+2p, p²+4p, p²+6p, ...} | p in P}

每个奇素数通过重复加法生成其倍数的流；所有这些流合并在一起为我们提供了所有奇怪的组合；和 primes 都是没有合数的奇数（和一个偶素数，2）。

In Python（可以作为可执行伪代码读取，希望如此），

def postponed_sieve():      # postponed sieve, by Will Ness,   
    yield 2; yield 3;       #   https://stackoverflow.com/a/10733621/849891
    yield 5; yield 7;       # original code David Eppstein / Alex Martelli 
    D = {}                  #   2002, http://code.activestate.com/recipes/117119
    ps = (p for p in postponed_sieve())  # a separate Primes Supply:
    p = ps.next() and ps.next()          #   (3) a Prime to add to dict
    q = p*p                              #   (9) when its sQuare is 
    c = 9                                #   the next Candidate
    while True:
        if c not in D:                # not a multiple of any prime seen so far:
            if c < q: yield c         #   a prime, or
            else:   # (c==q):         #   the next prime's square:
                add(D,c + 2*p,2*p)    #     (9+6,6 : 15,21,27,33,...)
                p=ps.next()           #     (5)
                q=p*p                 #     (25)
        else:                         # 'c' is a composite:
            s = D.pop(c)              #   step of increment
            add(D,c + s,s)            #   next multiple, same step
        c += 2                        # next odd candidate

def add(D,x,s):                          # make no multiple keys in Dict
    while x in D: x += s                 # increment by the given step
    D[x] = s  

一旦产生了质数，它就可以被遗忘。一个单独的主要供应是从同一生成器的单独调用中获取的，递归地，以维护字典。并且该供应商的主要供应也以递归方式取自另一个。每个只需要提供到其生产点的平方根，因此总体上需要很少的生成器（大约为log log N 生成器），并且它们的大小渐近无关紧要（sqrt(N)、sqrt( sqrt(N) ) 等)。

Answer 2

对于任何数字 k，您只需要尝试用质数进行除法，包括 √k。这是因为任何大于 √k 的素数都需要乘以小于 √k 的素数 。

证明：

√k * √k = k 所以 (a+√k) * √k > k (对于所有 0)。由此可见，(a+√k) 除以 k 且仅当存在小于 √k 的除数。

这通常用于极大地加快寻找素数的速度。

Answer 3

我编写了一个程序，可以无限制地按顺序生成素数，并用它对my blog 处的前十亿个素数求和。该算法使用埃拉托色尼的分段筛；在每个段计算额外的筛选素数，因此只要您有空间存储筛选素数，该过程就可以无限期地继续。这是伪代码：

function init(delta) # Sieve of Eratosthenes
  m, ps, qs := 0, [], []
  sieve := makeArray(2 * delta, True)
  for p from 2 to delta
    if sieve[p]
      m := m + 1; ps.insert(p)
      qs.insert(p + (p-1) / 2)
      for i from p+p to n step p
        sieve[i] := False
  return m, ps, qs, sieve

function advance(m, ps, qs, sieve, bottom, delta)
  for i from 0 to delta - 1
    sieve[i] := True
  for i from 0 to m - 1
    qs[i] := (qs[i] - delta) % ps[i]
  p := ps[0] + 2
  while p * p <= bottom + 2 * delta
    if isPrime(p) # trial division
      m := m + 1; ps.insert(p)
      qs.insert((p*p - bottom - 1) / 2)
    p := p + 2
  for i from 0 to m - 1
    for j from qs[i] to delta step ps[i]
      sieve[j] := False
  return m, ps, qs, sieve

这里ps是筛选出小于当前最大值的素数列表，qs是当前段中对应ps的最小倍数的偏移量。 advance 函数清除位数组，重置 qs，扩展 ps 和 qs 使用新的筛选素数，然后筛选下一段。

function genPrimes()
  bottom, i, delta := 0, 1, 50000
  m, ps, qs, sieve := init(delta)
  yield 2
  while True
    if i == delta # reset for next segment
      i, bottom := -1, bottom + 2 * delta
      m, ps, qs, sieve := \textbackslash
        advance(m, ps, qs, sieve, bottom, delta)
    else if sieve[i] # found prime
      yield bottom + 2*i + 1
    i := i + 1

段大小2 * delta 任意设置为 100000。此方法需要 O(sqrt(n)) 空间用于筛分素数加上恒定空间用于筛子。

使用轮子生成候选对象并测试候选对象的素性比较慢，但节省空间。

function genPrimes()
  w, wheel := 0, [1,2,2,4,2,4,2,4,6,2,6,4,2,4, \
       6,6,2,6,4,2,6,4,6,8,4,2,4,2,4,8,6,4,6, \
       2,4,6,2,6,6,4,2,4,6,2,6,4,2,4,2,10,2,10]
  p := 2; yield p
  repeat
    p := p + wheel[w]
    if w == 51 then w := 4 else w := w + 1
    if isPrime(p) yield p

从筛子开始，当筛子变得太大时切换到轮子可能会很有用。更好的是继续使用一些固定的筛选素数集进行筛选，一旦该集合变得太大，则仅报告那些通过素数测试的值bottom + 2*i + 1。