The odds, i.e. 2-coprimes, are generated by “rolling the wheel” [2], i.e. by repeated additions of 2, starting from the initial value of 3 (similarly from 5, 7, 9, …),
n=3; n+=2; n+=2; n+=2; ... # wheel = [2]
3 5 7 9
The 2-3-coprimes are generated by repeated additions of 2, then 4, and again 2, then 4, and so on:
n=5; n+=2; n+=4; n+=2; n+=4; ... # wheel = [2,4]
5 7 11 13 17
Here we do need to know where to start adding the differences from, 2 or 4, depending on the initial value. For 5, 11, 17, …, it’s 2 (i.e. 0-th element of the wheel); for 7, 13, 19, …, it’s 4 (i.e. 1-st element).
How can we know where to start? The point to the wheel optimization is that we work only on this sequence of coprimes (in this example, 2-3-coprimes). So in the part of the code where we get the recursively generated primes, we will also maintain the rolling wheel stream, and advance it until we see that next prime in it. The rolling sequence will need to produce two results – the value and the wheel position. Thus when we see the prime, we also get the corresponding wheel position, and we can start off the generation of its multiples starting from that position on the wheel. We multiply everything by p of course, to start from 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]; ...
So each entry in the dict will have to maintain its current value, its base prime, and its current wheel position (wrapping around to 0 for circularity, when needed).
To produce the resulting primes, we roll another coprimes sequence, and keep only those elements of it that are not in the dict, just as in the reference code.
update: after a few iterations on codereview (big thanks to the contributors there!) I’ve arrived at this code, using itertools as much as possible, for speed:
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()
Unlike the above description, this code directly calculates where to start rolling the wheel for each prime, to generate its multiples