python - Agregar factorización de rueda a un tamiz indefinido

primes infinite

Question

Estoy modificando un tamiz indefinido de Eratóstenes desde here por lo que utiliza la factorización de la rueda para omitir más compuestos que su forma actual de simplemente verificar todas las probabilidades.

He resuelto cómo generar los pasos a seguir para alcanzar todos los huecos a lo largo de la rueda. A partir de ahí, pensé que podría sustituir los + 2 por estos pasos de rueda, pero está causando que el tamiz se salte los cebos. Aquí está el código:

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

Estoy usando esto para verificarlo:

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

Cuando configuro el tamaño de la rueda en 0, obtengo la suma correcta de 24133 para los primeros cien números primos, pero cuando uso cualquier otro tamaño de rueda, termino con números primos faltantes, sumas incorrectas y compuestos. Incluso una rueda 2-3 (que utiliza pasos alternativos de 2 y 4) hace que el tamiz pierda el cebado. ¿Qué estoy haciendo mal?

Answer 1

Esta es la versión que se me ocurrió. No es tan limpio como el de Ness , pero funciona. Lo estoy publicando para que haya otro ejemplo sobre cómo usar la factorización de la rueda en caso de que alguien venga. Me queda la posibilidad de elegir qué tamaño de rueda usar, pero es fácil establecer una más permanente: solo genera el tamaño que deseas y pégalo en el código.

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)

dejar la opción para el tamaño de la rueda también le permite ver qué tan rápido las ruedas más grandes no hacen mucho. A continuación se muestra el código de prueba de cuánto tiempo se tarda en generar la rueda del tamaño seleccionado y qué tan rápido es el tamiz con esa rueda.

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)

Así es como se ejecutó en mi computadora usando PyPy (compatible con Python 2.7) cuando está configurado para generar diez millones de primos:

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

Son posibles ruedas más grandes, pero puede ver que se vuelven bastante largas de instalar. También existe la ley de rendimientos decrecientes a medida que las ruedas se hacen más grandes, no tiene mucho sentido pasar la rueda 2-13 ya que en realidad no lo hacen mucho más rápido. También terminé encontrando un error de memoria más allá de la rueda 2-23 (que tenía unos 36 millones de números en su lista de gaps ).

Answer 2

Las probabilidades, es decir, 2 coprimes, se generan al "hacer girar la rueda" [2] , es decir, mediante adiciones repetidas de 2, comenzando desde el valor inicial de 3 (de manera similar a partir de 5, 7, 9, ...),

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

Los 2-3 coprimos se generan mediante adiciones repetidas de 2, luego 4, y nuevamente 2, luego 4, y así sucesivamente:

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

Aquí necesitamos saber por dónde comenzar a sumar las diferencias, 2 o 4, dependiendo del valor inicial. Para 5, 11, 17, ..., es 2 (es decir, el elemento número 0 de la rueda); para 7, 13, 19, ..., es 4 (es decir, primer elemento).

¿Cómo podemos saber por dónde empezar? El punto para la optimización de la rueda es que solo trabajamos en esta secuencia de coprimos (en este ejemplo, 2-3 coprimos). Entonces, en la parte del código donde obtenemos los primos generados recursivamente, también mantendremos el flujo de la rueda rodante y avanzaremos hasta que veamos el próximo cebado en él. La secuencia de balanceo necesitará producir dos resultados: el valor y la posición de la rueda. Por lo tanto, cuando vemos la prima, también obtenemos la posición correspondiente de la rueda, y podemos comenzar la generación de sus múltiplos a partir de esa posición en la rueda. Multiplicamos todo por p por supuesto, para comenzar desde 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];    ...

Por lo tanto, cada entrada en el dict tendrá que mantener su valor actual, su prima base y su posición actual de la rueda (ajustando a 0 para circularidad, cuando sea necesario).

Para producir los primos resultantes, rodamos otra secuencia coprimes, y conservamos solo aquellos elementos que no están en el dict, como en el código de referencia.

actualización: después de algunas iteraciones en codereview (¡ muchas gracias a los colaboradores allí!) He llegado a este código, usando itertools tanto como sea posible, para la velocidad:

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

A diferencia de la descripción anterior, este código calcula directamente dónde comenzar a girar la rueda para cada cebado, para generar sus múltiplos