program story

Python에서 소수의 효율적인 무한 생성기를 구현하는 방법은 무엇입니까?

inputbox 2020. 12. 6. 21:32
반응형

Python에서 소수의 효율적인 무한 생성기를 구현하는 방법은 무엇입니까?


이것은 숙제가 아닙니다. 저는 단지 궁금합니다.

여기에서 INFINITE가 핵심 단어입니다.

나는 그것을 for p in primes(). 나는 이것이 Haskell의 내장 함수라고 믿습니다.

따라서 대답은 "그냥 체를하세요"만큼 순진 할 수 없습니다.

우선, 얼마나 많은 연속 소수가 소비 될지 모릅니다. 한 번에 100 개를 만들 수 있다고 가정 해 보겠습니다. 소수 공식의 빈도와 동일한 Sieve 접근 방식을 사용 하시겠습니까?

비 동시 접근 방식을 선호합니다.

읽어 주셔서 감사합니다 (그리고 쓰기;))!


"내가 더 본다면…"

erat2요리 책 기능은 더 빨라질 수 있습니다 (약 20-25 %).

erat2a

import itertools as it
def erat2a( ):
    D = {  }
    yield 2
    for q in it.islice(it.count(3), 0, None, 2):
        p = D.pop(q, None)
        if p is None:
            D[q*q] = q
            yield q
        else:
            # old code here:
            # x = p + q
            # while x in D or not (x&1):
            #     x += p
            # changed into:
            x = q + 2*p
            while x in D:
                x += 2*p
            D[x] = p

not (x&1)체크를 검증 x홀수이다. 그러나 둘 다 홀수 이기 때문에 단계의 절반을 추가하면 이상한 테스트와 함께 방지됩니다.qp2*p

erat3

약간의 erat2공상에 신경 쓰지 않는다면 다음 변경 사항으로 35-40 % 속도를 높일 수 있습니다 (주의 : itertools.compress함수 때문에 Python 2.7+ 또는 Python 3+ 필요 ).

import itertools as it
def erat3( ):
    D = { 9: 3, 25: 5 }
    yield 2
    yield 3
    yield 5
    MASK= 1, 0, 1, 1, 0, 1, 1, 0, 1, 0, 0, 1, 1, 0, 0,
    MODULOS= frozenset( (1, 7, 11, 13, 17, 19, 23, 29) )

    for q in it.compress(
            it.islice(it.count(7), 0, None, 2),
            it.cycle(MASK)):
        p = D.pop(q, None)
        if p is None:
            D[q*q] = q
            yield q
        else:
            x = q + 2*p
            while x in D or (x%30) not in MODULOS:
                x += 2*p
            D[x] = p

erat3함수는 모든 소수 (2, 3, 5 제외) 모듈로 30이 결과로 8 개의 숫자, 즉 MODULOSfrozenset에 포함 된 숫자 만 얻는다는 사실을 활용합니다. 따라서 초기 세 개의 소수를 산출 한 후 7부터 시작 하여 후보자 사용합니다.
후보 필터링은 itertools.compress함수를 사용합니다 . "마법"은 MASK순서에 있습니다. 7 개부터 시작하여 가능한 모든 후보에 대해 MASK15 개 요소 ( itertools.islice함수에서 선택한 30 개 숫자마다 15 개의 홀수 )가 1있습니다.주기는 itertools.cycle함수에서 지정한대로 반복됩니다 .
후보 필터링의 도입에는 or (x%30) not in MODULOS검사 라는 또 다른 수정이 필요합니다 . 그만큼erat2알고리즘은 모든 홀수를 처리했습니다. 이제 erat3알고리즘이 r30 후보 만 처리하므로 모든 D.keys()후보가 이러한 "거짓"후보 만 될 수 있는지 확인해야합니다 .

벤치 마크

결과

Atom 330 Ubuntu 9.10 서버, 버전 2.6.4 및 3.1.1+ :

$ testit
up to 8192
==== python2 erat2 ====
100 loops, best of 3: 18.6 msec per loop
==== python2 erat2a ====
100 loops, best of 3: 14.5 msec per loop
==== python2 erat3 ====
Traceback (most recent call last):
AttributeError: 'module' object has no attribute 'compress'
==== python3 erat2 ====
100 loops, best of 3: 19.2 msec per loop
==== python3 erat2a ====
100 loops, best of 3: 14.1 msec per loop
==== python3 erat3 ====
100 loops, best of 3: 11.7 msec per loop

