如何在Python中实现有效的素数无限生成器？

Question

这不是作业,我只是好奇.

INFINITE是这里的关键词.

我希望在primes()中使用它作为p.我相信这是Haskell中的内置函数.

所以,答案不能像"Just do a Sieve"那样天真.

首先,您不知道将消耗多少连续素数.好吧,假设你可以一次编制100个.您是否会使用相同的Sieve方法以及素数公式的频率？

我更喜欢非并发方法.

感谢您阅读(和写作;))!

Answer 1

"如果我再看到......"

erat2食谱中的功能可以进一步加快(约20-25%):

erat2a

import itertools as it
def erat2a( ):
    D = {  }
    yield 2
    for q in it.islice(it.count(3), 0, None, 2):
        p = D.pop(q, None)
        if p is None:
            D[q*q] = q
            yield q
        else:
            # old code here:
            # x = p + q
            # while x in D or not (x&1):
            #     x += p
            # changed into:
            x = q + 2*p
            while x in D:
                x += 2*p
            D[x] = p

该not (x&1)检查验证x为奇数.然而,由于这两个 q和p是奇数,通过添加2*p下列步骤一半与测试古怪沿着避免.

erat3

如果一个人不介意一点额外的幻想,erat2可以通过以下更改加速35-40%(注意:由于该itertools.compress功能,需要Python 2.7+或Python 3+ ):

import itertools as it
def erat3( ):
    D = { 9: 3, 25: 5 }
    yield 2
    yield 3
    yield 5
    MASK= 1, 0, 1, 1, 0, 1, 1, 0, 1, 0, 0, 1, 1, 0, 0,
    MODULOS= frozenset( (1, 7, 11, 13, 17, 19, 23, 29) )

    for q in it.compress(
            it.islice(it.count(7), 0, None, 2),
            it.cycle(MASK)):
        p = D.pop(q, None)
        if p is None:
            D[q*q] = q
            yield q
        else:
            x = q + 2*p
            while x in D or (x%30) not in MODULOS:
                x += 2*p
            D[x] = p

该erat3函数利用了这样的事实:所有素数(除了MODULOS2,3,5 除外)模数为30只导致8个数:冻结集中包含的数.因此,在产生最初的三个素数后,我们从7开始,只与候选人一起工作.
候选过滤使用该itertools.compress功能; "神奇"在MASK序列中; MASK有15个元素(每30个数字中有15个奇数,由itertools.islice函数选择),1每个可能的候选者都有一个元素,从7开始.循环按itertools.cycle函数指定重复.
候选过滤的引入需要另一种修改:or (x%30) not in MODULOS检查.该erat2算法处理所有奇数; 既然erat3算法只处理r30候选者,我们需要确保所有人D.keys()只能是这样的假候选者.

基准

结果

在Atom 330 Ubuntu 9.10服务器上,版本2.6.4和3.1.1+:

$ testit
up to 8192
==== python2 erat2 ====
100 loops, best of 3: 18.6 msec per loop
==== python2 erat2a ====
100 loops, best of 3: 14.5 msec per loop
==== python2 erat3 ====
Traceback (most recent call last):
…
AttributeError: 'module' object has no attribute 'compress'
==== python3 erat2 ====
100 loops, best of 3: 19.2 msec per loop
==== python3 erat2a ====
100 loops, best of 3: 14.1 msec per loop
==== python3 erat3 ====
100 loops, best of 3: 11.7 msec per loop

在AMD Geode LX Gentoo家庭服务器上,Python 2.6.5和3.1.2:

$ testit
up to 8192
==== python2 erat2 ====
10 loops, best of 3: 104 msec per loop
==== python2 erat2a ====
10 loops, best of 3: 81 msec per loop
==== python2 erat3 ====
Traceback (most recent call last):
…
AttributeError: 'module' object has no attribute 'compress'
==== python3 erat2 ====
10 loops, best of 3: 116 msec per loop
==== python3 erat2a ====
10 loops, best of 3: 82 msec per loop
==== python3 erat3 ====
10 loops, best of 3: 66 msec per loop

基准代码

甲primegen.py模块包含erat2,erat2a和erat3功能.以下是测试脚本:

