40

エラトステネスの篩と Python 3.1を使用して素数ジェネレーターを作成しました。このコードは、ideone.comで 0.32 秒で正しく適切に実行され、最大 1,000,000 の素数を生成します。

# from bitstring import BitString

def prime_numbers(limit=1000000):
    '''Prime number generator. Yields the series
    2, 3, 5, 7, 11, 13, 17, 19, 23, 29 ...    
    using Sieve of Eratosthenes.
    '''
    yield 2
    sub_limit = int(limit**0.5) 
    flags = [False, False] + [True] * (limit - 2)   
#     flags = BitString(limit)
    # Step through all the odd numbers
    for i in range(3, limit, 2):       
        if flags[i] is False:
#        if flags[i] is True:
            continue
        yield i
        # Exclude further multiples of the current prime number
        if i <= sub_limit:
            for j in range(i*3, limit, i<<1):
                flags[j] = False
#                flags[j] = True

問題は、1,000,000,000 までの数値を生成しようとすると、メモリが不足することです。

    flags = [False, False] + [True] * (limit - 2)   
MemoryError

ご想像のとおり、10 億のブール値 ( Python ではそれぞれ1 バイト4 または 8 バイト (コメントを参照)) を割り当てることは実際には不可能なので、bitstringを調べました。各フラグに 1 ビットを使用すると、メモリ効率が大幅に向上すると考えました。ただし、プログラムのパフォーマンスは大幅に低下しました。素数が 1,000,000 までの場合、実行時間は 24 秒です。これはおそらくビット文字列の内部実装によるものです。

上記のコード スニペットのように、3 行をコメントまたはコメント解除して、BitString を使用するように変更した内容を確認できます。

私の質問は、ビット文字列の有無にかかわらず、プログラムを高速化する方法はありますか?

編集: 投稿する前にコードを自分でテストしてください。当然、既存のコードよりも実行速度が遅い回答は受け入れられません。

もう一度編集します。

私のマシンのベンチマークのリストをまとめました。

4

11 に答える 11

34

あなたのバージョンにはいくつかの小さな最適化があります。if flags[i] is False:True と False の役割を逆にすることで、「 」を「 」に変更できますif flags[i]:。また、2 番目のステートメントの開始値は、代わりに にrangeすることができます。あなたの元のバージョンは、私のシステムでは 0.166 秒かかります。これらの変更により、以下のバージョンは私のシステムで 0.156 秒かかりました。i*ii*3

def prime_numbers(limit=1000000):
    '''Prime number generator. Yields the series
    2, 3, 5, 7, 11, 13, 17, 19, 23, 29 ...
    using Sieve of Eratosthenes.
    '''
    yield 2
    sub_limit = int(limit**0.5)
    flags = [True, True] + [False] * (limit - 2)
    # Step through all the odd numbers
    for i in range(3, limit, 2):
        if flags[i]:
            continue
        yield i
        # Exclude further multiples of the current prime number
        if i <= sub_limit:
            for j in range(i*i, limit, i<<1):
                flags[j] = True

ただし、これはメモリの問題には役立ちません。

C 拡張の世界に移り、開発版のgmpyを使用しました。(免責事項: 私はメンテナーの 1 人です。) 開発バージョンは gmpy2 と呼ばれ、xmpz と呼ばれる変更可能な整数をサポートします。gmpy2 と次のコードを使用すると、実行時間は 0.140 秒になります。1,000,000,000 の制限の実行時間は 158 秒です。

import gmpy2

def prime_numbers(limit=1000000):
    '''Prime number generator. Yields the series
    2, 3, 5, 7, 11, 13, 17, 19, 23, 29 ...
    using Sieve of Eratosthenes.
    '''
    yield 2
    sub_limit = int(limit**0.5)
    # Actual number is 2*bit_position + 1.
    oddnums = gmpy2.xmpz(1)
    current = 0
    while True:
        current += 1
        current = oddnums.bit_scan0(current)
        prime = 2 * current + 1
        if prime > limit:
            break
        yield prime
        # Exclude further multiples of the current prime number
        if prime <= sub_limit:
            for j in range(2*current*(current+1), limit>>1, prime):
                oddnums.bit_set(j)