AMD Geode LX Gentoo 홈 서버, Python 2.6.5 및 3.1.2 :

$ testit
up to 8192
==== python2 erat2 ====
10 loops, best of 3: 104 msec per loop
==== python2 erat2a ====
10 loops, best of 3: 81 msec per loop
==== python2 erat3 ====
Traceback (most recent call last):
AttributeError: 'module' object has no attribute 'compress'
==== python3 erat2 ====
10 loops, best of 3: 116 msec per loop
==== python3 erat2a ====
10 loops, best of 3: 82 msec per loop
==== python3 erat3 ====
10 loops, best of 3: 66 msec per loop

벤치 마크 코드

primegen.py모듈을 포함 erat2, erat2aerat3기능. 다음은 테스트 스크립트입니다.

#!/bin/sh
max_num=${1:-8192}
echo up to $max_num
for python_version in python2 python3
do
    for function in erat2 erat2a erat3
    do
        echo "==== $python_version $function ===="
        $python_version -O -m timeit -c \
        -s  "import itertools as it, functools as ft, operator as op, primegen; cmp= ft.partial(op.ge, $max_num)" \
            "next(it.dropwhile(cmp, primegen.$function()))"
    done
done

영업 이익은 요청 때문에 효율적인 구현, 여기에 상당한 개선의 활성 상태로 2002 코드 (여기에 볼 데이빗 엡스타인 / 알렉스 마르 텔리에 의해 그의 대답은 ) 는 사각이 중 볼 때까지 사전에서 소수의 정보를 기록하지 않습니다 후보자 . n 개의 소수가 생성 된 경우 공간 복잡도 O (n ) 대신 O (sqrt (n)) 이하로 낮 춥니 다 ( π (sqrt (n log n)) ~ 2 sqrt (n log n) / log (n log n) ~ 2 sqrt (n / log n) ). 결과적으로 시간 복잡성도 개선 됩니다 . , 더 빠르게 실행됩니다 .

단계과 함께 각 기본 프라임 (즉, 현재 생산 지점의 sqrt 아래)의 현재 배수 사전으로 "슬라이딩 체"를 만듭니다 .

from itertools import count
                                         # ideone.com/aVndFM
def postponed_sieve():                   # postponed sieve, by Will Ness      
    yield 2; yield 3; yield 5; yield 7;  # original code David Eppstein, 
    sieve = {}                           #   Alex Martelli, ActiveState Recipe 2002
    ps = postponed_sieve()               # a separate base Primes Supply:
    p = next(ps) and next(ps)            # (3) a Prime to add to dict
    q = p*p                              # (9) its sQuare 
    for c in count(9,2):                 # the Candidate
        if c in sieve:               # c's a multiple of some base prime
            s = sieve.pop(c)         #     i.e. a composite ; or
        elif c < q:  
             yield c                 # a prime
             continue              
        else:   # (c==q):            # or the next base prime's square:
            s=count(q+2*p,2*p)       #    (9+6, by 6 : 15,21,27,33,...)
            p=next(ps)               #    (5)
            q=p*p                    #    (25)
        for m in s:                  # the next multiple 
            if m not in sieve:       # no duplicates
                break
        sieve[m] = s                 # original test entry: ideone.com/WFv4f

(에서 볼 수 있듯이 여기에 이전, 원래의 코드는 법인 (法人)의 변화에 편집 한 대답 으로 팀 피터스 (Tim Peters) 아래). 관련 토론 이것을 참조하십시오 .

유사한 2-3-5-7 휠 기반 코드 가 ~ 2.15 배 더 빠르게 실행됩니다 (이론적 개선에 매우 근접 함 3/2 * 5/4 * 7/6 = 2.1875).


후손을 위해 다음은 Will Ness의 Python 3 용 아름다운 알고리즘의 재 작성입니다. 일부 변경이 필요합니다 (반복자에는 더 이상 .next()메소드가 없지만 새로운 next()내장 함수가 있습니다). 다른 변경 사항은 재미를위한 것입니다 (새로운 내용을 사용하면 원본의 yield from <iterable>4 개 yield문이 대체 됩니다. 더 많은 내용은 가독성을위한 것입니다 (나는 ;-) 한 글자 변수 이름을 남용하는 것을 좋아하지 않습니다).

