我正在从此处修改一个不确定的Eratosthenes筛子,因此与目前仅检查所有赔率的形式相比,它使用车轮分解来跳过更多合成。
我已经弄清楚了如何产生步骤以达到车轮上所有的间隙。从那里我可以算出我可以用+2代替这些滚轮,但这会导致筛子漏掉素数。这是代码:
from itertools import count, cycle def dvprm(end): "finds primes by trial division. returns a list" primes=[2] for i 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" for i in seq:factor*=i return factor def wheelGaps(primes): """returns list of steps to each wheel gap that start from the last value in primes""" strtPt= primes.pop(-1)#where the wheel starts whlCirm= prod(primes)# wheel's circumference #spokes are every number that are divisible by primes (composites) gaps=[]#locate where the non-spokes are (gaps) for i 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 the steps needed to jump to each gap (beginning from the start of the wheel) steps=[]#last step returns to start of wheel for i,j in enumerate(gaps): if i==0:continue steps.append(j - gaps[i-1]) return steps def wheel_setup(num): "builds initial data for sieve" initPrms=dvprm(num)#initial primes from the "roughing" pump gaps = wheelGaps(initPrms[:])#get the gaps c= initPrms.pop(-1)#prime that starts the wheel return initPrms, gaps, c def wheel_psieve(lvl=0, initData=None): '''postponed prime generator with wheels Refs: http://stackoverflow.com/a/10733621 http://stackoverflow.com/a/19391111''' whlSize=11#wheel size, 1 higher prime than # 5 gives 2- 3 wheel 11 gives 2- 7 wheel # 7 gives 2- 5 wheel 13 gives 2-11 wheel #set to 0 for no wheel if lvl:#no need to rebuild the gaps, just pass them down the levels initPrms, gaps, c = initData else:#but if its the top level then build the gaps if whlSize>4: initPrms, gaps, c = wheel_setup(whlSize) else: initPrms, gaps, c= dvprm(7), [2], 9 #toss out the initial primes for p in initPrms: yield p cgaps=cycle(gaps) compost = {}#found composites to skip ps=wheel_psieve(lvl+1, (initPrms, gaps, c)) p=next(ps)#advance lower level to appropriate square while p*p < c: p=next(ps) psq=p*p while True: step1 = next(cgaps)#step to next value step2=compost.pop(c, 0)#step to next multiple if not step2: #see references for 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
我正在用它来检查它:
def test(num=100): found=[] for i,p in enumerate(wheel_psieve(), 1): if i>num:break found.append(p) print sum(found) return found
当我将车轮尺寸设置为0时,对于前一百个质数,我得到正确的总和24133,但是当我使用任何其他车轮尺寸时,我最终会丢失质数,错误的总和和合成。即使是2-3轮(使用2和4的交替步骤)也会使筛子漏注。我究竟做错了什么?
赔率(即2互质数)是通过 “滚轮” 产生的[2],即通过从初始值3(类似于5、7、9 …)开始重复添加2,
[2]
n=3; n+=2; n+=2; n+=2; ... # wheel = [2] 3 5 7 9
通过重复相加2,然后4,再重复2,然后4等生成2-3-coprimes:
n=5; n+=2; n+=4; n+=2; n+=4; ... # wheel = [2,4] 5 7 11 13 17
在这里,我们确实需要知道从哪里开始添加与2或4之间的差异,具体取决于初始值。对于5、11、17 …,它是2(即车轮的第0个元素);对于7、13、19,…,它是4(即第一元素)。
我们怎么知道从哪里开始?滚轮优化的意义在于,我们仅处理此序列互质数(在此示例中为2-3-coprimes)。因此,在获得递归生成的质数的代码部分中,我们还将维护滚轮流,并对其进行推进,直到看到其中的下一个质数为止。滚动顺序将需要产生 两个 结果- 数值和车轮位置。因此,当我们看到质数时,我们也获得了相应的车轮位置,并且可以从车轮上的该位置开始生成其倍数。p当然,我们将所有内容相乘,从开始p*p:
p
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]; ...
因此,字典中的每个条目都必须保持其当前值,其基本质数和其当前的车轮位置(在需要时,将其环绕度调整为0)。
为了产生结果质数,我们滚动另一个互质数序列,并仅保留字典中未包含的那些元素,就像参考代码中那样。
更新: 在codereview上进行了几次迭代(非常感谢那里的贡献者!)我已经到达此代码,并尽可能使用itertools来提高速度:
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()
与上面的描述不同,此代码直接计算每个质数从何处开始滚动轮子,以生成其倍数