最適化を推進し、明快さを犠牲にすると、次のコードで 0.107 秒と 123 秒の実行時間が得られます。

import gmpy2

def prime_numbers(limit=1000000):
    '''Prime number generator. Yields the series
    2, 3, 5, 7, 11, 13, 17, 19, 23, 29 ...
    using Sieve of Eratosthenes.
    '''
    yield 2
    sub_limit = int(limit**0.5)
    # Actual number is 2*bit_position + 1.
    oddnums = gmpy2.xmpz(1)
    f_set = oddnums.bit_set
    f_scan0 = oddnums.bit_scan0
    current = 0
    while True:
        current += 1
        current = f_scan0(current)
        prime = 2 * current + 1
        if prime > limit:
            break
        yield prime
        # Exclude further multiples of the current prime number
        if prime <= sub_limit:
            list(map(f_set,range(2*current*(current+1), limit>>1, prime)))

編集: この演習に基づいて、 gmpy2 を accept に変更しましたxmpz.bit_set(iterator)。次のコードを使用すると、1,000,000,000 未満のすべての素数の実行時間は、Python 2.7 で 56 秒、Python 3.2 で 74 秒です。(コメントに記載されているように、xrangeよりも高速ですrange。)

import gmpy2

try:
    range = xrange
except NameError:
    pass

def prime_numbers(limit=1000000):
    '''Prime number generator. Yields the series
    2, 3, 5, 7, 11, 13, 17, 19, 23, 29 ...
    using Sieve of Eratosthenes.
    '''
    yield 2
    sub_limit = int(limit**0.5)
    oddnums = gmpy2.xmpz(1)
    f_scan0 = oddnums.bit_scan0
    current = 0
    while True:
        current += 1
        current = f_scan0(current)
        prime = 2 * current + 1
        if prime > limit:
            break
        yield prime
        if prime <= sub_limit:
            oddnums.bit_set(iter(range(2*current*(current+1), limit>>1, prime)))

編集#2:もう1回試してください!gmpy2 を受け入れるように変更しましたxmpz.bit_set(slice)。次のコードを使用すると、1,000,000,000 未満のすべての素数の実行時間は、Python 2.7 と Python 3.2 の両方で約 40 秒です。

from __future__ import print_function
import time
import gmpy2

def prime_numbers(limit=1000000):
    '''Prime number generator. Yields the series
    2, 3, 5, 7, 11, 13, 17, 19, 23, 29 ...
    using Sieve of Eratosthenes.
    '''
    yield 2
    sub_limit = int(limit**0.5)
    flags = gmpy2.xmpz(1)
    # pre-allocate the total length
    flags.bit_set((limit>>1)+1)
    f_scan0 = flags.bit_scan0
    current = 0
    while True:
        current += 1
        current = f_scan0(current)
        prime = 2 * current + 1
        if prime > limit:
            break
        yield prime
        if prime <= sub_limit:
            flags.bit_set(slice(2*current*(current+1), limit>>1, prime))

start = time.time()
result = list(prime_numbers(1000000000))
print(time.time() - start)

編集 #3: xmpz のビット レベルでのスライスを適切にサポートするように gmpy2 を更新しました。パフォーマンスに変化はありませんが、非常に優れた API です。少し調整して、時間を約 37 秒に短縮しました。(gmpy2 2.0.0b1 の変更点については、編集 #4 を参照してください。)

from __future__ import print_function
import time
import gmpy2

def prime_numbers(limit=1000000):
    '''Prime number generator. Yields the series
    2, 3, 5, 7, 11, 13, 17, 19, 23, 29 ...
    using Sieve of Eratosthenes.
    '''
    sub_limit = int(limit**0.5)
    flags = gmpy2.xmpz(1)
    flags[(limit>>1)+1] = True
    f_scan0 = flags.bit_scan0
    current = 0
    prime = 2
    while prime <= sub_limit:
        yield prime
        current += 1
        current = f_scan0(current)
        prime = 2 * current + 1
        flags[2*current*(current+1):limit>>1:prime] = True
    while prime <= limit:
        yield prime
        current += 1
        current = f_scan0(current)
        prime = 2 * current + 1

start = time.time()
result = list(prime_numbers(1000000000))
print(time.time() - start)