원본보다 훨씬 빠르지 만 알고리즘상의 이유로는 아닙니다. 속도 향상은 대부분 원본의 add()기능 을 제거하고 대신 인라인으로 수행하기 때문입니다.

def psieve():
    import itertools
    yield from (2, 3, 5, 7)
    D = {}
    ps = psieve()
    next(ps)
    p = next(ps)
    assert p == 3
    psq = p*p
    for i in itertools.count(9, 2):
        if i in D:      # composite
            step = D.pop(i)
        elif i < psq:   # prime
            yield i
            continue
        else:           # composite, = p*p
            assert i == psq
            step = 2*p
            p = next(ps)
            psq = p*p
        i += step
        while i in D:
            i += step
        D[i] = step

이것은 원래 내 코드가 아니지만 게시 할 가치가 있습니다. 원본은 여기에서 찾을 수 있습니다 : http://code.activestate.com/recipes/117119/

def gen_primes():
  D = {}
  q = 2  # first integer to test for primality.

  while True:
    if q not in D:
      # not marked composite, must be prime  
      yield q 

      #first multiple of q not already marked
      D[q * q] = [q] 
    else:
      for p in D[q]:
        D.setdefault(p + q, []).append(p)
      # no longer need D[q], free memory
      del D[q]

    q += 1

발전기이므로 다른 것과 같이 사용하십시오.

primes = gen_primes()
for p in primes:
  print p

내 데스크톱에서 1 백만 개의 소수를 생성하여 세트에 넣는 데 1.62 초가 걸립니다.


DO가 분할 된 세그먼트의 크기가 사용할 수있는 메모리 또는 비트 세트의 최대 크기에 의해 결정된다 체.

각 세그먼트에 대해 일정한 간격 [n; n + segment_size)를 비트 세트로 설정하고 상한의 제곱근 아래의 모든 소수를 체질합니다.

조밀 한 숫자 세트로 작업하기 때문에 비트 세트를 사용하면 해시 테이블이나 트리 데이터 구조보다 적은 메모리를 사용합니다.


이를 수행하는 또 다른 방법 :

import itertools
def primeseq():
    prime = [2]
    num = 0
    yield 2
    for i in itertools.count(3, 2):
        is_prime = True
        for num in prime:
            if i % num == 0:
                is_prime = False
                break
            elif num ** 2 > i: 
                break
        if is_prime:
            prime.append(i)
            yield i

그리고 또 다른 대답은 erat3여기 대답 보다 메모리 효율적입니다 .

import heapq

def heapprimegen():
    hp= []
    yield 2
    yield 3
    cn= 3
    nn, inc= 3, 6
    while 1:
        while cn < nn:
            yield cn
            heapq.heappush(hp, (3*cn, 2*cn))
            cn+= 2
        cn= nn+2
        nn, inc= heapq.heappushpop(hp, (nn+inc, inc))

사전이 아닌 소수 배수의 힙 (목록)을 유지합니다. 분명히 속도를 잃습니다.


다음은 Python2로 작성되었지만 Python3로 쉽게 조정되는 매우 빠른 무한 생성기입니다. 이를 사용하여 10 ** 9까지 소수를 추가하려면 다음을 사용하십시오.

from itertools import takewhile
from functools import partial
from operator import gt
print (sum(takewhile(partial(gt, 10**9), prime_gen_inf())))

분할 된 체로, 빠르지 만 Will Ness의 알고리즘보다 분명히 덜 우아합니다.

from operator import mul
from functools import reduce
def prod(x): return reduce(mul, x, 1)


def build_sieve(wheel):
    w = prod(wheel)
    w_phi = prod([p-1 for p in wheel])
    rems = [a for a in range(w) if all(a % p for p in wheel)]
    assert len(rems) == w_phi
    inv = {a:pow(a, w_phi - 1, w) for a in rems}
    try:
        known_p = wheel + rems[1 : rems.index(rems[1]*rems[1])]
    except ValueError:
        known_p = wheel + rems[1:]
    return wheel, w, w_phi, rems, inv, known_p