#!/bin/sh
max_num=${1:-8192}
echo up to $max_num
for python_version in python2 python3
do
    for function in erat2 erat2a erat3
    do
        echo "==== $python_version $function ===="
        $python_version -O -m timeit -c \
        -s  "import itertools as it, functools as ft, operator as op, primegen; cmp= ft.partial(op.ge, $max_num)" \
            "next(it.dropwhile(cmp, primegen.$function()))"
    done
done

Answer 2

由于OP要求有效实现,所以这里是David Eppstein/Alex Martelli(在他的回答中看到)对活动状态2002代码的重大改进:不要在字典中记录素数的信息,直到它的正方形被看到候选人.对于产生的n个素数(π(sqrt(n log n)) ~ 2 sqrt(n log n)/ log(n log n)〜,将空间复杂度降低到O(sqrt(n))以下而不是O(n).2 sqrt(n/log n)).因此,时间复杂性也得到改善,即运行得更快.

创建一个"滑动筛"作为每个基本素数的当前倍数的字典(即低于当前生产点的sqrt),以及它们的步长值:

from itertools import count
                                         # ideone.com/aVndFM
def postponed_sieve():                   # postponed sieve, by Will Ness      
    yield 2; yield 3; yield 5; yield 7;  # original code David Eppstein, 
    sieve = {}                           #   Alex Martelli, ActiveState Recipe 2002
    ps = postponed_sieve()               # a separate base Primes Supply:
    p = next(ps) and next(ps)            # (3) a Prime to add to dict
    q = p*p                              # (9) its sQuare 
    for c in count(9,2):                 # the Candidate
        if c in sieve:               # c's a multiple of some base prime
            s = sieve.pop(c)         #     i.e. a composite ; or
        elif c < q:  
             yield c                 # a prime
             continue              
        else:   # (c==q):            # or the next base prime's square:
            s=count(q+2*p,2*p)       #    (9+6, by 6 : 15,21,27,33,...)
            p=next(ps)               #    (5)
            q=p*p                    #    (25)
        for m in s:                  # the next multiple 
            if m not in sieve:       # no duplicates
                break
        sieve[m] = s                 # original test entry: ideone.com/WFv4f

(年纪大了,原来这里的代码编辑为包括变化中看到了答案由蒂姆·彼得斯,下同).另见这个相关的讨论.

类似的2-3-5-7基于轮的代码运行速度快〜2.15倍(这非常接近理论上的改进3/2 * 5/4 * 7/6 = 2.1875).

Answer 3

对于后人来说,这里是Will Ness为Python 3编写的漂亮算法的重写.需要进行一些更改(迭代器不再有.next()方法,但有一个新的next()内置函数).其他变化是为了好玩(使用new yield from <iterable>替换yield原始中的四个语句.更多是为了可读性(我不是过度使用的粉丝;-) 1个字母的变量名称).

它明显快于原始速度,但不是出于算法原因.加速主要是因为删除原始的add()功能,而是内联.

def psieve():
    import itertools
    yield from (2, 3, 5, 7)
    D = {}
    ps = psieve()
    next(ps)
    p = next(ps)
    assert p == 3
    psq = p*p
    for i in itertools.count(9, 2):
        if i in D:      # composite
            step = D.pop(i)
        elif i < psq:   # prime
            yield i
            continue
        else:           # composite, = p*p
            assert i == psq
            step = 2*p
            p = next(ps)
            psq = p*p
        i += step
        while i in D:
            i += step
        D[i] = step

Answer 4

这不是我原来的代码,但值得张贴.原文可以在这里找到:http://code.activestate.com/recipes/117119/

def gen_primes():
  D = {}
  q = 2  # first integer to test for primality.

  while True:
    if q not in D:
      # not marked composite, must be prime  
      yield q 

      #first multiple of q not already marked
      D[q * q] = [q] 
    else:
      for p in D[q]:
        D.setdefault(p + q, []).append(p)
      # no longer need D[q], free memory
      del D[q]

    q += 1

它是一个发电机,所以像其他任何一样使用它.

primes = gen_primes()
for p in primes:
  print p

在我的桌面上生成并放入一组100万个素数需要1.62秒.

Answer 5

做一个分段筛,其中段的大小由可用内存或bitset的最大大小决定.

对于每个段,表示某个区间中的数字[n; n + segment_size)作为位集和筛子,所有素数低于上限的平方根.

使用位集比使用散列表或树数据结构使用更少的内存,因为您正在处理密集的数字集.