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

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

  • Mod

    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.

    0
  • C
    Mod

    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.

    0
  • V

    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).

    0
  • ?

    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";
}
    0

  • 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.

    0
  • M

    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.

    0
  • J

    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.

    0
  • M

    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? 🤡

    0
  • C
    Mod

    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

    0
  • ?
     

    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.

    0
  • C
    Mod

    Eratosthenes schrieb:

    @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.

    Grenzen < 16 werden offenbar nicht korrekt behandelt, diese Fälle sind allerdings auch ziemlich uninteressant.

    Ich bin noch dabei, das Ganze etwas eleganter zu schreiben.

    0
  • ?

    Wenn das uninteressant ist - gut. Ich meine nur, dass ein korrektes Programm besser ist als ein schnelles. Weiteres Beispiel:
    Die Summe der Primzahlen bis 540 ist 23592. Dein Programm liefert aber 24121. Auch 24121 kommt in der Folge der Summen nicht vor. Interessant ist auch limit = 500. Da liefert dein Programm mal das richtige Ergebnis 21536. An dem niedrigen Limit < 16 liegt es dann wohl doch nicht.

    0
  • C
    Mod

    Code (2 bugs) korrigiert. Probleme traten auf, wenn der ganzzahlige Rest der Wurzel der Grenze zufällig eine Primzahl war.

    0
  • R

    Das mit der Summenbildung ist interessant. Ich habe mal versucht, daß in mein Programm mit einzubauen, also eine zusätzliche Variable deklariert/initialisiert und unten in der Schleife mit eingebaut, also

    unsigned long long int x=0;

    und

    for(unsigned long long int j=2;j<vb.size();j++)
{
if(vb[j]==(false))
      {
        x+=j;
      }
}

    Leider funktioniert die Addition nur bis zu einem Wert von ca. N=100 000,
    darüberhinaus scheint es fehlerhaft zu rechnen, z.B. bei 10^8 wird als Summe
    279209790387276 ausgegeben, richtig wäre aber 2556408900. Leider finde ich den Fehler nicht, hat jemand eine Idee, woran es liegen könnte?

    0
  • ?

    Ja. Ich habe das auch in dein Programm eingebaut. Ergebnis: Die erste Zahl 279209790387276 ist richtig. Das liefert auch das einfache Programm oben (SeppJ) und auch das von Camper (dessen Korrekturen leider nicht vorgestellt wurden).
    Folgerung: Du hast irgenwo an anderer Stelle noch eine Änderung vorgenommen und die vergessen oder übersehen.

    0
  • ?

    Ergänzung zum letzten Satz: bei der Kontrollrechnung.
    Vermutlich irgendwo int statt long.

    0
  • R

    Eratosthenes schrieb:

    Ja. Ich habe das auch in dein Programm eingebaut. Ergebnis: Die erste Zahl 279209790387276 ist richtig. Das liefert auch das einfache Programm oben (SeppJ) und auch das von Camper (dessen Korrekturen leider nicht vorgestellt wurden).
    Folgerung: Du hast irgenwo an anderer Stelle noch eine Änderung vorgenommen und die vergessen oder übersehen.

    Hm, also das Programm von SeppJ ( http://www.c-plusplus.net/forum/286485-10 ) als auch die Optimierung von Dir geben bei mir als Summe der Primzahlen bis 10^8 =
    2556408900 aus.
    Ich habe den Sourcecode mit Copy und Paste übernommen und kompiliert. (Codeblocks)

    0
  • R

    Erathostenes schrieb:

    Ergänzung zum letzten Satz: bei der Kontrollrechnung.
    Vermutlich irgendwo int statt long.

    Ok. Evtl. in Zeile 22 ( http://www.c-plusplus.net/forum/286485-20 )?

    Statt

    unsigned long anti_optimizer = 0;

    unsigned long long anti_optimizer = 0;
    0
  • C
    Mod

    Hier mal die überarbeitete Version meines Siebes.
    Das Programm akzeptiert 2 Argumente:
    1. die Grenze bis zu der (einschließlich) gesiebt werden soll
    2. eine Zahl p
    p bestimmt die Größe der Einzelsiebe ( Grenze/p ), als auch die Menge der Zahlen die überhaupt gesiebt werden müssen. Die optimale Wahl hängt deshalb von der Grenze ab - das Ergebnis sollte allerdings für jede beliebige Wahl korrekt sein. Im Allgemeinen sind solche p günstig für die die Eulersche φ-Funktion relativ klein ist, also etwa 2,6,30,210,2310,30030 usw.

    #include <vector>
#include <cmath>
#include <cassert>
#include <algorithm>
#include <iostream>
#include <ctime>
#include <boost/lexical_cast.hpp>
#include <boost/cstdint.hpp>

using std::cout;
using std::vector;
using std::sqrt;
using std::make_pair;
using std::pair;
using std::size_t;
using std::sort;
using boost::uint32_t;
using boost::uint64_t;
using boost::int64_t;
#ifdef M32
typedef uint32_t uintfast_t;
#else
typedef uint64_t uintfast_t;
#endif

typedef vector<char> sieve_type;

uint64_t sum = 0;
uint64_t prime_count = 0;
void add_prime(uint64_t v)
{
		++prime_count;
        sum += v;
}

uint32_t gcd(uint32_t a, uint32_t b)
{
	return b == 0 ? a : gcd( b, a % b );
}

pair<int32_t,int32_t> euclid(uint32_t a, uint32_t b)
{
    if ( b == 0 )
        return make_pair( 1, 0 );
    int32_t q = static_cast<int32_t>( a / b );
    pair<int32_t,int32_t> e = euclid( b, a % b );
	return make_pair( e.second, e.first - q * e.second );
}

struct sort_by_remainder : std::binary_function< pair<uint32_t,uint32_t>, pair<uint32_t,uint32_t>, bool >
{
	bool operator()(const pair<uint32_t,uint32_t>& lhs, const pair<uint32_t,uint32_t>& rhs) const
	{
		return lhs.second < rhs.second || ( lhs.second == rhs.second && lhs.first < rhs.first );
	}
};

vector<pair<uint32_t,uint32_t> > simple_sieve(uint32_t limit, uint32_t p)
{
	// save all coprimes <= limit with p
	assert( p != 0 );
    vector<pair<uint32_t,uint32_t> > result;
	uint32_t q = p;
	if ( limit >= 2 )
	{
		add_prime( 2 );
		if ( q % 2 != 0 )
            result.push_back( make_pair( 2 / p, 2 % p ) );
		else while ( q % 2 == 0 )
			q /= 2;
    }
    sieve_type sieve( limit );
    for ( uint32_t i = 3; i * i <= limit; i += 2 )
        if ( sieve[ i - 1 ] == 0 )
            for ( uint32_t j = i * i; j <= limit; j += 2 * i )
                sieve[ j - 1 ] = 1;
    for ( uint32_t i = 3; i <= limit; i += 2 )
        if ( sieve[ i - 1 ] == 0 )
		{
			add_prime( i );
			if ( q % i != 0 )
	            result.push_back( make_pair( i / p, i % p ) );
			else while ( q % i == 0 )
				q /= i;
        }
	if ( q != 1 )
		add_prime( q );
	sort( result.begin(), result.end(), sort_by_remainder() );
    return result;
}

vector<uint32_t> prime_remainders(uint32_t p)
{
    vector<uint32_t> result;
    for ( uint32_t i = 0; i < p; ++i )
		if ( gcd( p, i ) == 1 )
			result.push_back( i );
    return result;
}

vector<uint32_t> small_coeff(uint32_t p, const vector<uint32_t>& remainders)
{
	// find smallest k >= 0 for each remainder r with k * r = 1 mod p
	vector<uint32_t> result( p );
	for ( vector<uint32_t>::const_iterator it = remainders.begin(); it != remainders.end(); ++it )
		result[ *it ] = ( p + euclid( p, *it ).second ) % p;
	return result;
}

template<typename func>
void sieve_residue_class(uint64_t limit, uint32_t lower, uint32_t p, sieve_type& sieve, const vector<pair<uint32_t,uint32_t> >& small_primes, const vector<uint32_t>& coeff, func f, uintfast_t remainder)
{
	uintfast_t sieve_limit = static_cast<uintfast_t>( ( limit - remainder ) / p + 1 );
	std::fill_n( &sieve.front(), static_cast<size_t>( sieve_limit ), 0 );
	uintfast_t sec = 0, w = 0;
    for ( vector<pair<uint32_t,uint32_t> >::const_iterator it = small_primes.begin(); it != small_primes.end(); ++it )
    {
		if ( it->second != sec )
		{
			sec = it->second;
			w = static_cast<uint32_t>( static_cast<uint64_t>( remainder ) * coeff[ static_cast<size_t>( sec ) ] % p );
			if ( w < sec )
				w += p;
		}
		uintfast_t first = it->first;
		uintfast_t v = first * p + sec;
		uintfast_t x = static_cast<uintfast_t>( static_cast<uint64_t>( v ) * ( first * p + w ) / p );

        for ( ; x < sieve_limit; x += v )
            sieve[ static_cast<size_t>( x ) ] = 1;
    }
    for ( uintfast_t i = ( p + lower - remainder ) / p; i < sieve_limit; ++i )
        if ( sieve[ static_cast<size_t>( i ) ] == 0 )
            f( static_cast<uint64_t>( i ) * p + remainder );
}

template<typename func>
void sieve_by_remainders(uint64_t limit, uint32_t lower, uint32_t p, const vector<uint32_t>& remainders, const vector<pair<uint32_t,uint32_t> >& small_primes, const vector<uint32_t>& coeff, func f)
{
	uint32_t sieve_size = static_cast<uint32_t>( limit / p + 1 );
	sieve_type sieve( sieve_size );
	for ( vector<uint32_t>::const_iterator it = remainders.begin(); it != remainders.end(); ++it )
		sieve_residue_class( limit, lower, p, sieve, small_primes, coeff, f, *it );
}

int main(int argc, char** argv)
{
    std::time_t t = std::clock();
    uint64_t limit = argc > 1 ? boost::lexical_cast<uint64_t>(argv[1]) : 1000;
	if ( limit == 0 )
		limit = 1;
    uint32_t p = argc > 2 ? boost::lexical_cast<uint32_t>(argv[2]) : 30030;
	if ( p == 0 || p > limit )
		p = static_cast<uint32_t>( limit );
    uint32_t small_limit = static_cast<uint32_t>( sqrt( static_cast<double>( limit ) ) );
	vector<pair<uint32_t,uint32_t> > small_primes = simple_sieve( small_limit, p );
	vector<uint32_t> remainders = prime_remainders( p );
	vector<uint32_t> coeff = small_coeff( p, remainders );
	sieve_by_remainders( limit, small_limit, p, remainders, small_primes, coeff, add_prime );
    cout << "\nlimit: " << limit
		 << "\nsum: " << sum
		 << "\ncount: " << prime_count
		 << "\nc/l: " << (double)prime_count/(limit-1)
		 << "\nTime: " << (std::clock()-t)/(double)CLOCKS_PER_SEC
         << "\nsieve-size: " << limit / p + 1
         << "\np: " << p
         << "\nresidues: " << remainders.size() << '\n';
}

    Für 32-bit Targets ist es empfehlenswert, das Makor M32 zu setzen, damit bestimmte Berechnungen nur mit 32bit Integern ausgeführt werden.
    Damit erhalte ich folgende Zeiten (atom 330, gcc-4.5.2, 64bit):
    1e7 210 30ms
    1e8 2310 370 ms
    1e9 2310 4.34 s
    1e10 30030 43.7 s
    1e11 210210 900 s

    0
  • R

    Nicht schlecht, der Code, leider habe ich die Boostbibliotheken nicht installiert (Codeblocks), um ihn testen können. werde aber mal versuchen, die genannten Optimierungsvorschläge in mein Programm mit einzubauen.

    Danke soweit erstmal @all , habe viele Anregungen bekommen, was man noch besser machen kann. 🙂

    @Camper:
    Den Programmcode verstehe ich leider nur zu ca. 20% oder so 😕

    Hast du bei Einzelsieben die möglichen Algorithmen alle umgesetzt?,
    falls nicht, ließe sich da evtl. noch mehr optimieren.

    0
Anmelden zum Antworten