#Adjust the chunk variable based on your computer's architecture.
#
#Adjust the line with #! if you don't need "true" infinite.  If you don't need
#primes larger than 1<<32, use array('H', []), if 1<<64 use 'L', if 1<<128 (in
#Python3) use 'Q', otherwise use empty list [].
#To save memory, comment out the lines with #*, and uncomment the commented-out
#lines 
import itertools
from itertools import islice, count, compress, izip
chain_f = itertools.chain.from_iterable
from array import array
def prime_gen_inf(chunk=250000, sieve_info = build_sieve([2,3,5,7])):
    """    Indefinitely yields primes    """
    wheel, w, w_phi, rems, inv, known_p = sieve_info
    for p in known_p: yield p
    new_n = 0;
    while True:
        size = min(chunk, (p * p - new_n) / w)
        sieve = bytearray([1]) * size * w_phi
        n, new_n = new_n, new_n + size * w
        if not n:
            zero = bytearray([0])
            seen = len(known_p) - len(wheel) + 1
            sieve[:seen:1] = zero * seen
            p_gen = islice(prime_gen_inf(), len(wheel), None)
            new_p = next(p_gen)
            ps = []                                         #! array('H', [])
            p_invs = bytearray([])                                         #*
        while new_p * new_p < new_n:
            ps.append(new_p)
            p_invs.append(inv[new_p % w])                                  #*
            new_p = next(p_gen)
        for p, p_inv, modp in izip(ps, p_invs, [-n % p for p in ps]):      #*
            s = [(modp + p * (p_inv * (r - modp) % w)) / w for r in rems]  #*
        #for p in ps:
        #    s = [(-n%p + p * (inv[p%w] * (r - -n%p) % w)) / w for r in rems]
            for i, start in enumerate(s):
                slice_size = ((size - start - 1) / p + 1)
                sieve[i + start * w_phi :: p * w_phi] = zero * slice_size
        for p in compress(chain_f(izip(*[count(n+r, w) for r in rems])), sieve):
            yield p

다음은 Haskell에서 수행되는 방식에 대해 좀 더 사실적인 생성기입니다. 알려진 소수의 합성물에 대해 필터링 한 다음 나머지 소수를 목록에 추가합니다.

def gen_primes():
    primes = []
    i = 2
    while True:
        prime = True
        for p in primes:
            if not (i % p):
                prime = False
                break
        if prime:
            yield i
            primes.append(i)
        i += 1

몇 번 전에 무한 소수 생성기에 대한 기사를 썼습니다.

http://stacktrace.it/2008/01/progetto-eulero-problema-3/

이탈리아어로되어 있지만 Google을 사용하면 성가신 번역이있을 수 있습니다. http://tinyurl.com/yzpyeom


다음은 복잡한 힙 기반 구현으로, 다른 힙 기반 구현보다 훨씬 빠르지는 않지만 (다른 답변에서 속도 비교 참조) 훨씬 적은 메모리를 사용합니다.

이 구현은 동일한 숫자 요소를 포함하는 두 개의 힙 (tu 및 wv)을 ​​사용합니다. 각 요소는 int 쌍입니다. 모든 소수를 발견하기 위해 q**2(여기서 q소수이다), 각 힙 대부분에 포함 할 2*pi(q-1)요소 pi(x)보다는 긍정적 인 소수의 수가 크지 않다 x. 따라서 총 정수 수는 최대 4*pi(floor(sqrt(n)))입니다. (우리는 힙에 절반만큼 많은 것을 푸시함으로써 메모리에서 2의 요소를 얻을 수 있지만 알고리즘이 느려질 것입니다.)

