Добавление факторизации колес к неопределенному ситу

Im, изменяя неопределенное сито Eratoshenes от здесь, поэтому он использует факторизацию колес, чтобы пропускать больше композитов, чем его текущая форма просто проверки всех коэффициентов.

Ive разработал, как создать шаги, которые нужно предпринять, чтобы достичь всех пробелов вдоль колеса. Оттуда я решил, что могу просто заменить +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 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, я получаю правильную сумму 24133 для первых 100 простых чисел, но когда я использую любой другой размер колеса, я получаю недостающие простые числа, неправильные суммы и композиты. Даже 2-3 колесо (которое использует чередующиеся шаги 2 и 4) делает пропущенные пропуски сита. Что я делаю неправильно?

Ответы

Ответ 1

Коэффициенты, т.е. 2-совремы, генерируются путем "качения колеса" [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,..., it 2 (то есть 0-й элемент колеса); для 7, 13, 19,..., it 4 (то есть 1-й элемент).

Как мы можем знать, с чего начать? Точка к оптимизации колес заключается в том, что мы работаем только с этой последовательностью взаимных чисел (в этом примере 2-3-coprimes). Таким образом, в части кода, где мы получаем рекурсивно сгенерированные простые числа, мы также будем поддерживать поток катящегося колеса и продвигаем его, пока не увидим, что в нем будет следующее простое. Последовательность качения должна произвести два результата - значение и положение колеса. Таким образом, когда мы видим премьер, мы также получаем соответствующее положение колеса, и мы можем начать генерировать его кратность, начиная с этой позиции на колесе. Конечно, мы умножаем все на 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];    ...

Таким образом, каждая запись в dict должна поддерживать текущее значение, его базовое число и текущее положение колеса (при необходимости округляя до 0 для округлости).

Чтобы получить полученные простые числа, мы катим другую совпадающую последовательность и сохраняем только те ее элементы, которые не указаны в dict, так же, как и в ссылочном коде.

update: после нескольких итераций на codereview (большое спасибо вкладчикам!) Я пришел к этому коду, используя как можно больше 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)))
    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)))   # start from
    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]
            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] = [email protected]

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

В отличие от вышеприведенного описания, этот код непосредственно вычисляет, где начать катить колесо для каждого штриха, чтобы генерировать его кратные

Ответ 2

Это версия, с которой я столкнулся. Это не так чисто, как 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 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 wheel
        ix = cyc(ix + 1)  # and advance c 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 было около 36 миллионов номеров).