Primzahlenberechnung - Sieb des Eratosthenes - Was kann noch optimiert werden?

[manni66]

vector<bool> ist eine Spezialisierung von vector, um pro Bool-Wert auch tatsächlich nur ein Bit zu verwenden. Theoretisch könnte ein vector<char> schneller sein. Praktisch muss man es nachmessen.

[EOutOfResources]

manni66 schrieb:

Theoretisch könnte ein vector<char> schneller sein.

Wegen was? Der Bitmaske?

[SeppJ]

EOutOfResources schrieb:

manni66 schrieb:

Theoretisch könnte ein vector<char> schneller sein.

Wegen was? Der Bitmaske?

Weil vector<bool> eine unnötig komplizierte Spezialisierung ist, während ein vector<char> der übliche schnelle vector wäre.

[SeppJ]

SeppJ schrieb:

EOutOfResources schrieb:

manni66 schrieb:

Theoretisch könnte ein vector<char> schneller sein.

Wegen was? Der Bitmaske?

Weil vector<bool> eine unnötig komplizierte Spezialisierung ist, während ein vector<char> der übliche schnelle vector wäre.

Was noch wichtig ist: Solche Überlegungen müssen nicht stimmen. Man kann zum Beispiel leicht vergessen, dass ein vector<bool> viel bessere Lokalität hat. Wenn der ganze vector (oder große Teile) in den CPU-Cache passen, dann ist das ein gewaltiger Vorteil. Ich hab's mal ausprobiert:

#include <vector>
#include <iostream>
#include <cstddef>

using namespace std;

typedef vector<bool> container;  // Hier entsprechend auf char oder bool setzen

int main()
{
  const size_t max_n = 100000000;
  container numbers(max_n, 1);
  for (unsigned n=2; n<max_n; ++n)
    if (numbers[n])
      for (unsigned i=2*n; i<max_n; i+=n)
        numbers[i]=0;
  unsigned long anti_optimizer = 0;
  for (unsigned n=2; n<max_n; ++n)
    if (numbers[n])
      anti_optimizer += n;
  cout<<"Summe aller Primzahlen bis "<<max_n<<" ist "<<anti_optimizer<<".\n";
}

Die vector<char> Version rechnet bei mir (Core 2, g++, -O3 -xHOST) knapp 8 Sekunden, der vector<bool> etwas über 5 Sekunden. Also ganz klarer Sieg für den kleinen, aber komplizierten, Datentyp.

[volkard]

SeppJ schrieb:

Wenn der ganze vector (oder große Teile) in den CPU-Cache passen, dann ist das ein gewaltiger Vorteil.
...
Die vector<char> Version rechnet bei mir (Core 2, g++, -O3 -xHOST) knapp 8 Sekunden, der vector<bool> etwas über 5 Sekunden. Also ganz klarer Sieg für den kleinen, aber komplizierten, Datentyp.

Jup. Der vector paßt wohl gerade rein? Rein interessehalber ie steht's bei max_n=1e9?

[*Rewind*]

SeppJ schrieb:

Die vector<char> Version rechnet bei mir (Core 2, g++, -O3 -xHOST) knapp 8 Sekunden, der vector<bool> etwas über 5 Sekunden. Also ganz klarer Sieg für den kleinen, aber komplizierten, Datentyp.

ich weiß nicht mit welchem Code du getestet hast, aber ich hab's mal mit deinem probiert (max_n = 100000000;). Ergebnis:

vector<char> container; -> 9s
vector<bool> container; -> 19s

Edit: hab die Angaben korrigiert.

[camper]

Auf meinem Rechner (Atom 330) mit n=100000000:
char: 10.35s
bool: 12.75s

Die Lokalität mit bool ist nicht von besonders großem Vorteil, wenn man bedenkt, dass bool für jeden einzelnen Wert sowohl eine Lese- als auch eine Schreibvorgang impliziert. Wärend mit char nur geschrieben wird - hier könnte man ggf. auch low-level optimieren, z.B. mit unsigned arbeiten und ungeordnet schreiben.

[volkard]

[Rewind] schrieb:

probiert (max_n = 1000000;). Ergebnis:
vector<char> container; -> 2-3s
vector<bool> container; -> 32-33s

Naja, -O3 muß schon sein. Am besten auch -march=native.
So lahm kann Dein Rechner gar nicht sein.

Bei mir:
Bis 1e5
bool: 0.0014
char: 0.0015

Bis 1e6
bool: 0.14
char: 0.18

