patternpythonMinor
Wheel based unbounded sieve of Eratosthenes in Python
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):
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
Is this just because of the increased code complexity? Is there a way to make it run faster, without ch
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
-
Turning the wheel step by step from prime to prime seems wasteful. Instead of this code
you could do directly
-
When a wheel is involved I find it hard to resist using
you could do this
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
My version of full solution. Instead of the index
Edited to use
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 = ibaseyou 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 thisi = 0
while True:
c += wheel[i] ; i = (i+1) % wsizeyou could do this
for step in cycle(wheel):
c += stepThe 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 = ibasei = 0
while True:
c += wheel[i] ; i = (i+1) % wsizefor step in cycle(wheel):
c += stepfrom 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.