[PYTHON] 무기한 체에 휠 분해능 추가하기
PYTHON무기한 체에 휠 분해능 추가하기
여기에서 에라 토 스테 네스 (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으로 설정하면 처음 100 개의 소수에만 24133의 정확한 합계를 얻지 만 다른 휠 크기를 사용할 때 누락 된 소수, 잘못된 합계 및 복합체로 끝납니다. 2-3 바퀴 (2와 4의 교체 단계를 사용)조차도 체가 소량을 잃게 만듭니다. 내가 뭘 잘못하고 있죠?
해결법
-
==============================
1."바퀴를 굴리며"[2], 즉 3의 초기 값 (5, 7, 9, ...과 유사하게)에서 시작하여 2의 반복 된 덧셈에 의해 오즈, 즉 2 복사가 생성되고,
"바퀴를 굴리며"[2], 즉 3의 초기 값 (5, 7, 9, ...과 유사하게)에서 시작하여 2의 반복 된 덧셈에 의해 오즈, 즉 2 복사가 생성되고,
n=3; n+=2; n+=2; n+=2; ... # wheel = [2] 3 5 7 92-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 차 요소)입니다.
어디에서 시작할 지 어떻게 알 수 있습니까? 휠 최적화에 대한 요점은 우리가이 일련의 coprimes (이 예에서는 2-3 coprimes)에서만 작업한다는 것입니다. 그래서 재귀 적으로 생성 된 소수를 얻는 코드 부분에서 우리는 롤링 휠 스트림을 유지하고 그 다음 소수를 볼 때까지 계속 진행할 것입니다. 롤링 시퀀스는 두 가지 결과, 즉 값과 휠 위치를 생성해야합니다. 따라서 프라임을 볼 때 해당 휠 위치를 얻습니다. 휠의 해당 위치에서 시작하여 배수를 생성 할 수 있습니다. 우리는 모든 것을 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으로 줄 바꿈).
결과의 소수를 생성하기 위해, 다른 코드 시퀀스를 굴려 참조 코드와 마찬가지로 dict에없는 요소 만 유지합니다.
update : codereview에 대한 몇 번의 반복 작업 (참여자들 덕분에 큰 덕을 보았습니다!) :이 코드를 가능한 한 많이 사용하여 속도에 도달했습니다 :
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()위의 설명과 달리이 코드는 각 자릿수에 대해 바퀴를 굴리는 위치를 직접 계산하여 배수를 생성합니다
-
==============================
2.이것은 제가 생각해 낸 버전입니다. 네스만큼 깨끗하지는 않지만 작동합니다. 휠체어를 사용하는 방법에 대한 또 다른 예가 있습니다. 나는 사용할 바퀴의 크기를 선택할 수있는 능력을 남겨 뒀지 만 좀 더 영구적 인 것을 손쉽게 잡을 수 있습니다. 원하는 크기를 생성하여 코드에 붙여 넣기 만하면됩니다.
이것은 제가 생각해 낸 버전입니다. 네스만큼 깨끗하지는 않지만 작동합니다. 휠체어를 사용하는 방법에 대한 또 다른 예가 있습니다. 나는 사용할 바퀴의 크기를 선택할 수있는 능력을 남겨 뒀지 만 좀 더 영구적 인 것을 손쉽게 잡을 수 있습니다. 원하는 크기를 생성하여 코드에 붙여 넣기 만하면됩니다.
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: http://stackoverflow.com/a/10733621 http://stackoverflow.com/a/19391111""" 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)다음은 1 천만 소수를 생성 할 때 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큰 바퀴도 가능하지만 설치 시간이 길어지는 것을 알 수 있습니다. 바퀴가 커지면 수익이 줄어드는 법칙도 있습니다. 실제로 2-13 휠을 지나치는 것은 그리 어렵지 않습니다. 또한 2-23 휠 (간격 목록에 약 3,600 만 개의 숫자가 있음)을 지나서 메모리 오류가 발생했습니다.
from https://stackoverflow.com/questions/30553925/adding-wheel-factorization-to-an-indefinite-sieve by cc-by-sa and MIT license
'PYTHON' 카테고리의 다른 글
| [PYTHON] Django 전환, 코드 블록을위한 언어 전환, 하나의 언어로 번역 수행 (0) | 2018.11.16 |
|---|---|
| [PYTHON] 파이썬에서 XML을 JSON으로 변환하는 방법은 무엇입니까? [복제] (0) | 2018.11.16 |
| [PYTHON] 파이썬 numpy : datetime64 [ns]를 datetime64 [D] (Numba와 함께 사용)로 변환 할 수 없습니다. (0) | 2018.11.16 |
| [PYTHON] Windows 용 Python C 확장 모듈 빌드하기 (0) | 2018.11.16 |
| [PYTHON] 다른 변수를 기반으로 주문을 보존하여 collect_list (0) | 2018.11.16 |