当前位置: 首页 > 知识库问答 >
问题:

将砂轮因式分解添加到不确定筛中

蒋栋
2023-03-14

我正在从这里修改埃拉托斯特尼的无限筛子,因此它使用车轮分解来跳过比当前仅检查所有赔率的形式更多的复合材料。

我已经想出了如何生成到达轮子上所有间隙的步骤。从那里,我想我可以用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个答案

贲培
2023-03-14

这是我想到的版本。它不像Ness那样干净,但是很有效。我把它贴出来,这样就有了另一个如何使用轮子因式分解的例子,以防有人路过。我可以选择使用什么样的轮子尺寸,但是很容易确定一个更永久的尺寸——只需生成您想要的尺寸并粘贴到代码中。

from itertools import count

def wpsieve():
    """prime number generator
    call this function instead of roughing or turbo"""
    whlSize = 11
    initPrms, gaps, c = wheel_setup(whlSize)

    for p in initPrms:
        yield p

    primes = turbo(0, (gaps, c))

    for p, x in primes:
        yield p

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(int(j - gaps[i-1]))
    return steps

def wheel_setup(num):
    "builds initial data for sieve"
    initPrms = roughing(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 roughing(end):
    "finds primes by trial division (roughing pump)"
    primes = [2]
    for i in range(3, end + 1, 2):
        if all(map(lambda x: i%x, primes)):
            primes.append(i)
    return primes

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

    gaps, c = initData

    yield (c, 0)

    compost = {}  # found composites to skip
    # store as current value: (base prime, wheel index)

    ps = turbo(lvl + 1, (gaps, c))

    p, x = next(ps)
    psq = p*p
    gapS = len(gaps) - 1

    ix = jx = kx = 0  # indices for cycling the wheel

    def cyc(x): return 0 if x > gapS else x  # wheel cycler

    while True:
        c += gaps[ix]  # add next step on c's wheel
        ix = cyc(ix + 1)  # and advance c's index

        bp, jx = compost.pop(c, (0,0))  # get base prime and its wheel index

        if not bp:

            if c < psq:  # prime
                yield c, ix  # emit index for above recursive level
                continue
            else:
                jx = kx  # swap indices as a new prime comes up
                bp = p
                p, kx = next(ps)
                psq = p*p

        d = c + bp * gaps[jx]  # calc new multiple
        jx = cyc(jx + 1)

        while d in compost:
            step = bp * gaps[jx]
            jx = cyc(jx + 1)
            d += step

        compost[d] = (bp, jx)

保留车轮尺寸的选项还可以让您看到较大的车轮不会起多大作用。下面是测试代码,说明生成选定尺寸的砂轮需要多长时间,以及筛子与该轮子的速度有多快。

import time
def speed_test(num, whlSize):

    print('-'*50)

    t1 = time.time()
    initPrms, gaps, c = wheel_setup(whlSize)
    t2 = time.time()

    print('2-{} wheel'.format(initPrms[-1]))
    print('setup time: {} sec.'.format(round(t2 - t1, 5)))

    t3 = time.time()      
    prm = initPrms[:]
    primes = turbo(0, (gaps, c))
    for p, x in primes:
        prm.append(p)
        if len(prm) > num:
            break
    t4 = time.time()

    print('run time  : {} sec.'.format(len(prm), round(t4 - t3, 5)))
    print('prime sum : {}'.format(sum(prm)))

for w in [5, 7, 11, 13, 17, 19, 23, 29]:
    speed_test(1e7-1, w)

当设置为生成一千万个素数时,它是如何使用PyPy(兼容Python 2.7)在我的计算机上运行的:

2- 3 wheel
setup time:  0.0 sec.
run time  : 18.349 sec.
prime sum : 870530414842019
--------------------------------------------------
2- 5 wheel
setup time:  0.001 sec.
run time  : 13.993 sec.
prime sum : 870530414842019
--------------------------------------------------
2- 7 wheel
setup time:  0.001 sec.
run time  :  7.821 sec.
prime sum : 870530414842019
--------------------------------------------------
2- 11 wheel
setup time:  0.03 sec.
run time  :  6.224 sec.
prime sum : 870530414842019
--------------------------------------------------
2- 13 wheel
setup time:  0.011 sec.
run time  :  5.624 sec.
prime sum : 870530414842019
--------------------------------------------------
2- 17 wheel
setup time:  0.047 sec.
run time  :  5.262 sec.
prime sum : 870530414842019
--------------------------------------------------
2- 19 wheel
setup time:  1.043 sec.
run time  :  5.119 sec.
prime sum : 870530414842019
--------------------------------------------------
2- 23 wheel
setup time: 22.685 sec.
run time  :  4.634 sec.
prime sum : 870530414842019

更大的轮子是可能的,但你可以看到他们变得相当长的设置。随着轮子变大,回报率也在递减——没有太多的点可以超过2-13轮子,因为它们并没有使它变得更快。我还遇到了一个2-23轮的内存错误(它的gaps列表中有大约3600万个数字)。

孟英光
2023-03-14

赔率,即2-coprimes,是通过“滚动轮子”[2]产生的,即通过重复添加2,从初始值3开始(类似地从5,7,9,…),

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-st 元素)。

