在不确定的筛子上增加车轮分解 [英] Adding wheel factorization to an indefinite sieve

问题描述

我正在此处修改一个不确定的Eratosthenes筛子，因此它使用车轮分解来跳过比当前形式更多的合成只是检查所有赔率.

I’m modifying an indefinite sieve of Eratosthenes from here so it uses wheel factorization to skip more composites than its current form of just checking all odds.

我已经弄清楚了如何生成要沿着方向盘上的所有间隙采取的步骤.从那里，我认为我可以用+2代替这些滚轮，但这会导致筛子漏掉素数.这是代码:

I’ve worked out how to generate the steps to take to reach all the gaps along the wheel. From there I figured I could just substitute the +2’s for these wheel steps but it’s causing the sieve to skip primes. Here's the code:

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

我正在使用它进行检查:

I'm using this to check it:

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时，对于前一百个质数，我得到正确的总和24133，但是当我使用任何其他车轮尺寸时，我最终会丢失质数，错误的总和和复合物.即使是2-3轮(使用2和4的交替步骤)也会使筛子漏注.我在做什么错了?

When I set the wheel size to 0, I get the correct sum of 24133 for the first one hundred primes, but when I use any other wheel size, I end up with missing primes, incorrect sums and composites. Even a 2-3 wheel (which uses alternate steps of 2 and 4) makes the sieve miss primes. What am I doing wrong?

推荐答案

赔率(即2互质数)是通过滚轮" [2]产生的，即重复加2 ，从初始值3(类似于5、7、9，...)开始，

The odds, i.e. 2-coprimes, are generated by "rolling the wheel" [2], i.e. by repeated additions of 2, starting from the initial value of 3 (similarly from 5, 7, 9, ...),

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

通过重复相加2，然后4，再重复2，然后4等生成2-3-coprimes:

The 2-3-coprimes are generated by repeated additions of 2, then 4, and again 2, then 4, and so on:

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(即第一个元素).

Here we do need to know where to start adding the differences from, 2 or 4, depending on the initial value. For 5, 11, 17, ..., it's 2 (i.e. 0-th element of the wheel); for 7, 13, 19, ..., it's 4 (i.e. 1-st element).

我们怎么知道从哪里开始?滚轮优化的意义在于，我们仅处理该序列互质数(在本示例中为2-3-coprimes).因此，在代码的一部分中，我们获得了递归生成的质数，我们还将维护滚轮流，并对其进行推进，直到看到其中的下一个质数为止.滚动序列将需要产生两个结果-值和车轮位置.因此，当我们看到质数时，我们还获得了相应的车轮位置，并且可以从车轮上的该位置开始生成其倍数.我们当然将所有内容乘以p，从p*p开始:

How can we know where to start? The point to the wheel optimization is that we work only on this sequence of coprimes (in this example, 2-3-coprimes). So in the part of the code where we get the recursively generated primes, we will also maintain the rolling wheel stream, and advance it until we see that next prime in it. The rolling sequence will need to produce two results - the value and the wheel position. Thus when we see the prime, we also get the corresponding wheel position, and we can start off the generation of its multiples starting from that position on the wheel. We multiply everything by p of course, to start from 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为圆度).

So each entry in the dict will have to maintain its current value, its base prime, and its current wheel position (wrapping around to 0 for circularity, when needed).

要产生最终的质数，我们滚动另一个互质数序列，并仅保留字典中未包含的那些元素，就像参考代码中那样.

To produce the resulting primes, we roll another coprimes sequence, and keep only those elements of it that are not in the dict, just as in the reference code.

更新:经过几次迭代关于codereview (感谢贡献者那里！)为了速度，我尽可能地使用itertools到达了这段代码:

update: after a few iterations on codereview (big thanks to the contributors there!) I've arrived at this code, using itertools as much as possible, for speed:

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()

与上面的描述不同，此代码直接计算每个素数从何处开始滚动轮子，以生成其倍数

Unlike the above description, this code directly calculates where to start rolling the wheel for each prime, to generate its multiples

这篇关于在不确定的筛子上增加车轮分解的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！