Segmentiertes Primzahlsieb

Code-Stücke können hier veröffentlicht werden.
Antworten
oliver1978
User
Beiträge: 1
Registriert: Freitag 22. Dezember 2017, 12:06

Freitag 22. Dezember 2017, 12:43

Seid gegrüßt 8)

ich versuche mich momentan an einem segmentierten Primzahlsieb nach Eratosthenes und habe hierfür zwei Dateien:

(A) __primes.py

Code: Alles auswählen

from __future__ import division
from math import sqrt as _sqrt_
from math import ceil as _ceil_

#
def segmented_sieve(limit):
    # credits to Kim Walisch

    L1CACHE = 32768

    sqrt = int(_sqrt_(limit))
    segment_size = max(sqrt, L1CACHE)

    #if limit < 2:
    #    cnt = 0
    #else:
    #    cnt = 1

    s = 3
    n = 3

    # generate small primes
    is_prime = [1] * (sqrt + 1)
    i = 2
    while i * i <= sqrt:
        if is_prime[i] == 1:
            j = i * i
            while j <= sqrt:
                is_prime[j] = 0
                j += i
        i += 1

    # using lists instead of vectors
    primes = []
    nextp = []
    result = [2]

    for low in range(0, limit + 1, segment_size):
        sieve = [1] * segment_size

        high = min(low + segment_size - 1, limit)

        while s * s <= high:
            if is_prime[s]:
                primes.append(s)
                nextp.append(s * s - low)
            s += 2

        i = 0
        while i < len(primes):
            j = nextp[i]
            k = primes[i] * 2
            while j < segment_size:
                sieve[j] = 0
                j += k
            nextp[i] = j - segment_size
            i += 1

        while n <= high:
            if sieve[n - low] == 1:
                result.append(n - low)
                #cnt += 1
            n += 2

    del (primes)
    del (nextp)
    return result

#
def primes_in_segment(a, b):
    if a > b: a, b = b, a

    primes = segmented_sieve(int(_sqrt_(b)))
    segment = [True] * (b - a + 1)

    for p in primes:
        i = int(p * _ceil_(a / p))
        if (i & 1) == 0: i += p
        while i <= b:
            segment[i-a] = False
            i += p + p

    if (a <= 2) and (b >= 2):
        segprimes = [2]
    else:
        segprimes = []

    i = (a & 1) ^ 1
    while i <= (b - a):
        if segment[i]:
            segprimes.append(i+a)
        i += 2

    del (primes)
    del (segment)
    return segprimes
(B) main.py

Code: Alles auswählen

from __primes import segmented_sieve, primes_in_segment

#
if __name__ == "__main__":

    A = 10**9
    B = A + 5 * 10**6

    P = primes_in_segment(A, B)
    print (len(P), P[:10])
Es soll hier wirklich um Zahlengrößen gehen, wo ich mit einem sequentiellen Sieb abstinke. Beispielhaft suche ich hier die Primzahlen zwischen 1.000.000.000 und 1.005.000.000 und lasse mir die Anzahl der gefundenen Primzahlen und die ersten zehn davon anzeigen.

So weit, so gut. Wenn ich A beispielsweise auf 10**14 erhöhe, oder auch nur auf 10**10, kommt als Ergebnis gar nichts mehr. Ich bekomme dann eine leere Liste mit Länge Null, was sonst.

Ich weiß, die ein oder andere Sache im Quelltext ist nicht ganz pythonisch...aber vielleicht hat jemand eine Idee, wie ich das Problem beheben kann :?:

Ich benutze Python v3.6.0 unter Windows 7 Pro x64.
Sirius3
User
Beiträge: 8637
Registriert: Sonntag 21. Oktober 2012, 17:20

Freitag 22. Dezember 2017, 13:49

@oliver1978: laß das mit den Unterstrichen sein, Wurzel heißt `sqrt` und nicht anders. Zeile 19/20: s und n werden definiert, aber hier gar nicht verwendet. Variableninitalisierung sollte so nah an der Verwendung sein wie möglich. Zeile 23: für boolsche Werte gibt es True und False, nicht 1 und 0. `del` braucht man nicht und kann weg.

Hast Du schon versucht, mit Debugausgaben nachzuvollziehen, was in den einzelnen Schleifen so passiert?
Antworten