HiveBrain v1.2.0
Get Started
← Back to all entries
patternpythonMinor

Wheel based unbounded sieve of Eratosthenes in Python

Submitted by: @import:stackexchange-codereview··
0
Viewed 0 times
sieveunboundedpythonbasedwheeleratosthenes

Problem

Answering this SO question recently, I've developed the following code (based on a well-known Active Code recipe, as discussed here):

wheel = [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]
wsize = 48

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

def wsieve():       # wheel-sieve, by Will Ness. cf. ideone.com/WFv4f
    yield 11        # cf. https://stackoverflow.com/a/10733621/849891     
    mults = {}      #     https://stackoverflow.com/a/19391111/849891
    ps = wsieve() 
    p = next(ps)   
    psq, c, i = p*p, 11, 0                   # 13 = 11 + wheel[0]
    cbase, ibase = 11, 0  
    while True:
        c += wheel[i] ; i = (i+1) % wsize    # 17 = 13 + wheel[1]
        if c in mults:
            (j,pbase) = mults.pop(c)         # add(mults, NEXT c, j, pbase)
        elif c < psq:              
            yield c ; continue   
        else:          # (c==psq)       
            while not cbase == p:
                cbase += wheel[ibase]
                ibase = (ibase+1) % wsize    # ibase - initial offset into wheel, for p
            j, pbase = ibase, p              # add(mults, NEXT c, ibase, p)
            p = next(ps) ; psq = p*p
        m = c + pbase*wheel[j] ; j = (j+1) % wsize    # mults(p) = map (*p)
        while m in mults:                             #              roll(wheel,ibase,p)
            m += pbase*wheel[j] ; j = (j+1) % wsize    
        mults[m] = (j,pbase)


Comparing its performance on ideone with its odds-based predecessor reveals it runs about 1.5x faster. Empirical complexity seems okay (~ n1.1, for n primes produced).

Theoretically, it's supposed to be 1 / (2/3 4/5 6/7) = 2.1875x times faster, if I'm not mistaken (indeed, taking ratio of circumference over number of teeth, (210/48) / (2/1) = 2.1875). But it is slower than that.

Is this just because of the increased code complexity? Is there a way to make it run faster, without ch

Solution

I thought of two things that achieve a 27 % speedup and simplify wsieve a bit, at the expense of some precomputations.

-
Turning the wheel step by step from prime to prime seems wasteful. Instead of this code

while not cbase == p:
    cbase += wheel[ibase]
    ibase = (ibase+1) % wsize   
j = ibase


you could do directly j = spoke_index[p % 210] having precomputed a suitable spoke_index dictionary with 48 entries. The variables cbase and ibase can be eliminated.

-
When a wheel is involved I find it hard to resist using itertools.cycle. Instead of this

i = 0                   
while True:
    c += wheel[i] ; i = (i+1) % wsize


you could do this

for step in cycle(wheel):
    c += step


The other uses of the wheel are more difficult to replace because there is no direct way to start the cycle at an arbitrary position. However, we can precompute all 48 rotations of wheels and put them in a dictionary for easy lookup.

My version of full solution. Instead of the index j I put a cyclic iterator in mults.

from itertools import cycle

CIRCUMFERENCE = 2*3*5*7
BASE_PRIMES = (2,3,5,7)
NEXT_PRIME = 11

def wheel(start):
    result = []
    i = start
    for j in range(i + 1, i + 1 + CIRCUMFERENCE):
        if all(j % k for k in BASE_PRIMES):
            result.append(j - i)
            i = j
    return result

def rotated_wheels():
    result = {}
    i = 1
    while i < CIRCUMFERENCE:
        result[i] = wheel(i)
        i = i + result[i][0]
    return result

def primes(): 
    yield from BASE_PRIMES
    yield from wsieve()

def wsieve(wheels=rotated_wheels()):       
    yield NEXT_PRIME          
    mults = {}      
    ps = wsieve() 
    p = next(ps)   
    psq, c = p*p, p             
    cwheel = cycle(wheels[c])
    for step in cwheel:
        c += step
        if c in mults:
            (mwheel, pbase) = mults.pop(c)         
        elif c < psq:              
            yield c 
            continue   
        else:          # (c==psq)     
            mwheel = cycle(wheels[p % CIRCUMFERENCE])
            pbase = p             
            p = next(ps) ; psq = p*p
        m = c
        for mstep in mwheel:
            m += pbase * mstep
            if m not in mults:
                break
        mults[m] = (mwheel, pbase)


Edited to use for instead of while.

Code Snippets

while not cbase == p:
    cbase += wheel[ibase]
    ibase = (ibase+1) % wsize   
j = ibase
i = 0                   
while True:
    c += wheel[i] ; i = (i+1) % wsize
for step in cycle(wheel):
    c += step
from itertools import cycle

CIRCUMFERENCE = 2*3*5*7
BASE_PRIMES = (2,3,5,7)
NEXT_PRIME = 11

def wheel(start):
    result = []
    i = start
    for j in range(i + 1, i + 1 + CIRCUMFERENCE):
        if all(j % k for k in BASE_PRIMES):
            result.append(j - i)
            i = j
    return result

def rotated_wheels():
    result = {}
    i = 1
    while i < CIRCUMFERENCE:
        result[i] = wheel(i)
        i = i + result[i][0]
    return result

def primes(): 
    yield from BASE_PRIMES
    yield from wsieve()

def wsieve(wheels=rotated_wheels()):       
    yield NEXT_PRIME          
    mults = {}      
    ps = wsieve() 
    p = next(ps)   
    psq, c = p*p, p             
    cwheel = cycle(wheels[c])
    for step in cwheel:
        c += step
        if c in mults:
            (mwheel, pbase) = mults.pop(c)         
        elif c < psq:              
            yield c 
            continue   
        else:          # (c==psq)     
            mwheel = cycle(wheels[p % CIRCUMFERENCE])
            pbase = p             
            p = next(ps) ; psq = p*p
        m = c
        for mstep in mwheel:
            m += pbase * mstep
            if m not in mults:
                break
        mults[m] = (mwheel, pbase)

Context

StackExchange Code Review Q#92365, answer score: 9

Revisions (0)

No revisions yet.