다른 dict 및 heap 기반 접근 방식 (예 : erat2b, heap_prime_gen_squares 및 heapprimegen)은 약`2 * pi (n) '정수를 저장합니다. 왜냐하면 소수를 찾을 때마다 힙 또는 dict를 확장하기 때문입니다. 비교로 1_000_000 개의 소수를 찾기 위해이 구현은 4141 개 미만의 정수를 저장하고 다른 구현은 1_000_000 개 이상의 정수를 저장합니다.

import heapq

def heap_prime_gen_smallmem():
    yield 2
    yield 3
    f = 5
    fmar3 = 2
    q = 7
    q6 = 7 * 6
    qmar3 = 4
    tu = [(25, 30), (35, 30)]
    vw = [(25, 30), (35, 30)]
    while True:
        qmar3 += 2   
        if qmar3 == 6:  
            qb = q + 4
            q6b = q6 + 24
            qmar3 = 2
        else:
            qb = q + 2
            q6b = q6 + 12
        if q < tu[0][0]:
            d = q * q
            while f < d:
                a, b = vw[0]
                if f < a: 
                    yield f   
                else:
                    a, b = vw[0]
                    heapq.heapreplace(vw, (a + b, b))
                    a, b = vw[0]
                    while f >= a:
                        heapq.heapreplace(vw, (a + b, b))
                        a, b = vw[0]   
                fmar3 += 2
                if fmar3 == 6:
                    f += 4
                    fmar3 = 2
                else:
                    f += 2
            c = q * qb   
            heapq.heappush(tu, (d, q6))
            heapq.heappush(tu, (c, q6))
            heapq.heappush(vw, (d, q6))
            heapq.heappush(vw, (c, q6))
        else:
            a, b = tu[0]
            heapq.heapreplace(tu, (a + b, b))
            a, b = tu[0]  
            while q >= a:
                heapq.heapreplace(tu, (a + b, b))
                a, b = tu[0]
        q = qb
        q6 = q6b

다음은 dict 대신 힙을 사용하는 간단하지만 느리지는 않습니다.

import heapq

def heap_prime_gen_squares(): 
    yield 2  
    yield 3  
    h = [(9, 6)]
    n = 5
    while True:
        a, b = h[0]
        while n < a:
            yield n
            heapq.heappush(h, (n * n, n << 1))
            n += 2
        heapq.heapreplace(h, (a + b, b))  # Replace h[0], which is still (a, b).

처음 1 백만 소수에 대한 사용자 시간의 속도 측정 값 (작을수록 좋습니다) :

  • postponed_sieve (dict 기반) : 8.553 초
  • erat2b (dict 기반) : 9.513s
  • erat2a (dict 기반) : 10.313s
  • heap_prime_gen_smallmem (힙 기반) : 23.935s
  • heap_prime_gen_squares (힙 기반) : 27.302s
  • heapprimegen (dict 기반) : 145.029s

따라서 dict 기반 접근 방식이 가장 빠른 것 같습니다.


I know the post is old, but I came by myself across this question... The following code is based on a very simple idea: a growing sieve of Eratosthenes. This solution is really slower than the best ones here, but it is easy to grasp and designed to be readable...

I used integers to store the results of the sieve. In binary format, an integer is a list of 0s and 1s, 0 at position i if i is not a prime, 1 if it may be a prime. The requisite infinity is a result of the fact that Python 3 integers are unbounded.

def primes():
    container, size = 1 << 2, 3 # we start with 0b100 (from right to left: 0 and 1 are not primes, 2 is
    last_prime = 1
    while True:
        prime = next((j for j in range(last_prime+1, size) if container & 1 << j), None) # find the next prime
        while not prime:
            container, size = expand(container, size, 2**16) # add 65536 cells and sieve the container
            prime = next((j for j in range(last_prime+1, size) if container & 1 << j), None)
        yield prime
    last_prime = prime

How to expand the container? Just add a bunch of 1s at the left of the container (in binary format) and sieve them. This is identical to the standard sieve, with a slight difference. In the standard sieve, if we find a prime i, we start to cross the cells at i*i, with a step of i.

Here, this may have been done for the first part of container. We just need to start at the beginnig of the new part of the container if it is farther than i*i.

def expand(container, size, n):
    new_size = size + n
    container += (1 << (new_size + 1) - 1) - (1 << size) # add n 1's
    for i in range(2, new_size):
        if container & (1 << i): # i is a prime
            t = sum(1 << j for j in range(max(i, size // i)*i, new_size, i)) # set 1 for all mutiple
            container &= ~t # cross the cells

    return container, new_size

Test for a million primes:

import itertools
assert 78498 == len(list(itertools.takewhile(lambda p: p<1000000, primes())))

참고URL : https://stackoverflow.com/questions/2211990/how-to-implement-an-efficient-infinite-generator-of-prime-numbers-in-python

반응형