python - Adding wheel factorization to an indefinite sieve -
i’m modifying indefinite sieve of eratosthenes here uses wheel factorization skip more composites current form of checking odds.
i’ve worked out how generate steps take reach gaps along wheel. there figured substitute +2’s these wheel steps it’s causing sieve skip primes. here's code:
from itertools import count, cycle def dvprm(end): "finds primes trial division. returns list" primes=[2] 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" in seq:factor*=i return factor def wheelgaps(primes): """returns list of steps each wheel gap start last value in primes""" strtpt= primes.pop(-1)#where wheel starts whlcirm= prod(primes)# wheel's circumference #spokes every number divisible primes (composites) gaps=[]#locate non-spokes (gaps) 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 steps needed jump each gap (beginning start of wheel) steps=[]#last step returns start of wheel i,j in enumerate(gaps): if i==0:continue steps.append(j - gaps[i-1]) return steps def wheel_setup(num): "builds initial data sieve" initprms=dvprm(num)#initial primes "roughing" pump gaps = wheelgaps(initprms[:])#get gaps c= initprms.pop(-1)#prime starts wheel return initprms, gaps, c def wheel_psieve(lvl=0, initdata=none): '''postponed prime generator wheels refs: http://stackoverflow.com/a/10733621 http://stackoverflow.com/a/19391111''' whlsize=11#wheel size, 1 higher prime # 5 gives 2- 3 wheel 11 gives 2- 7 wheel # 7 gives 2- 5 wheel 13 gives 2-11 wheel #set 0 no wheel if lvl:#no need rebuild gaps, pass them down levels initprms, gaps, c = initdata else:#but if top level build gaps if whlsize>4: initprms, gaps, c = wheel_setup(whlsize) else: initprms, gaps, c= dvprm(7), [2], 9 #toss out initial primes p in initprms: yield p cgaps=cycle(gaps) compost = {}#found composites skip ps=wheel_psieve(lvl+1, (initprms, gaps, c)) p=next(ps)#advance lower level appropriate square while p*p < c: p=next(ps) psq=p*p while true: step1 = next(cgaps)#step next value step2=compost.pop(c, 0)#step next multiple if not step2: #see references 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 i'm using check it:
def test(num=100): found=[] i,p in enumerate(wheel_psieve(), 1): if i>num:break found.append(p) print sum(found) return found when set wheel size 0, correct sum of 24133 first 1 hundred primes, when use other wheel size, end missing primes, incorrect sums , composites. 2-3 wheel (which uses alternate steps of 2 , 4) makes sieve miss primes. doing wrong?
the odds, i.e. 2-coprimes, generated "rolling wheel" [2], i.e. repeated additions of 2, starting initial value of 3 (similarly 5, 7, 9, ...),
n=3; n+=2; n+=2; n+=2; ... # wheel = [2] 3 5 7 9 the 2-3-coprimes generated repeated additions of 2, 4, , again 2, 4, , on:
n=5; n+=2; n+=4; n+=2; n+=4; ... # wheel = [2,4] 5 7 11 13 17 here need know start adding differences from, 2 or 4, depending on initial value. 5, 11, 17, ..., it's 2 (i.e. 0-th element of wheel); 7, 13, 19, ..., it's 4 (i.e. 1-st element).
how can know start? point wheel optimization work on sequence of coprimes (in example, 2-3-coprimes). in part of code recursively generated primes, maintain rolling wheel stream, , advance until see next prime in it. rolling sequence need produce two results - value , wheel position. when see prime, corresponding wheel position, , can start off generation of multiples starting position on wheel. multiply p of course, start p*p:
for (i, p) # (wheel position, summated value) in enumerated roll of wheel: when p next prime: multiples of p m = p*p; # map (p*) (roll wheel-at-i p) m += p*wheel[i]; m += p*wheel[i+1]; ... so each entry in dict have maintain current value, base prime, , current wheel position (wrapping around 0 circularity, when needed).
to produce resulting primes, roll coprimes sequence, , keep elements of not in dict, in reference code.
update: after few iterations on codereview (big contributors there!) i've arrived @ code, using itertools as possible, speed:
from itertools import accumulate, chain, cycle, count def wsieve(): # wheel-sieve, 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))) 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))) # start mults = {} c in cs: # candidates, coprime 210, 11 if c in mults: wheel = mults.pop(c) elif c < psq: yield c continue else: # c==psq: map (p*) (roll wh p) = roll (wh*p) (p*p) = d[(p-11) % 210] wheel = accumulate( chain([psq], cycle([p*d d in wh11[i:] + wh11[:i]]))) next(wheel) p = next(ps) psq = p**2 m in wheel: # pop, save in m, , advance if m not in mults: break mults[m] = wheel # mults[143] = wheel@187 def primes(): yield (2, 3, 5, 7) yield wsieve() unlike above description, code directly calculates start rolling wheel each prime, generate multiples
Comments
Post a Comment