Bis 1e7
bool: 1.08
char: 1.84

Bis 1e8
bool: 13.77
char: 21.11

Und der Prozessor ist nicht neu. AMD Sempron 3000+ @ 1800MhZ.

[SeppJ]

[Rewind] schrieb:

ich weiß nicht mit welchem Code du getestet hast, aber ich hab's mal mit deinem probiert (max_n = 1000000;). Ergebnis:

vector<char> container; -> 2-3s
vector<bool> container; -> 32-33s

Das widerspricht sich ja auch nicht. Je nach Compiler und System können schon so große Unterschiede auftauchen. Anscheinend hat deine Standardbib einen sehr schlechten vector<bool>, meine einen sehr guten. Darf ich fragen, welcher Compiler?

volkard schrieb:

Jup. Der vector paßt wohl gerade rein? Rein interessehalber ie steht's bei max_n=1e9?

Eigentlich nicht, der Prozessor müsste irgendwie so um 2-4MB Cache haben. Der vector müsste aber mindestens 12.5 MB groß sein.

Abgesehen davon, dass ich gerade an einem anderen Rechner sitze, dürfte mein Core2 Rechner jedoch Probleme mit 1 Milliarde char bekommen, von der RAM-Größe her, daher unvergleichbar.

Ich sitze gerade an einem i7 mit viel mehr RAM und vergleichbar großem Cache, da probiere ich mal (gleicher Compiler wie oben) mit 1 Milliarde:
vector<bool>: gut 24 Sekunden
vector<char>: knapp 26 Sekunden

Sehr interessant. Ich lasse mal im Hintergrund mit 8 Milliarden laufen und melde mich in 10-20 Minuten mit den Ergebnissen.

edit: Ach, dazu habe ich nicht die Geduld. Dafür die Ergebnisse von dem i7 Rechner mit 100 Millionen:
bool: 1.75 Sekunden
char: 2 Sekunden

Da das auch ungefähr das gleiche ist wie bei 1 Milliarde, breche ich mal den 8 Milliarden test ab, da wird sicherlich nichts großarti anderes rauskommen.

[*Rewind*]

Wie bereits in der letzten Antwort dazugeschrieben, die Zeitangaben habe ich korrigiert (9s für char und 19s für bool mit n=100mil). Ich verwende momentan VS2008 PE.

[SeppJ]

Was mir gerade auffällt: Wenn man sich mal Zwischenergebnisse ausgeben lässt, sieht man, dass ein sehr großer Teil der Zeit für das Erstellen des vectors draufgeht. Kein Wunder, dass der vector<bool> so gut abschneidet, denn das geht natürlich deutlich schneller wenn man weniger Speicher braucht.

Auch wieder so ein Aspekt, an den man gar nicht gedacht hätte, ohne dass man es ausprobiert. War wohl doch nicht die Lokalität.

Immer schön vorsichtig sein, bei der Mikrooptimierung.

[camper]

volkard schrieb:

Bis 1e8
bool: 13.77
char: 21.11

Und der Prozessor ist nicht neu. AMD Sempron 3000+ @ 1800MhZ.

Dass ist langsamer als mein atom 330 (1.6GHz) - sollte eigentlich nicht passieren. Sind denn auch alle Checks deaktiviert? Visual C++ hat da ja eine Menge standardmäßig aktiv.

[volkard]

camper schrieb:

volkard schrieb:

Bis 1e8
bool: 13.77
char: 21.11

Und der Prozessor ist nicht neu. AMD Sempron 3000+ @ 1800MhZ.

Dass ist langsamer als mein atom 330 (1.6GHz) - sollte eigentlich nicht passieren. Sind denn auch alle Checks deaktiviert? Visual C++ hat da ja eine Menge standardmäßig aktiv.

Ja. NDEBUG ist definiert, und -O3 -march=native.
gcc 4.5.0 (MinGW).

[Eratosthenes]

SeppJ hat eine wichtige Optimierung vergessen, die im ersten Beispiel
von redrew99 noch drin war: Die Quadratwurzel als Grenze.
Der folgende Code ist noch schneller:

#include <vector>
#include <iostream>
#include <cstddef>
#include <cmath>
using namespace std;
#define NDEBUG

typedef vector<bool> container;  // Hier entsprechend auf char oder bool setzen