編集 #4: gmpy2 2.0.0b1 で、前の例を壊すいくつかの変更を加えました。gmpy2 は、True を 1 ビットの無限ソースを提供する特別な値として扱わなくなりました。代わりに -1 を使用してください。

from __future__ import print_function
import time
import gmpy2

def prime_numbers(limit=1000000):
    '''Prime number generator. Yields the series
    2, 3, 5, 7, 11, 13, 17, 19, 23, 29 ...
    using Sieve of Eratosthenes.
    '''
    sub_limit = int(limit**0.5)
    flags = gmpy2.xmpz(1)
    flags[(limit>>1)+1] = 1
    f_scan0 = flags.bit_scan0
    current = 0
    prime = 2
    while prime <= sub_limit:
        yield prime
        current += 1
        current = f_scan0(current)
        prime = 2 * current + 1
        flags[2*current*(current+1):limit>>1:prime] = -1
    while prime <= limit:
        yield prime
        current += 1
        current = f_scan0(current)
        prime = 2 * current + 1

start = time.time()
result = list(prime_numbers(1000000000))
print(time.time() - start)

編集 #5: gmpy2 2.0.0b2 にいくつかの拡張を行いました。セットまたはクリアされたすべてのビットを反復できるようになりました。実行時間が ~30% 向上しました。

from __future__ import print_function
import time
import gmpy2

def sieve(limit=1000000):
    '''Returns a generator that yields the prime numbers up to limit.'''

    # Increment by 1 to account for the fact that slices do not include
    # the last index value but we do want to include the last value for
    # calculating a list of primes.
    sieve_limit = gmpy2.isqrt(limit) + 1
    limit += 1

    # Mark bit positions 0 and 1 as not prime.
    bitmap = gmpy2.xmpz(3)

    # Process 2 separately. This allows us to use p+p for the step size
    # when sieving the remaining primes.
    bitmap[4 : limit : 2] = -1

    # Sieve the remaining primes.
    for p in bitmap.iter_clear(3, sieve_limit):
        bitmap[p*p : limit : p+p] = -1

    return bitmap.iter_clear(2, limit)

if __name__ == "__main__":
    start = time.time()
    result = list(sieve(1000000000))
    print(time.time() - start)
    print(len(result))
于 2010-05-25T03:41:49.593 に答える
10

OK、これが私の2番目の答えですが、速度が本質であるため、bitarrayモジュールについて言及する必要があると思いました-それはbitstringの宿敵ですが:)。これは、C の拡張機能である (純粋な Python よりも高速である) だけでなく、スライスの割り当てもサポートしているため、このケースに最適です。ただし、Python 3 ではまだ利用できません。

これを最適化しようとさえしていません。ビット文字列バージョンを書き直しただけです。私のマシンでは、素数が 100 万未満の場合、0.16 秒かかります。

10 億の場合、それは完璧に実行され、2 分 31 秒で完了します。

import bitarray

def prime_bitarray(limit=1000000):
    yield 2
    flags = bitarray.bitarray(limit)
    flags.setall(False)
    sub_limit = int(limit**0.5)
    for i in xrange(3, limit, 2):
        if not flags[i]:
            yield i
            if i <= sub_limit:
                flags[3*i:limit:i*2] = True
于 2010-05-24T20:53:33.490 に答える
8

さて、これは今夜行った(ほぼ完全な)包括的なベンチマークで、どのコードが最も速く実行されるかを確認します。うまくいけば、誰かがこのリストを役に立つと思うでしょう. 私のマシンで完了するのに 30 秒以上かかるものはすべて省略しました。

情報を提供してくれたすべての人に感謝します。私はあなたの努力から多くの洞察を得ました。あなたもそうであることを願っています。

私のマシン: AMD ZM-86、2.40 Ghz デュアルコア、4 GB の RAM。これは HP Touchsmart Tx2 ラップトップです。ペーストビンにリンクしている可能性がありますが、自分のマシンで次のすべてのベンチマークを行ったことに注意してください。

ビルドできるようになったら、gmpy2 ベンチマークを追加します。

すべてのベンチマークは Python 2.6 x86 でテストされています

1,000,000 までの素数 n のリストを返す: ( Python ジェネレーターを使用)

