Sieve of Eratosthenes: Difference between revisions
→Fast infinite generator using a wheel: added a very large wheel with links to ideone.com
GordonBGood (talk | contribs) (→Fast infinite generator using a wheel: added links to test implementations on ideone.com...) |
GordonBGood (talk | contribs) (→Fast infinite generator using a wheel: added a very large wheel with links to ideone.com) |
||
Line 3,697:
===Fast infinite generator using a wheel===
Although theoretically over three times faster than odds-only, the following code using a 2/3/5/7 wheel is only about 1.5 times faster than the above odds-only code due to the extra overheads in code complexity
{{works with|Python|2.6+, 3.x}}
<lang python>def primes():
Line 3,729:
n += gaps[ni]; ni = nni # advance on the wheel
for p, pi in wheel_prime_pairs(): yield p # strip out indexes</lang>
Further gains of about 1.5 times in speed can be made using the same code by only changing the tables and a few constants for a further constant factor gain of about 1.5 times in speed by using a 2/3/5/7/1/13/17 wheel (with the gaps list 92160 elements long) computed for a slight constant overhead time as per the [http://ideone.com/4Ld26g test link for Python 2.7] and [http://ideone.com/72Dmyt test link for Python 3.x]. Further wheel factorization will not really be worth it as the gains will be small (if any and not losses) and the gaps table huge - it is already too big for efficient use by 32-bit Python 3 and the wheel should likely be stopped at 13:
<lang python>def primes():
whlPrms = [2,3,5,7,11,13,17] # base wheel primes
for p in whlPrms: yield p
def makeGaps():
buf = [True] * (3 * 5 * 7 * 11 * 13 * 17 + 1) # all odds plus extra for o/f
for p in whlPrms:
if p < 3:
continue # no need to handle evens
strt = (p * p - 19) >> 1 # start position (divided by 2 using shift)
while strt < 0: strt += p
buf[strt::p] = [False] * ((len(buf) - strt - 1) // p + 1) # cull for p
whlPsns = [i + i for i,v in enumerate(buf) if v]
return [whlPsns[i + 1] - whlPsns[i] for i in range(len(whlPsns) - 1)]
gaps = makeGaps() # big wheel gaps
def wheel_prime_pairs():
yield (19,0); bps = wheel_prime_pairs() # additional primes supply
p, pi = next(bps); q = p * p # adv to get 11 sqr'd is 121 as next square to put
sieve = {}; n = 23; ni = 1 # into sieve dict; init cndidate, wheel ndx
while True:
if n not in sieve: # is not a multiple of previously recorded primes
if n < q: yield (n, ni) # n is prime with wheel modulo index
else:
npi = pi + 1 # advance wheel index
if npi > 92159: npi = 0
sieve[q + p * gaps[pi]] = (p, npi) # n == p * p: put next cull position on wheel
p, pi = next(bps); q = p * p # advance next prime and prime square to put
else:
s, si = sieve.pop(n)
nxt = n + s * gaps[si] # move current cull position up the wheel
si = si + 1 # advance wheel index
if si > 92159: si = 0
while nxt in sieve: # ensure each entry is unique by wheel
nxt += s * gaps[si]
si = si + 1 # advance wheel index
if si > 92159: si = 0
sieve[nxt] = (s, si) # next non-marked multiple of a prime
nni = ni + 1 # advance wheel index
if nni > 92159: nni = 0
n += gaps[ni]; ni = nni # advance on the wheel
for p, pi in wheel_prime_pairs(): yield p # strip out indexes
</lang>
=={{header|R}}==
| |||