int main()
{
   const size_t max_n = 100000000;

   container numbers(max_n, 1);
   size_t grenze = sqrt(max_n); // impliziter cast

   for (unsigned n=2; n<grenze; ++n)     // geaendert!
      if (numbers[n])
         for (unsigned i=2*n; i<max_n; i+=n)
            numbers[i]=0;

   unsigned long anti_optimizer = 0;
   for (unsigned n=2; n<max_n; ++n)
      if (numbers[n])
         anti_optimizer += n;
   cout<<"Summe aller Primzahlen bis "<< max_n
       << " (ausschließlich)  ist "<<anti_optimizer<<".\n";
}

[DocShoe]

Du musst aber schon gegen sqrt( max_n ) +1 laufen (oder zumindest gegen ceil( sqrt( max_n ) ) ), sonst fehlen dir ein paar Überprüfungen.

Edit:
Ist im Original aber drin.

[Michael E.]

DocShoe schrieb:

Du musst aber schon gegen sqrt( max_n ) +1 laufen (oder zumindest gegen ceil( sqrt( max_n ) ) ), sonst fehlen dir ein paar Überprüfungen.

Edit:
Ist im Original aber drin.

Nö, floor(sqrt(max_n)) reicht. Außer wenn der sqrt-Algorithmus Mist baut und die Rechenungenauigkeit dazu führt, dass floor(sqrt(max_n)) ungleich der "richtigen Wurzel" von max_n (abgerundet) ist.

[Jockelx]

Michael E. schrieb:

DocShoe schrieb:

Du musst aber schon gegen sqrt( max_n ) +1 laufen (oder zumindest gegen ceil( sqrt( max_n ) ) ), sonst fehlen dir ein paar Überprüfungen.

Edit:
Ist im Original aber drin.

Nö, floor(sqrt(max_n)) reicht. Außer wenn der sqrt-Algorithmus Mist baut und die Rechenungenauigkeit dazu führt, dass floor(sqrt(max_n)) ungleich der "richtigen Wurzel" von max_n (abgerundet) ist.

Nö, es wird ja auf echt kleiner abgefragt.
Für max_n=5 ist grenze=2 und die Schleife wird überhaupt nicht durchlaufen und 4 somit als Primzahl erkannt.

[Michael E.]

Jockelx schrieb:

Nö, es wird ja auf echt kleiner abgefragt.

Da hast du Recht. Ich hab mir den Code gar nicht durchgelesen. Wer prüft auch schon mit strikt kleiner?

[camper]

Optimierungsmethoden gibt es beim Sieben einige.
Ist etwa das Sieb so groß, dass es nicht mehr in den Cache passt, könnte man in Erwägung ziehen, jeweils nacheinander nur einen Teilbereich zu sieben.
Oder man reduziert die Zahlenmenge die überhaupt betrachtet werden muss.

Alle Primzahlen ausßer 2, 3 haben bei der Division durch 6 entweder den Rest 1 oder 5. Es genügt also, nur diese Restklassen zu betrachten - damit ist die Zahlenmenge schon mal auf ein Drittel reduziert.
Oder man eliminiert auch noch die 5 und betrachtet nur noch die Restklassen
{1,7,11,13,17,19,23,29} mod 30 - das ist dann nur noch 4/15 der Ausgangsmenge, usw.

Jede Restklasse kann separat gesiebt werden, diese Einzelsiebe sind ganz nebenbei kleiner, was das Cacheproblem reduziert. Zudem biete sich hier eine Parallelisierung an.

Schließlich kann man noch den Siebalgorithmus selbst ändern, wenn man z.B. das Sieb von Atkins benutzt.

Ich habe mal die 2. Methode implementiert. Der Code ist relativ schnell hingeschrieben, deshalb nicht sehr leserlich und wahrscheinlich auch nicht besonders optimal.

#include <vector>
#include <cmath>
#include <cassert>
#include <algorithm>
#include <functional>
#include <numeric>
#include <iostream>
#include <ctime>
#include <boost/utility.hpp>
#include <boost/tuple/tuple.hpp>
#include <utility>
#include <boost/lexical_cast.hpp>

using std::cout;
using std::vector;
using std::sqrt;
using std::accumulate;
using std::multiplies;
using std::make_pair;
using std::pair;

typedef unsigned long long number_type;
typedef vector<char> sieve_type;
number_type limit = 0;

pair<long long, long long> euclid(number_type a, number_type b)
{
    if ( b == 0 )
        return make_pair( 1, 0 );
    long long q = a / b;
    long long s, t;
    boost::tie( s, t ) = euclid( b, a % b );
    return make_pair( t, s - q * t );
}