我们如何知道从哪里开始?轮子优化的要点是我们只处理这个共素序列(在这个例子中,2-3个共素)。因此,在我们获得递归生成的素数的代码部分,我们还将维护滚动的车轮流,并推进它,直到我们看到其中的下一个素数。滚动序列将需要产生两个结果-值和车轮位置。因此,当我们看到质数时,我们也得到相应的轮位置,并且我们可以从轮上的该位置开始产生它的倍数。我们用< code>p乘以所有的东西,当然,从< code>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];    ...

因此,dict中的每个条目都必须保持其当前值、基本素数和当前轮位置(如果需要,将其环绕为0表示圆度)。

为了生成得到的素数,我们滚动另一个共素数序列,并仅保留其中不在字典中的那些元素,就像在参考代码中一样。

更新:在codereview上进行了几次迭代之后(非常感谢那里的贡献者!为了速度,我尽可能多地使用迭代工具,得出了这段代码:

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

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

 类似资料:
  • 问题内容: 我正在从此处修改一个不确定的Eratosthenes筛子,因此与目前仅检查所有赔率的形式相比,它使用车轮分解来跳过更多合成。 我已经弄清楚了如何产生步骤以达到车轮上所有的间隙。从那里我可以算出我可以用+2代替这些滚轮,但这会导致筛子漏掉素数。这是代码: 我正在用它来检查它: 当我将车轮尺寸设置为0时,对于前一百个质数,我得到正确的总和24133,但是当我使用任何其他车轮尺寸时,我最终会

  • 我在PrimeFaces3.5中使用JSF。我使用而不使用属性,而是显式地使用和创建行和列,如展示案例(http://www.primefaces.org/showcase/ui/panelgrid.jsf)所示。 现在我需要在一些CSS类的帮助下对一行进行不同的样式设置。但是要么我错过了它,要么就是没有办法向添加类?!我甚至可以设置属性,但在呈现的输出中会忽略它...

  • 我正在实现一个相当快的质数生成器,我得到了一些不错的结果,在埃拉托斯特尼的筛子上进行了一些优化。特别是,在算法的初步部分,我以这种方式跳过2和3的所有倍数: 这里是一个根据埃拉托色尼筛的布尔数组。我认为这是一种只考虑质数2和3的轮式因式分解,按照模式2、4、2、4递增。. 我想做的是实现一个更大的轮子,也许考虑素数2,3和5。 我已经阅读了很多关于它的文档,但我没有看到任何使用埃拉托斯特尼筛子的实

  • 我已经搜索了又搜索,除了我称之为“hack方法”的方法之外,找不到其他方法将自定义分类添加到自定义管理菜单中。 然后我注册我的帖子类型并确保它们使用 这可以工作,自定义帖子类型显示在我的自定义菜单中。 但是自定义分类法不接受同一属性的字符串,只接受true或false。 因此,要添加它,您必须创建一个子菜单页 这是一种“黑客”方式。 还有别的办法吗?如果不修改WordPress核心,我可以覆盖re

  • 问题内容: 我需要能够向ELF文件添加任意部分。我无法在该程序中使用GPL代码,因此BFD成为不可能。我可以使用libelf / gelf来阅读节,但是这些文档很少,因此我无法弄清楚如何添加节。有人知道怎么做这个吗?我宁愿不编写自己的ELF代码。 问题答案: 关于ELF文件标头的问题有几个(可能)相关答案。提到的接受的答案用于将部分添加到ELF文件,并且BSD bintools 声称具有BSD许可

  • 我想在gradle.properties文件中添加ndk.abifilters属性。现在我有了这栋楼里的房子。这是我建筑的一部分 这是我的gradle.properties文件的一部分