Sebastian の numpy ジェネレーター バージョン (更新) - 121 ミリ秒 @

マークのふるい + ホイール- 154 ミリ秒

スライスを使用した Robert のバージョン- 159 ミリ秒

スライスを使用した改良版 - 205 ミリ秒

列挙型の Numpy ジェネレーター- 249 ミリ秒 @

マークの基本的なふるい- 317 ミリ秒

私の元のソリューションに対する casevh の改善- 343 ミリ秒

私の変更された numpy ジェネレーター ソリューション- 407 ミリ秒

問題の私のオリジナルの方法- 409ミリ秒

Bitarray ソリューション- 414 ms @

bytearray を使用した純粋な Python - 1394 ミリ秒 @

Scott の BitString ソリューション- 6659 ms @

「@」は、このメソッドが私のマシン設定で最大 n < 1,000,000,000 を生成できることを意味します。

さらに、ジェネレーターが必要なく、一度にリスト全体が必要な場合は、次のようにします。

RosettaCode の numpy ソリューション- 32 ミリ秒 @

(numpy ソリューションは、最大 10 億を生成することもできます。これには 61.6259 秒かかりました。メモリが 1 回スワップされたため、2 倍の時間が発生したと思われます。)

于 2010-05-25T10:22:40.187 に答える
6

関連する質問: python で N 以下のすべての素数を一覧表示する最速の方法

こんにちは、Python で10^9までの素数をできるだけ早く生成するコードを探していますが、これはメモリの問題のために困難です。これまでのところ、最大10^6 & 10^7 (私の古いマシンではそれぞれ 0,171 秒と 1,764 秒のクロッキング) の素数を生成するためにこれを思いつきましたが、これは (少なくとも私のコンピューターでは) 「私のおそらく、他のコードでn//i-i +1類似の代わりに (以下を参照) を使用しているためです。len(flags[i2::i<<1])私が間違っている場合は修正してください。改善のための提案は大歓迎です。

def primes(n):
    """ Returns  a list of primes < n """
    sieve = [True] * n
    for i in xrange(3,int(n**0.5)+1,2):
        if sieve[i]:
            sieve[i*i::2*i]=[False]*((n-i*i-1)/(2*i)+1)
    return [2] + [i for i in xrange(3,n,2) if sieve[i]]

彼のコードの 1 つで、Xavier はflags[i2::i<<1]とを使用してlen(flags[i2::i<<1])います。サイズを明示的に計算しました。これは、 の間i*i..n(n-i*i)//2*iの倍数があるという事実を使用して、2*iを数えたいからです。これにより、コードが高速になります。上記のコードに基づく基本的なジェネレーターは次のようになります。i*i1len(sieve[i*i::2*i])(n-i*i)//(2*i) +1

def primesgen(n):
    """ Generates all primes <= n """
    sieve = [True] * n
    yield 2
    for i in xrange(3,int(n**0.5)+1,2):
        if sieve[i]:
            yield i
            sieve[i*i::2*i] = [False]*((n-i*i-1)/(2*i)+1)
    for i in xrange(i+2,n,2):
        if sieve[i]: yield i

