Verbesserung des reinen Python-Primzahlensiebs durch die Rekursionsformel

Question

Ich versuche, die Champion-Lösung im Primzahl-Thread weiter zu optimieren, indem ich die komplexe Formel für die Länge der Teilliste herausnehme. len() der gleichen Teilfolge ist zu langsam, da len teuer ist und die Erzeugung der Teilfolge teuer ist. Dies scheint die Funktion etwas zu beschleunigen, aber ich konnte die Division noch nicht entfernen, obwohl ich die Division nur innerhalb der Bedingungsanweisung durchführe. Natürlich könnte ich versuchen, die Längenberechnung zu vereinfachen, indem ich die Optimierung der Startmarkierung für n anstelle von n*n herausnehme...

Ich ersetzte Division / durch Ganzzahldivision //, um mit Python 3 oder

from __future__ import division

Es wäre auch interessant, ob diese Rekursionsformel helfen könnte, die Numpy-Lösung zu beschleunigen, aber ich habe nicht viel Erfahrung mit Numpy.

Wenn Sie Psyco für den Code aktivieren, sieht die Sache jedoch ganz anders aus, und der Atkins-Sieb-Code wird schneller als diese spezielle Slicing-Technik.

import cProfile

def rwh_primes1(n):
    # http://stackoverflow.com/questions/2068372/fastest-way-to-list-all-primes-below-n-in-python/3035188#3035188
    """ Returns  a list of primes < n """
    sieve = [True] * (n//2)
    for i in xrange(3,int(n**0.5)+1,2):
        if sieve[i//2]:
            sieve[i*i//2::i] = [False] * ((n-i*i-1)//(2*i)+1)
    return [2] + [2*i+1 for i in xrange(1,n/2) if sieve[i]]

def primes(n):
    # http://stackoverflow.com/questions/2068372/fastest-way-to-list-all-primes-below-n-in-python/3035188#3035188
    # recurrence formula for length by amount1 and amount2 Tony Veijalainen 2010
    """ Returns  a list of primes < n """
    sieve = [True] * (n//2)
    amount1 = n-10
    amount2 = 6

    for i in xrange(3,int(n**0.5)+1,2):
        if sieve[i//2]:
             ## can you make recurrence formula for whole reciprocal?
            sieve[i*i//2::i] = [False] * (amount1//amount2+1)
        amount1-=4*i+4
        amount2+=4

    return [2] + [2*i+1 for i in xrange(1,n//2) if sieve[i]]

numprimes=1000000
print('Profiling')
cProfile.Profile.bias = 4e-6
for test in (rwh_primes1, primes):
    cProfile.run("test(numprimes)")

Profiling (kein großer Unterschied zwischen den Versionen)

         3 function calls in 0.191 CPU seconds

   Ordered by: standard name

   ncalls  tottime  percall  cumtime  percall filename:lineno(function)
        1    0.006    0.006    0.191    0.191 <string>:1(<module>)
        1    0.185    0.185    0.185    0.185 myprimes.py:3(rwh_primes1)
        1    0.000    0.000    0.000    0.000 {method 'disable' of '_lsprof.Profiler' objects}

         3 function calls in 0.192 CPU seconds

   Ordered by: standard name

   ncalls  tottime  percall  cumtime  percall filename:lineno(function)
        1    0.006    0.006    0.192    0.192 <string>:1(<module>)
        1    0.186    0.186    0.186    0.186 myprimes.py:12(primes)
        1    0.000    0.000    0.000    0.000 {method 'disable' of '_lsprof.Profiler' objects}

Interessanterweise wird durch die Erhöhung des Limits auf 10**8 und das Hinzufügen von Timing-Dekoratoren zu Funktionen das Profiling entfernt:

rwh_primes1 took 23.670 s
primes took 22.792 s
primesieve took 10.850 s

Interessanterweise ist die Zeit um die Hälfte kürzer als bei der Version mit Zahlenliste, wenn man keine Liste von Primzahlen erzeugt, sondern das Sieb selbst zurückgibt.

Answer 1

Sie könnten eine Radoptimierung durchführen. Vielfache von 2 und 3 sind keine Primzahlen, also werden sie nicht gespeichert. Sie könnten dann bei 5 beginnen und Vielfache von 2 und 3 überspringen, indem Sie in Schritten von 2,4,2,4,2,4 usw. erhöhen.

Nachstehend finden Sie einen C++-Code dafür. Hoffentlich hilft das.

void sieve23()
{
    int lim=sqrt(MAX);
    for(int i=5,bit1=0;i<=lim;i+=(bit1?4:2),bit1^=1)
    {
        if(!isComp[i/3])
        {
            for(int j=i,bit2=1;;)
            {
                j+=(bit2?4*i:2*i);
                bit2=!bit2;
                if(j>=MAX)break;
                isComp[j/3]=1;
            }
        }
    }
}

Answer 2

Wenn Sie sich für C++ entscheiden, um die Geschwindigkeit zu erhöhen, habe ich das Python-Sieb nach C++ portiert. Die vollständige Diskussion ist hier zu finden: Portierung des optimierten Siebs von Eratosthenes von Python nach C++ .

Auf Intel Q6600, Ubuntu 10.10, kompiliert mit g++ -O3 und bei N=100000000 dauert dies 415 ms.

#include <vector>
#include <boost/dynamic_bitset.hpp>

// http://vault.embedded.com/98/9802fe2.htm - integer square root
unsigned short isqrt(unsigned long a) {
    unsigned long rem = 0;
    unsigned long root = 0;

    for (short i = 0; i < 16; i++) {
        root <<= 1;
        rem = ((rem << 2) + (a >> 30));
        a <<= 2;
        root++;

        if (root <= rem) {
            rem -= root;
            root++;
        } else root--;

    }

    return static_cast<unsigned short> (root >> 1);
}

// https://stackoverflow.com/questions/2068372/fastest-way-to-list-all-primes-below-n-in-python/3035188#3035188
// https://stackoverflow.com/questions/5293238/porting-optimized-sieve-of-eratosthenes-from-python-to-c/5293492
template <class T>
void primesbelow(T N, std::vector<T> &primes) {
    T i, j, k, sievemax, sievemaxroot;

    sievemax = N/3;
    if ((N % 6) == 2) sievemax++;

    sievemaxroot = isqrt(N)/3;

    boost::dynamic_bitset<> sieve(sievemax);
    sieve.set();
    sieve[0] = 0;

    for (i = 0; i <= sievemaxroot; i++) {
        if (sieve[i]) {
            k = (3*i + 1) | 1;
            for (j = k*k/3; j < sievemax; j += 2*k) sieve[j] = 0;
            for (j = (k*k+4*k-2*k*(i&1))/3; j < sievemax; j += 2*k) sieve[j] = 0;
        }
    }

    primes.push_back(2);
    primes.push_back(3);

    for (i = 0; i < sievemax; i++) {
        if (sieve[i]) primes.push_back((3*i+1)|1);
    }

}