将轮分解添加到不定筛

我正在修改埃拉托色尼的不定筛here https://stackoverflow.com/a/10733621因此，它使用轮分解来跳过比当前仅检查所有赔率的形式更多的组合。

我已经弄清楚了如何生成到达轮子上所有间隙所需的步骤。从那里我想我可以用 +2 代替这些轮子步骤，但它会导致筛子跳过素数。这是代码：

from itertools import count, cycle

def dvprm(end):
    "finds primes by trial division. returns a list"
    primes=[2]
    for i in range(3, end+1, 2):
        if all(map(lambda x:i%x, primes)):
            primes.append(i)
    return primes

def prod(seq, factor=1):
    "sequence -> product"
    for i in seq:factor*=i
    return factor

def wheelGaps(primes):
    """returns list of steps to each wheel gap
    that start from the last value in primes"""
    strtPt= primes.pop(-1)#where the wheel starts
    whlCirm= prod(primes)# wheel's circumference

    #spokes are every number that are divisible by primes (composites)
    gaps=[]#locate where the non-spokes are (gaps)
    for i in xrange(strtPt, strtPt+whlCirm+1, 2):
        if not all(map(lambda x:i%x,primes)):continue#spoke 
        else: gaps.append(i)#non-spoke

    #find the steps needed to jump to each gap (beginning from the start of the wheel)
    steps=[]#last step returns to start of wheel
    for i,j in enumerate(gaps):
        if i==0:continue
        steps.append(j - gaps[i-1])
    return steps

def wheel_setup(num):
    "builds initial data for sieve"
    initPrms=dvprm(num)#initial primes from the "roughing" pump
    gaps = wheelGaps(initPrms[:])#get the gaps
    c= initPrms.pop(-1)#prime that starts the wheel

    return initPrms, gaps, c

def wheel_psieve(lvl=0, initData=None):
    '''postponed prime generator with wheels
    Refs:  http://stackoverflow.com/a/10733621
           http://stackoverflow.com/a/19391111'''

    whlSize=11#wheel size, 1 higher prime than
    # 5 gives 2- 3 wheel      11 gives 2- 7 wheel 
    # 7 gives 2- 5 wheel      13 gives 2-11 wheel
    #set to 0 for no wheel

    if lvl:#no need to rebuild the gaps, just pass them down the levels
        initPrms, gaps, c = initData
    else:#but if its the top level then build the gaps
        if whlSize>4:
            initPrms, gaps, c = wheel_setup(whlSize) 
        else:
            initPrms, gaps, c= dvprm(7), [2], 9

    #toss out the initial primes
    for p in initPrms:
        yield p

    cgaps=cycle(gaps)
    compost = {}#found composites to skip

    ps=wheel_psieve(lvl+1, (initPrms, gaps, c))

    p=next(ps)#advance lower level to appropriate square
    while p*p < c:
        p=next(ps)
    psq=p*p

    while True:
        step1 = next(cgaps)#step to next value

        step2=compost.pop(c, 0)#step to next multiple

        if not step2:

            #see references for details
            if c < psq:
                yield c
                c += step1
                continue

            else:
                step2=2*p
                p=next(ps)
                psq=p*p

        d = c + step2
        while d in compost:
            d+= step2
        compost[d]= step2

        c += step1

我用这个来检查它：

def test(num=100):
    found=[]
    for i,p in enumerate(wheel_psieve(), 1):
        if i>num:break
        found.append(p)

    print sum(found)
    return found

当我将轮大小设置为 0 时，我得到前 100 个素数的正确总和 24133，但是当我使用任何其他轮大小时，我最终会得到素数丢失、总和和合成错误。即使是 2-3 轮（使用 2 和 4 的交替步骤）也会使筛子错过素数。我究竟做错了什么？

赔率，即 2-互质数，由下式生成“滚轮子” [2]，即从初始值 3 开始（类似地从 5、7、9、...）开始重复添加 2，

n=3; n+=2; n+=2; n+=2; ...           # wheel = [2]
  3     5     7     9

2-3-互质是通过重复添加 2、然后 4、再添加 2、然后 4 等生成的：

n=5; n+=2; n+=4; n+=2; n+=4; ...     # wheel = [2,4]
  5     7    11    13    17

这里我们确实需要知道从哪里开始添加 2 或 4 的差异，具体取决于初始值。对于 5、11、17...，它是 2（即轮盘的第 0 个元素）；对于 7、13、19...，它是 4（即第 1 个元素）。

我们怎样才能知道从哪里开始呢？轮优化的要点是我们只处理这个互质序列（在本例中为 2-3-互质）。因此，在获得递归生成素数的代码部分中，我们还将维护滚轮流，并推进它，直到看到其中的下一个素数。轧制顺序需要产生two结果 - 值和车轮位置。因此，当我们看到质数时，我们也得到了相应的轮子位置，并且我们可以从轮子上的该位置开始生成其倍数。我们将所有内容乘以p当然，要从p*p:

for (i, p) # the (wheel position, summated value) 
           in enumerated roll of the wheel:
  when p is the next prime:
    multiples of p are m =  p*p;       # map (p*) (roll wheel-at-i from p)
                       m += p*wheel[i]; 
                       m += p*wheel[i+1];    ...

因此，字典中的每个条目都必须保持其当前值、其基素数和当前轮位置（在需要时，环绕到 0 以实现循环）。

为了产生结果素数，我们滚动另一个互素序列，并只保留其中不在字典中的元素，就像参考代码中一样。

update:经过几次迭代进行代码审查 https://codereview.stackexchange.com/q/92365/9064（非常感谢那里的贡献者！）为了提高速度，我尽可能多地使用 itertools 得到了这段代码：

from itertools import accumulate, chain, cycle, count
def wsieve():  # wheel-sieve, by Will Ness.    ideone.com/mqO25A

    wh11 = [ 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]
    cs = accumulate(chain([11], cycle(wh11)))    # roll the wheel from 11
    yield(next(cs))       # cf. ideone.com/WFv4f,
    ps = wsieve()         # codereview.stackexchange.com/q/92365/9064
    p = next(ps)          # 11
    psq = p**2            # 121
    D = dict(zip(accumulate(chain([0], wh11)), count(0)))  # wheel roll lookup dict
    mults = {}
    for c in cs:          # candidates, coprime with 210, from 11
        if c in mults:
            wheel = mults.pop(c)
        elif c < psq:
            yield c
            continue
        else:    # c==psq:  map (p*) (roll wh from p) = roll (wh*p) from (p*p)
            i = D[(p-11) % 210]                 # look up wheel roll starting point
            wheel = accumulate( chain( [psq], 
                             cycle( [p*d for d in wh11[i:] + wh11[:i]])))
            next(wheel)
            p = next(ps)
            psq = p**2
        for m in wheel:   # pop, save in m, and advance
            if m not in mults:
                break
        mults[m] = wheel  # mults[143] = wheel@187

def primes():
    yield from (2, 3, 5, 7)
    yield from wsieve()

与上面的描述不同，这段代码直接计算每个素数从哪里开始滚动，以生成它的倍数