少し変更すると、上記のコードのわずかに遅いバージョンを記述できます。これは、ふるいの半分のサイズで開始しsieve = [True] * (n//2)、同じn. それがメモリの問題にどのように反映されるかわかりません。実装の例として、サイズの半分のふるいから始まる numpy rosetta コードの修正バージョン (おそらくより高速) を示します。

import numpy
def primesfrom3to(n):
    """ Returns a array of primes, 3 <= p < n """
    sieve = numpy.ones(n/2, dtype=numpy.bool)
    for i in xrange(3,int(n**0.5)+1,2):
        if sieve[i/2]: sieve[i*i/2::i] = False
    return 2*numpy.nonzero(sieve)[0][1::]+1

より高速でよりメモリに関するジェネレーターは次のようになります。

import numpy
def primesgen1(n):
""" Input n>=6, Generates all primes < n """
sieve1 = numpy.ones(n/6+1, dtype=numpy.bool)
sieve5 = numpy.ones(n/6  , dtype=numpy.bool)
sieve1[0] = False
yield 2; yield 3;
for i in xrange(int(n**0.5)/6+1):
    if sieve1[i]:
        k=6*i+1; yield k;
        sieve1[ ((k*k)/6) ::k] = False
        sieve5[(k*k+4*k)/6::k] = False
    if sieve5[i]:
        k=6*i+5; yield k;
        sieve1[ ((k*k)/6) ::k] = False
        sieve5[(k*k+2*k)/6::k] = False
for i in xrange(i+1,n/6):
        if sieve1[i]: yield 6*i+1
        if sieve5[i]: yield 6*i+5
if n%6 > 1:
    if sieve1[i+1]: yield 6*(i+1)+1

または、もう少しコードを追加します。

import numpy
def primesgen(n):
    """ Input n>=30, Generates all primes < n """
    size = n/30 + 1
    sieve01 = numpy.ones(size, dtype=numpy.bool)
    sieve07 = numpy.ones(size, dtype=numpy.bool)
    sieve11 = numpy.ones(size, dtype=numpy.bool)
    sieve13 = numpy.ones(size, dtype=numpy.bool)
    sieve17 = numpy.ones(size, dtype=numpy.bool)
    sieve19 = numpy.ones(size, dtype=numpy.bool)
    sieve23 = numpy.ones(size, dtype=numpy.bool)
    sieve29 = numpy.ones(size, dtype=numpy.bool)
    sieve01[0] = False
    yield 2; yield 3; yield 5;
    for i in xrange(int(n**0.5)/30+1):
        if sieve01[i]:
            k=30*i+1; yield k;
            sieve01[     (k*k)/30::k] = False
            sieve07[(k*k+ 6*k)/30::k] = False
            sieve11[(k*k+10*k)/30::k] = False
            sieve13[(k*k+12*k)/30::k] = False
            sieve17[(k*k+16*k)/30::k] = False
            sieve19[(k*k+18*k)/30::k] = False
            sieve23[(k*k+22*k)/30::k] = False
            sieve29[(k*k+28*k)/30::k] = False
        if sieve07[i]:
            k=30*i+7; yield k;
            sieve01[(k*k+ 6*k)/30::k] = False
            sieve07[(k*k+24*k)/30::k] = False
            sieve11[(k*k+16*k)/30::k] = False
            sieve13[(k*k+12*k)/30::k] = False
            sieve17[(k*k+ 4*k)/30::k] = False
            sieve19[     (k*k)/30::k] = False
            sieve23[(k*k+22*k)/30::k] = False
            sieve29[(k*k+10*k)/30::k] = False
        if sieve11[i]:
            k=30*i+11; yield k;
            sieve01[     (k*k)/30::k] = False
            sieve07[(k*k+ 6*k)/30::k] = False
            sieve11[(k*k+20*k)/30::k] = False
            sieve13[(k*k+12*k)/30::k] = False
            sieve17[(k*k+26*k)/30::k] = False
            sieve19[(k*k+18*k)/30::k] = False
            sieve23[(k*k+ 2*k)/30::k] = False
            sieve29[(k*k+ 8*k)/30::k] = False
        if sieve13[i]:
            k=30*i+13; yield k;
            sieve01[(k*k+24*k)/30::k] = False
            sieve07[(k*k+ 6*k)/30::k] = False
            sieve11[(k*k+ 4*k)/30::k] = False
            sieve13[(k*k+18*k)/30::k] = False
            sieve17[(k*k+16*k)/30::k] = False
            sieve19[     (k*k)/30::k] = False
            sieve23[(k*k+28*k)/30::k] = False
            sieve29[(k*k+10*k)/30::k] = False
        if sieve17[i]:
            k=30*i+17; yield k;
            sieve01[(k*k+ 6*k)/30::k] = False
            sieve07[(k*k+24*k)/30::k] = False
            sieve11[(k*k+26*k)/30::k] = False
            sieve13[(k*k+12*k)/30::k] = False
            sieve17[(k*k+14*k)/30::k] = False
            sieve19[     (k*k)/30::k] = False
            sieve23[(k*k+ 2*k)/30::k] = False
            sieve29[(k*k+20*k)/30::k] = False
        if sieve19[i]:
            k=30*i+19; yield k;
            sieve01[     (k*k)/30::k] = False
            sieve07[(k*k+24*k)/30::k] = False
            sieve11[(k*k+10*k)/30::k] = False
            sieve13[(k*k+18*k)/30::k] = False
            sieve17[(k*k+ 4*k)/30::k] = False
            sieve19[(k*k+12*k)/30::k] = False
            sieve23[(k*k+28*k)/30::k] = False
            sieve29[(k*k+22*k)/30::k] = False
        if sieve23[i]:
            k=30*i+23; yield k;
            sieve01[(k*k+24*k)/30::k] = False
            sieve07[(k*k+ 6*k)/30::k] = False
            sieve11[(k*k+14*k)/30::k] = False
            sieve13[(k*k+18*k)/30::k] = False
            sieve17[(k*k+26*k)/30::k] = False
            sieve19[     (k*k)/30::k] = False
            sieve23[(k*k+ 8*k)/30::k] = False
            sieve29[(k*k+20*k)/30::k] = False
        if sieve29[i]:
            k=30*i+29; yield k;
            sieve01[     (k*k)/30::k] = False
            sieve07[(k*k+24*k)/30::k] = False
            sieve11[(k*k+20*k)/30::k] = False
            sieve13[(k*k+18*k)/30::k] = False
            sieve17[(k*k+14*k)/30::k] = False
            sieve19[(k*k+12*k)/30::k] = False
            sieve23[(k*k+ 8*k)/30::k] = False
            sieve29[(k*k+ 2*k)/30::k] = False
    for i in xrange(i+1,n/30):
            if sieve01[i]: yield 30*i+1
            if sieve07[i]: yield 30*i+7
            if sieve11[i]: yield 30*i+11
            if sieve13[i]: yield 30*i+13
            if sieve17[i]: yield 30*i+17
            if sieve19[i]: yield 30*i+19
            if sieve23[i]: yield 30*i+23
            if sieve29[i]: yield 30*i+29
    i += 1
    mod30 = n%30
    if mod30 > 1:
        if sieve01[i]: yield 30*i+1
    if mod30 > 7:
        if sieve07[i]: yield 30*i+7
    if mod30 > 11:
        if sieve11[i]: yield 30*i+11
    if mod30 > 13:
        if sieve13[i]: yield 30*i+13
    if mod30 > 17:
        if sieve17[i]: yield 30*i+17
    if mod30 > 19:
        if sieve19[i]: yield 30*i+19
    if mod30 > 23:
        if sieve23[i]: yield 30*i+23

Ps: このコードを高速化する方法について何か提案があれば、操作の順序を変更することから何かを事前に計算することまで、コメントしてください。

于 2010-05-25T21:02:28.383 に答える
4

これは私がしばらく前に書いたバージョンです。速度についてあなたのものと比較するのは興味深いかもしれません。ただし、スペースの問題については何もしません (実際、おそらくあなたのバージョンよりも悪いでしょう)。

from math import sqrt

def basicSieve(n):
    """Given a positive integer n, generate the primes < n."""
    s = [1]*n
    for p in xrange(2, 1+int(sqrt(n-1))):
        if s[p]:
            a = p*p
            s[a::p] = [0] * -((a-n)//p)
    for p in xrange(2, n):
        if s[p]:
            yield p 

ホイールを使用した高速バージョンがありますが、はるかに複雑です。元のソースはこちらです。

さて、これはホイールを使用したバージョンです。 wheelSieveメインのエントリーポイントです。

from math import sqrt
from bisect import bisect_left

def basicSieve(n):
    """Given a positive integer n, generate the primes < n."""
    s = [1]*n
    for p in xrange(2, 1+int(sqrt(n-1))):
        if s[p]:
            a = p*p
            s[a::p] = [0] * -((a-n)//p)
    for p in xrange(2, n):
        if s[p]:
            yield p

class Wheel(object):
    """Class representing a wheel.

    Attributes:
       primelimit -> wheel covers primes < primelimit.
       For example, given a primelimit of 6
       the wheel primes are 2, 3, and 5.
       primes -> list of primes less than primelimit
       size -> product of the primes in primes;  the modulus of the wheel
       units -> list of units modulo size
       phi -> number of units

    """
    def __init__(self, primelimit):
        self.primelimit = primelimit
        self.primes = list(basicSieve(primelimit))

        # compute the size of the wheel
        size = 1
        for p in self.primes:
            size *= p
        self.size = size

        # find the units by sieving
        units = [1] * self.size
        for p in self.primes:
            units[::p] = [0]*(self.size//p)
        self.units = [i for i in xrange(self.size) if units[i]]

        # number of units
        self.phi = len(self.units)

    def to_index(self, n):
        """Compute alpha(n), where alpha is an order preserving map
        from the set of units modulo self.size to the nonnegative integers.

        If n is not a unit, the index of the first unit greater than n
        is given."""
        return bisect_left(self.units, n % self.size) + n // self.size * self.phi

    def from_index(self, i):
        """Inverse of to_index."""

        return self.units[i % self.phi] + i // self.phi * self.size

def wheelSieveInner(n, wheel):
    """Given a positive integer n and a wheel, find the wheel indices of
    all primes that are less than n, and that are units modulo the
    wheel modulus.
    """

    # renaming to avoid the overhead of attribute lookups
    U = wheel.units
    wS = wheel.size
    # inverse of unit map
    UI = dict((u, i) for i, u in enumerate(U))
    nU = len(wheel.units)

    sqroot = 1+int(sqrt(n-1)) # ceiling of square root of n

    # corresponding index (index of next unit up)
    sqrti = bisect_left(U, sqroot % wS) + sqroot//wS*nU
    upper = bisect_left(U, n % wS) + n//wS*nU
    ind2 = bisect_left(U, 2 % wS) + 2//wS*nU

    s = [1]*upper
    for i in xrange(ind2, sqrti):
        if s[i]:
            q = i//nU
            u = U[i%nU]
            p = q*wS+u
            u2 = u*u
            aq, au = (p+u)*q+u2//wS, u2%wS
            wp = p * nU
            for v in U:
                # eliminate entries corresponding to integers
                # congruent to p*v modulo p*wS
                uvr = u*v%wS
                m = aq + (au > uvr)
                bot = (m + (q*v + u*v//wS - m) % p) * nU + UI[uvr]
                s[bot::wp] = [0]*-((bot-upper)//wp)
    return s

def wheelSieve(n, wheel=Wheel(10)):
    """Given a positive integer n, generate the list of primes <= n."""
    n += 1
    wS = wheel.size
    U = wheel.units
    nU = len(wheel.units)
    s = wheelSieveInner(n, wheel)
    return (wheel.primes[:bisect_left(wheel.primes, n)] +
            [p//nU*wS + U[p%nU] for p in xrange(bisect_left(U, 2 % wS)
             + 2//wS*nU, len(s)) if s[p]])

ホイールとは何かについて: (2 を除いて) すべての素数は奇数であるため、ほとんどのふるいはすべての偶数を見逃しています。同様に、もう少し進んで、すべての素数 (2 と 3 を除く) が 1 または 5 モジュロ 6 (== 2 * 3) に一致することに気付くことができるので、それらの数のエントリのみをふるいに保存するだけで済みます。 . 次のステップは、すべての素数 (2、3、および 5 を除く) が 1、7、11、13、17、19、23、29 (モジュロ 30) のいずれかに合同であることに注意することです (ここでは 30 == 2*3 *5) など。

于 2010-05-24T14:57:59.887 に答える
3

ビット文字列を使用して速度を向上させる方法の1つは、現在の数値の倍数を除外するときに「set」メソッドを使用することです。

したがって、重要なセクションは

for i in range(3, limit, 2):
    if flags[i]:
        yield i
        if i <= sub_limit:
            flags.set(1, range(i*3, limit, i*2))    

私のマシンでは、これは元のマシンの約3倍の速度で実行されます。

私の他の考えは、ビット文字列を使用して奇数のみを表すことでした。flags.findall([0])次に、ジェネレータを返すを使用して未設定のビットを見つけることができます。それがはるかに高速であるかどうかはわかりません(ビットパターンを効率的に見つけるのはそれほど簡単ではありません)。

[完全な開示:私はビットストリングモジュールを書いたので、ここでいくつかの誇りを持っています!]


比較として、ビット文字列メソッドから内臓を取り除いて、同じように実行しますが、範囲チェックなどは行いません。これにより、10億個の要素で機能する純粋なPythonの妥当な下限が得られると思います(アルゴリズムを変更しますが、これは不正行為だと思います!)

def prime_pure(limit=1000000):
    yield 2
    flags = bytearray((limit + 7) // 8)
    sub_limit = int(limit**0.5)
    for i in xrange(3, limit, 2):
        byte, bit = divmod(i, 8)
        if not flags[byte] & (128 >> bit):
            yield i
            if i <= sub_limit:
                for j in xrange(i*3, limit, i*2):
                    byte, bit = divmod(j, 8)
                    flags[byte] |= (128 >> bit)

私のマシンでは、これは100万要素に対して約0.62秒で実行されます。これは、ビット配列の回答の約4分の1の速度であることを意味します。それは純粋なPythonにとってはかなり合理的だと思います。

于 2010-05-24T13:53:25.380 に答える
2

検討したいオプションの1つは、C / C ++モジュールをコンパイルするだけなので、ビットをいじる機能に直接アクセスできます。AFAIKは、そのような性質のものを記述し、C /C++で実行された計算の完了時にのみ結果を返すことができます。これを入力しているので、速度のために配列をネイティブCとして格納するnumpyも見ることができます。ただし、それがビット文字列モジュールよりも高速になるかどうかはわかりません。

于 2010-05-24T14:29:45.837 に答える
2

Python の整数型は任意のサイズになる可能性があるため、一連のビットを表現するための巧妙なビット文字列ライブラリは必要ありません。単一の数値だけです。

これがコードです。1,000,000 の制限を簡単に処理し、あまり文句を言わずに 10,000,000 を処理することもできます。

def multiples_of(n, step, limit):
    bits = 1 << n
    old_bits = bits
    max = 1 << limit
    while old_bits < max:
        old_bits = bits
        bits += bits << step
        step *= 2
    return old_bits

def prime_numbers(limit=10000000):
    '''Prime number generator. Yields the series                                                                        
    2, 3, 5, 7, 11, 13, 17, 19, 23, 29 ...                                                                              
    using Sieve of Eratosthenes.                                                                                        
    '''
    yield 2
    sub_limit = int(limit**0.5)
    flags = ((1 << (limit - 2)) - 1) << 2
    # Step through all the odd numbers                                                                                  
    for i in xrange(3, limit, 2):
        if not (flags & (1 << i)):
            continue
        yield i
        # Exclude further multiples of the current prime number                                                         
        if i <= sub_limit:
            flags &= ~multiples_of(i * 3, i << 1, limit)

このmultiples_of関数は、すべての倍数を個別に設定するコストを回避します。最初のビットを設定し、それを最初のステップ ( i << 1) だけシフトし、結果をそれ自体に追加します。次に、ステップを 2 倍にして、次のシフトで両方のビットをコピーするというように、制限に達するまで繰り返します。これは、O(log(limit)) 時間の制限まで、数値のすべての倍数を設定します。

于 2010-05-24T13:44:01.440 に答える
1

エラトステネスのセグメント化されたふるいを使用できます。各セグメントで使用されたメモリは、次のセグメントで再利用されます。

于 2012-02-02T13:58:07.913 に答える
1

別の非常にばかげたオプションですが、非常に高速に利用できる素数の数の大きなリストが本当に必要な場合に役立ちます。たとえば、プロジェクト オイラーの問題に答えるツールとしてそれらが必要な場合 (問題がマインド ゲームとしてコードを最適化するだけの場合は関係ありません)。

遅いソリューションを使用してリストを生成し、それを python ソース ファイルに保存すると、次のprimes.pyものが含まれます。

primes = [ a list of a million primes numbers here ]

それらを使用するにはimport primes、ディスク IO の速度で最小限のメモリ フットプリントでそれらを取得するだけです。複雑さは素数の数です:-)

最適化が不十分なソリューションを使用してこのリストを生成したとしても、それは 1 回しか行われず、あまり問題になりません。

おそらく pickle/unpickle を使用してさらに高速化することもできますが、アイデアはわかります...

于 2010-05-27T00:36:23.097 に答える