vector<pair<number_type, number_type> > simple_sieve(number_type limit, number_type p)
{
    assert( limit > 2 );
    vector<pair<number_type,number_type> > result;
    result.push_back( make_pair( 2, 0 ) );
    sieve_type sieve( limit );
    for ( number_type i = 3; i * i < limit; i += 2 )
        if ( sieve[ i ] == 0 )
            for ( number_type j = i * i; j < limit; j += 2 * i )
                sieve[ j ] = 1;
    for ( number_type i = 3; i < limit; i += 2 )
        if ( sieve[ i ] == 0 )
        {
            long long t = euclid( p, i % p ).second;
            if ( t < 0 )
                t = p - ( -t % p );
            result.push_back( make_pair( i, t ) );
        }
    return result;
}

vector<number_type> candidate_residues(const vector<number_type>& filtered_primes, number_type size)
{
    sieve_type sieve( size );
    for ( vector<number_type>::const_iterator it = filtered_primes.begin(); it != filtered_primes.end(); ++it )
        for ( number_type i = *it; i < size; i += *it )
            sieve[ i ] = 1;
    vector<number_type> result;
    for ( number_type i = 1; i < size; ++i )
        if ( sieve[ i ] == 0 )
            result.push_back( i );
    return result;
}

template<typename func>
void residue_sieve(sieve_type& sieve, number_type residue, number_type size, const vector<pair<number_type,number_type> >& small_primes, number_type filtered, func f)
{
    number_type limit = sieve.size();
    sieve.assign( limit, 0 );
    for ( vector<pair<number_type,number_type> >::const_iterator it = small_primes.begin() + filtered; it != small_primes.end(); ++it )
    {
        number_type v = it->first;
        number_type w = residue * it->second % size;
        number_type x = ( v * w / size ) + ( w < v ? v : 0);
        for ( ; x < limit; x += v )
            sieve[ x ] = 1;
    }
    number_type i = residue == 1 ? 1 : 0;
    for ( ; i < limit; ++i )
        if ( sieve[ i ] == 0 )
            f( i * size + residue );
}

number_type sum = 0;
void add_prime(number_type v)
{
    if ( v < limit )
    {
        sum += v;
    }
}

int main(int argc, char** argv)
{
    limit = 1 + ( argc > 1 ? boost::lexical_cast<number_type>(argv[1]) : 1000 );
    number_type sieve_limit = argc > 2 ? boost::lexical_cast<number_type>(argv[2]) : 131072;
    std::time_t t = std::clock();
    number_type filtered[] = { 2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61, 67, 71, 73, 79, 83, 89, 97 };
    number_type num_filtered = 0;
    number_type p = 1;
    do
    {
        p *= filtered[ num_filtered++ ];
    } while ( limit / p >= sieve_limit );
    vector<number_type> filtered_primes( filtered, filtered + num_filtered );
    number_type small_limit = sqrt( limit - 1 ) + 1;
    vector<pair<number_type,number_type> > small_primes = simple_sieve( small_limit, p );
    vector<number_type> residues = candidate_residues( filtered_primes, p );
    for ( vector< number_type >::const_iterator it = filtered_primes.begin(); it != filtered_primes.end(); ++it )
        add_prime( *it );
    sieve_type sieve( ( limit + p - 1 ) / p, 0 );
    for ( vector<number_type>::const_iterator it = residues.begin(); it != residues.end(); ++it )
        residue_sieve( sieve, *it, p, small_primes, filtered_primes.size(), add_prime );
    cout << "\nlimit: " << limit-1 << "\tsum: " << sum << '\n';
    cout << "Time: " << (std::clock()-t)/(double)CLOCKS_PER_SEC << '\n';
    cout << "sieve-Size: " << sieve.size() << '\n';
    cout << "p: " << p << '\n';
    cout << "residues: " << residues.size() << '\n';
}

Zeit für Primzahlen bis 1e9 (Atom 330): 6.64s
Zeit mit dem Code von Eratosthenes: 77.46s
Also eine Steigerung um den Faktor: 11

[Eratosthenes]

for (unsigned n=2; n<= grenze; ++n)  ...

ist natürlich richtig.

@Camper: das ist zwar recht viel Code, aber leider falsch. So ist die Summe aller Primzahlen bis 10 nur 17 (=2+3+5+7). Dein Programm spuckt aber beim Aufruf

a.out 10

folgendes aus:

limit: 10       sum: 23

Eine Summe von 23 kann überhaupt bei der Summation von Primzahlen nicht vorkommen.