将轮分解添加到不定筛
Adding wheel factorization to an indefinite sieve
我正在修改 here 的埃拉托色尼无限筛法，因此它使用轮分解来跳过比当前仅检查所有可能性的形式更多的复合材料。
我已经弄清楚如何生成步数以到达轮子上的所有间隙。从那里我想我可以用 +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:
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-coprimes，由"rolling the wheel" [2]产生，即从初始值开始重复加2 3（类似于 5、7、9、...），
n=3; n+=2; n+=2; n+=2; ... # wheel = [2]
3 5 7 9
2-3 互质数通过重复加法 2、4、2、4 等生成：
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（即第 1 个元素）。
我们如何知道从哪里开始？轮优化的要点是我们只处理这个互素数序列（在这个例子中，2-3-互素数）。因此，在我们获得递归生成素数的代码部分，我们还将维护滚轮流，并推进它直到我们在其中看到下一个素数。滚动序列需要产生 两个 结果 - 值和车轮位置。因此，当我们看到素数时，我们也得到了相应的轮子位置，我们可以从轮子上的那个位置开始生成它的倍数。当然，我们将所有内容乘以 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]; ...
因此 dict 中的每个条目都必须保持其当前值、其基本质数和其当前轮子位置（在需要时环绕到 0 以实现循环）。
为了产生结果素数，我们滚动另一个互素数序列，并只保留其中不在字典中的元素，就像在参考代码中一样。
更新： 经过几次迭代 on 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()
不同于上面的描述，这段代码直接计算每个素数从哪里开始转动轮子，生成它的倍数
这是我想出的版本。它不像 ' 那样干净，但它可以工作。我发布了它，所以还有另一个关于如何使用轮子分解的示例，以防万一有人来。我保留了选择要使用的车轮尺寸的能力，但很容易确定一个更永久的尺寸 - 只需生成您想要的尺寸并将其粘贴到代码中即可。
from itertools import count
def wpsieve():
"""prime number generator
call this function instead of roughing or turbo"""
whlSize = 11
initPrms, gaps, c = wheel_setup(whlSize)
for p in initPrms:
yield p
primes = turbo(0, (gaps, c))
for p, x in primes:
yield p
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(int(j - gaps[i-1]))
return steps
def wheel_setup(num):
"builds initial data for sieve"
initPrms = roughing(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 roughing(end):
"finds primes by trial division (roughing pump)"
primes = [2]
for i in range(3, end + 1, 2):
if all(map(lambda x: i%x, primes)):
primes.append(i)
return primes
def turbo(lvl=0, initData=None):
"""postponed prime generator with wheels (turbo pump)
Refs:
"""
gaps, c = initData
yield (c, 0)
compost = {} # found composites to skip
# store as current value: (base prime, wheel index)
ps = turbo(lvl + 1, (gaps, c))
p, x = next(ps)
psq = p*p
gapS = len(gaps) - 1
ix = jx = kx = 0 # indices for cycling the wheel
def cyc(x): return 0 if x > gapS else x # wheel cycler
while True:
c += gaps[ix] # add next step on c's wheel
ix = cyc(ix + 1) # and advance c's index
bp, jx = compost.pop(c, (0,0)) # get base prime and its wheel index
if not bp:
if c < psq: # prime
yield c, ix # emit index for above recursive level
continue
else:
jx = kx # swap indices as a new prime comes up
bp = p
p, kx = next(ps)
psq = p*p
d = c + bp * gaps[jx] # calc new multiple
jx = cyc(jx + 1)
while d in compost:
step = bp * gaps[jx]
jx = cyc(jx + 1)
d += step
compost[d] = (bp, jx)
保留车轮尺寸选项还可以让您了解大车轮的作用有多快。下面是测试代码，用于生成选定尺寸的轮子需要多长时间以及该轮子的筛子有多快。
import time
def speed_test(num, whlSize):
print('-'*50)
t1 = time.time()
initPrms, gaps, c = wheel_setup(whlSize)
t2 = time.time()
print('2-{} wheel'.format(initPrms[-1]))
print('setup time: {} sec.'.format(round(t2 - t1, 5)))
t3 = time.time()
prm = initPrms[:]
primes = turbo(0, (gaps, c))
for p, x in primes:
prm.append(p)
if len(prm) > num:
break
t4 = time.time()
print('run time : {} sec.'.format(len(prm), round(t4 - t3, 5)))
print('prime sum : {}'.format(sum(prm)))
for w in [5, 7, 11, 13, 17, 19, 23, 29]:
speed_test(1e7-1, w)
以下是它在我的计算机上使用 PyPy（Python 2.7 兼容）运行 设置为生成一千万个素数的方式：
2- 3 wheel
setup time: 0.0 sec.
run time : 18.349 sec.
prime sum : 870530414842019
--------------------------------------------------
2- 5 wheel
setup time: 0.001 sec.
run time : 13.993 sec.
prime sum : 870530414842019
--------------------------------------------------
2- 7 wheel
setup time: 0.001 sec.
run time : 7.821 sec.
prime sum : 870530414842019
--------------------------------------------------
2- 11 wheel
setup time: 0.03 sec.
run time : 6.224 sec.
prime sum : 870530414842019
--------------------------------------------------
2- 13 wheel
setup time: 0.011 sec.
run time : 5.624 sec.
prime sum : 870530414842019
--------------------------------------------------
2- 17 wheel
setup time: 0.047 sec.
run time : 5.262 sec.
prime sum : 870530414842019
--------------------------------------------------
2- 19 wheel
setup time: 1.043 sec.
run time : 5.119 sec.
prime sum : 870530414842019
--------------------------------------------------
2- 23 wheel
setup time: 22.685 sec.
run time : 4.634 sec.
prime sum : 870530414842019
更大的轮子是可能的，但你可以看到它们变得相当长，无法安装。随着轮子变大，还有 returns 递减法则 - 越过 2-13 轮子没有多大意义，因为它们并没有真正让它变得更快。我也最终 运行 进入了 2-23 轮的记忆错误（它的 gaps 列表中有大约 3600 万个数字）。
我正在修改 here 的埃拉托色尼无限筛法，因此它使用轮分解来跳过比当前仅检查所有可能性的形式更多的复合材料。
我已经弄清楚如何生成步数以到达轮子上的所有间隙。从那里我想我可以用 +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:
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-coprimes，由"rolling the wheel" [2]产生，即从初始值开始重复加2 3（类似于 5、7、9、...），
n=3; n+=2; n+=2; n+=2; ... # wheel = [2]
3 5 7 9
2-3 互质数通过重复加法 2、4、2、4 等生成：
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（即第 1 个元素）。
我们如何知道从哪里开始？轮优化的要点是我们只处理这个互素数序列（在这个例子中，2-3-互素数）。因此，在我们获得递归生成素数的代码部分，我们还将维护滚轮流，并推进它直到我们在其中看到下一个素数。滚动序列需要产生 两个 结果 - 值和车轮位置。因此，当我们看到素数时，我们也得到了相应的轮子位置，我们可以从轮子上的那个位置开始生成它的倍数。当然，我们将所有内容乘以 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]; ...
因此 dict 中的每个条目都必须保持其当前值、其基本质数和其当前轮子位置（在需要时环绕到 0 以实现循环）。
为了产生结果素数，我们滚动另一个互素数序列，并只保留其中不在字典中的元素，就像在参考代码中一样。
更新： 经过几次迭代 on 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()
不同于上面的描述，这段代码直接计算每个素数从哪里开始转动轮子，生成它的倍数
这是我想出的版本。它不像
from itertools import count
def wpsieve():
"""prime number generator
call this function instead of roughing or turbo"""
whlSize = 11
initPrms, gaps, c = wheel_setup(whlSize)
for p in initPrms:
yield p
primes = turbo(0, (gaps, c))
for p, x in primes:
yield p
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(int(j - gaps[i-1]))
return steps
def wheel_setup(num):
"builds initial data for sieve"
initPrms = roughing(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 roughing(end):
"finds primes by trial division (roughing pump)"
primes = [2]
for i in range(3, end + 1, 2):
if all(map(lambda x: i%x, primes)):
primes.append(i)
return primes
def turbo(lvl=0, initData=None):
"""postponed prime generator with wheels (turbo pump)
Refs:
"""
gaps, c = initData
yield (c, 0)
compost = {} # found composites to skip
# store as current value: (base prime, wheel index)
ps = turbo(lvl + 1, (gaps, c))
p, x = next(ps)
psq = p*p
gapS = len(gaps) - 1
ix = jx = kx = 0 # indices for cycling the wheel
def cyc(x): return 0 if x > gapS else x # wheel cycler
while True:
c += gaps[ix] # add next step on c's wheel
ix = cyc(ix + 1) # and advance c's index
bp, jx = compost.pop(c, (0,0)) # get base prime and its wheel index
if not bp:
if c < psq: # prime
yield c, ix # emit index for above recursive level
continue
else:
jx = kx # swap indices as a new prime comes up
bp = p
p, kx = next(ps)
psq = p*p
d = c + bp * gaps[jx] # calc new multiple
jx = cyc(jx + 1)
while d in compost:
step = bp * gaps[jx]
jx = cyc(jx + 1)
d += step
compost[d] = (bp, jx)
保留车轮尺寸选项还可以让您了解大车轮的作用有多快。下面是测试代码，用于生成选定尺寸的轮子需要多长时间以及该轮子的筛子有多快。
import time
def speed_test(num, whlSize):
print('-'*50)
t1 = time.time()
initPrms, gaps, c = wheel_setup(whlSize)
t2 = time.time()
print('2-{} wheel'.format(initPrms[-1]))
print('setup time: {} sec.'.format(round(t2 - t1, 5)))
t3 = time.time()
prm = initPrms[:]
primes = turbo(0, (gaps, c))
for p, x in primes:
prm.append(p)
if len(prm) > num:
break
t4 = time.time()
print('run time : {} sec.'.format(len(prm), round(t4 - t3, 5)))
print('prime sum : {}'.format(sum(prm)))
for w in [5, 7, 11, 13, 17, 19, 23, 29]:
speed_test(1e7-1, w)
以下是它在我的计算机上使用 PyPy（Python 2.7 兼容）运行 设置为生成一千万个素数的方式：
2- 3 wheel
setup time: 0.0 sec.
run time : 18.349 sec.
prime sum : 870530414842019
--------------------------------------------------
2- 5 wheel
setup time: 0.001 sec.
run time : 13.993 sec.
prime sum : 870530414842019
--------------------------------------------------
2- 7 wheel
setup time: 0.001 sec.
run time : 7.821 sec.
prime sum : 870530414842019
--------------------------------------------------
2- 11 wheel
setup time: 0.03 sec.
run time : 6.224 sec.
prime sum : 870530414842019
--------------------------------------------------
2- 13 wheel
setup time: 0.011 sec.
run time : 5.624 sec.
prime sum : 870530414842019
--------------------------------------------------
2- 17 wheel
setup time: 0.047 sec.
run time : 5.262 sec.
prime sum : 870530414842019
--------------------------------------------------
2- 19 wheel
setup time: 1.043 sec.
run time : 5.119 sec.
prime sum : 870530414842019
--------------------------------------------------
2- 23 wheel
setup time: 22.685 sec.
run time : 4.634 sec.
prime sum : 870530414842019
更大的轮子是可能的，但你可以看到它们变得相当长，无法安装。随着轮子变大，还有 returns 递减法则 - 越过 2-13 轮子没有多大意义，因为它们并没有真正让它变得更快。我也最终 运行 进入了 2-23 轮的记忆错误（它的 gaps 列表中有大约 3600 万个数字）。