Mein Primzahlprogramm



  • Könnt ihr mal kurz ein Blick auf mein Programm werfen und schauen ob man am Vorgehen noch etwas verbessern könnte?
    Mir gehts um algorithmische Optimierungen und nicht darum mit Asm noch das letzte bischen rauszukitzeln.

    template< typename Listener >
    void findPrimes( unsigned to, Listener listener )
    {
    	if( !to ) return;
    
    	const size_t size = size_t( to / 2.0 + 0.5 );
    
    	vector< unsigned > primes;
    	primes.resize( size );
    
    	vector< bool > isPrime;
    	isPrime.resize( size );
    
    	primes[ 0 ] = 2;
    	isPrime[ 0 ] = true;
    
    	for( size_t i = 1, j = 3; i < size; ++i, j += 2 )
    	{
    		primes[ i ] = j;
    		isPrime[ i ] = true;
    	}
    
    	for( size_t i = 0; i < size; ++i )
    	{
    		for( size_t j = i + 1; j < size; ++j )
    		{
    			if( !isPrime[ j ] )
    				continue;
    			if( primes[ j ] % primes[ i ] == 0 )
    				isPrime[ j ] = false;
    		}
    	}
    
    	for( size_t i = 0; i < size; ++i )
    		if( isPrime[ i ] )
    			listener( primes[ i ] );
    }
    
    void printPrime( unsigned prime )
    {
    	cout << prime << " is a prime.\n";
    }
    
    int main()
    {
    	unsigned to;
    
    	while( ( cout << "Gebe das Limit ein: " ) && !( cin >> to ) )
    	{
    		cout << "\nSo nicht versuchs nochmal!" << endl;
    	}
    
    	findPrimes( to, printPrime );
    }
    


  • Naja, wenn du alle bisherigen Primzahlen kennst, reicht es, die Teilbarkeit der Zahl durch diese zu überprüfen.

    primes.push_back (2);
    
    for (unsigned j = 3; j <= to; j += 2)
    {
        std::vector<unsigned>::const_iterator i = primes.begin () + 1;
        for (; i != primes.end () && ((j % *i) != 0) ; ++i)
            ;
        if (i == primes.end ())
            primes.push_back (j);
    }
    

    Code ungetestet. Das bekommt man sicherlich auch noch irgendwie schöner mit der STL hin.



  • Die Idee ist gut, hab selber schon überlegt wie ich den Speicherverbrauch reduzieren kann und so vermeide ich gleich noch unnötige Tests.



  • Hab sie gleich mal umgesetzt:

    template< typename Listener >
    void findPrimes( unsigned to, Listener listener )
    {
    	if( to < 2 ) return;
    
    	vector< unsigned > primes;
    
    	primes.push_back( 2 );
    
    	for( size_t i = 3; i <= to; i += 2 )
    	{
    		vector< unsigned >::const_iterator it = primes.begin(), end = primes.end();
    
    		while( i % *it != 0 )
    		{
    			++it;
    			if( it == end )
    			{
    				primes.push_back( i );
    				break;
    			}
    		}
    	}
    
    	for_each( primes.begin(), primes.end(), listener );
    }
    
    void printPrime( unsigned prime )
    {
    	cout << prime << " is a prime.\n";
    }
    
    int main()
    {
    	unsigned to;
    
    	while( ( cout << "Gebe das Limit ein: " ) && !( cin >> to ) )
    	{
    		cout << "\nSo nicht versuchs nochmal!" << endl;
    	}
    
    	findPrimes( to, printPrime );
    }
    


  • Mir fällt gerade ein, dass ich den listener in dem if bereits benachrichtigen kann.



  • hab den code nur überflogen. Aber weist du dass eine zahl niemals einen Teiler haben kann der größer als sqrt(zahl) ist ? Ich konnte damit meinen Primzahl algorithmus um paar 100% beschleunigen.



  • walljumper schrieb:

    hab den code nur überflogen. Aber weist du dass eine zahl niemals einen Teiler haben kann der größer als sqrt(zahl) ist ? Ich konnte damit meinen Primzahl algorithmus um paar 100% beschleunigen.

    Weiß ich habs vorhin auch mal implementiert, nur ging es dann nicht, allerdings hab ich ganz vergessen, dass ich mein Programm geändert hatte es gibt gar nicht mehr die Primzahlen aus sondern wieviele es sind. Gut dass du mich nochmal dran erinnert hast 🙂



  • Jetzt ist es richtig schnell, er braucht für die Primzahlen bis 1mio nur noch ~10Sekunden.
    Mit der zweiten Version oben brauchte er ungefähr so lange für die Primzahlen bis 100.000 (welcher er jetzt quasi sofort hat).



  • hmm muss an dem vielen stl liegen oder den templates weil mein code macht eigentlich das selbe nur das ich bis 1mio nur ~0,4 sec brauch.

    #include <iostream>
    #include <math>
    #include <time>
    #pragma hdrstop
    using namespace std;
    //---------------------------------------------------------------------------
    
    #pragma argsused
    int main()
    {
     int Primes[78499];
     int Zahl=1;
     int stop;
     int Prim=0;
     bool isPrime=true;
     double time1,time2;
    
     time1= -clock();
     Primes[0]=2;
     while(Zahl<=1000000)
    	{
    	Zahl++;
    	isPrime=true;
    	stop = sqrt(Zahl);
    	 for(int i=0;Primes[i]<=stop;i++)
         	{
    			if(Zahl%Primes[i]==0)
                	{
    				isPrime=false;
    				break;
                    }
            }
    	 if(isPrime)
         	{
    		Prim++;
    		Primes[Prim]=Zahl;
    
            }
    	}
    	time2=clock();
    	cout<<Prim+1<<endl<<(time1+time2)/CLOCKS_PER_SEC;
    	getchar();
    
    	return 0;
    }
    


  • Jetzt wo du es sagst, ich hab ja noch das Debug-Build aktiviert 🙄



  • Gibts eine Möglichkeit abzuschätzen wieviele Primzahlen bis zur Obergrenze vorhanden sind?

    Mir geht es nicht darum es genau abzuschätzen, ich weiß, dass das nicht möglich ist, aber gibt es einen groben Richtwert, der in etwa hinkommt, so dass ich bereits am Anfang das Array näherungsweise auf die tatsächliche Größe bringen kann?

    Ich glaub so kann ich auch nochmal etwas Zeit gewinnen.



  • Für die Primzahlen bis 100mio brauchte mein Programm jetzt ~3Minuten 45Sekunden.



  • Anzahl von Primzahlen für bestimmte Grenzen:
    http://www.michael-holzapfel.de/themen/primzahlen/pz-anzahl.htm

    Hatte da auch mal was gemacht,w erde gleich mal gucken, wo ichs hab.

    Edit:
    Habs gefunden. Das braucht bis 10 Mio 470 ms. Basiert auf dem Sieb des Eratosthenes.

    #include <sstream>
    using std::istringstream;
    #include <iostream>
    using std::cout;
    using std::endl;
    #include <cmath>
    #include <ctime>
    using std::clock;
    //using std::CLOCKS_PER_SEC;
    
    long long countPrime(long long);
    
    int main(int argc, char *argv[])
    {
        if(argc < 2)
            return 2;
        istringstream isst;
        long long border, res;
        isst.str(argv[1]);
        isst >> border;
        long start = clock();
        res = countPrime(border);
        long end = clock();
        cout << res << " Primzahlen bis " << border << " in " << ((end-start)/static_cast<double>(CLOCKS_PER_SEC))*1000 << " ms bei " << (sizeof(char)*border+1)/(1024.0f*1024.0f) << " MB ";
    
    	return 0;
    }
    
    long long countPrime(long long border)
    {
        char *mat = new char[border+1];
        mat[0] = mat[1] = 0;
        for(long long i=2;i<border+1;i++)
        {
            mat[i] = 0;
            mat[++i] = 1;
        }
    
        long long s = static_cast<long long> (std::sqrt(border));
        long long p = 3;
        long long primeCount = border/2;
        bool foundNewP = false;
    
        while(p <= s)
        {
            for(long long i=p*3;i<=border;i+=p*2)
            {
                if(mat[i] != 0)
                {
                    primeCount--;
                    mat[i] = 0;
                }
            }
    
            foundNewP = false;
            for(long long i=p+1;i<=s;i++)
            {
                if(mat[i] != 0)
                {
                    p = i;
                    foundNewP = true;
                    break;
                }
            }
    
            if(!foundNewP)
                break;
        }
    
        delete[] mat;
    
        return primeCount;
    }
    

    Fehlen nur leider jegliche Kommentare, weil ich da einige Versuche gemacht hab und keine Lust immer alle zu kommentieren :>



  • hab den alg noch bissel verbessert die Änderungen sind markiert brauche jetzt nur noch ~0,18 sec

    //---------------------------------------------------------------------------
    #include <iostream>
    #include <math>
    #include <time>
    #pragma hdrstop
    using namespace std;
    //---------------------------------------------------------------------------
    const int lim=1000000;
    #pragma argsused
    int main()
    {
     double time1,time2;
     time1= -clock();
     unsigned int *Primes= new unsigned int[lim];
     bool *x= new bool[lim];   //   Änderung
     for(int i=0;i<lim;i++)    //   Änderung
    	x[i]=false;            //   Änderung
    
     unsigned int end=0;  //   Änderung
     unsigned int Zahl=1;
     unsigned int stop;
     unsigned int Prim=0;
     bool isPrime=true;
    
     Primes[0]=2;
     while(Zahl<=lim)
    	{
    	Zahl+=2;
    	if(x[Zahl]) continue;
    	isPrime=true;
    	stop = sqrt(Zahl);
    	 for(int i=0;Primes[i]<=stop;i++)
    		{
    			if(Zahl%Primes[i]==0)
    				{
    				isPrime=false;
    				break;
    				}
    		}
    	 if(isPrime)
    		{
    		Prim++;
    		Primes[Prim]=Zahl;
    		end=lim/Zahl;           	 	//   Änderung
    		for(unsigned int i=2;i<end;i++) //   Änderung
    			x[i*Zahl]=true;      		//   Änderung
    
    		}
    	}
    	time2=clock();
    	cout<<Prim+1<<endl<<(time1+time2)/CLOCKS_PER_SEC;
    	getchar();
    
    	return 0;
    }
    

    bis 100 mio brauch ich 126 sec
    bis 10 mio 4,5 sec



  • Ich frag mich warum eure beiden so viel schneller sind wie meiner.
    Testet meinen mal bitte bei euch und sagt mal wie lang er da braucht.

    unsigned calcNext( unsigned n, unsigned i )
    {
    	return ( n + i / n ) >> 1;
    }
    
    unsigned int isqrt( unsigned number ) 
    {
    	unsigned int n  = 1;
    	unsigned int n1 = calcNext( n, number );
    
    	while( abs( int( n1 - n ) ) > 1 ) 
    	{
    		n  = n1;
    		n1 = calcNext( n, number );
    	}
    	while( ( n1 * n1 ) > number )
    	{
    		n1 -= 1;
    	}
    	return n1;
    }
    
    template< typename Listener >
    void findPrimes( unsigned to, Listener listener )
    {
    	if( to < 2 ) return;
    
    	vector< unsigned > primes;
    	primes.reserve( to / 2 );
    
    	primes.push_back( 2 );
    	listener( 2 );
    
    	for( unsigned i = 3; i <= to; i += 2 )
    	{
    		vector< unsigned >::const_iterator it = primes.begin(), end = primes.end();
    
    		const unsigned sqr = unsigned( isqrt( double( i ) ) );
    
    		while( i % *it != 0 )
    		{
    			++it;
    			if( *it > sqr || it == end )
    			{
    				primes.push_back( i );
    				listener( i );
    				break;
    			}
    		}
    	}
    }
    
    unsigned printPrime( unsigned prime )
    {
    	static unsigned size = 0;
    
    	return ++size;
    }
    
    __int64 GetPerformanceFrequency()
    {
    	__int64 frequ = 0;
    	QueryPerformanceFrequency( reinterpret_cast< LARGE_INTEGER* >( &frequ ) );
    	return frequ;
    
    	}
    
    	__int64 GetTimeMS()
    	{
    		static __int64 frequ = GetPerformanceFrequency() / 1000;
    		__int64 perfcount;
    		QueryPerformanceCounter( reinterpret_cast< LARGE_INTEGER* >( &perfcount ) );
    		return perfcount / frequ;
    	}
    
    int main()
    {
    	unsigned to;
    
    	while( ( cout << "Gebe das Limit ein: " ) && !( cin >> to ) )
    	{
    		cout << "\nSo nicht versuchs nochmal!" << endl;
    	}
    
    	__int64 start = GetTimeMS();
    
    	findPrimes( to, printPrime );
    
    	__int64 end = GetTimeMS();
    
    	cout << printPrime( 1 ) - 1 << endl << "Duration: " << end - start << "ms";
    
    }
    


  • Sehe gerade, ich hab noch nen alten Überrest im Code gehabt:

    wurzel sollte so aussehen:
    
    const unsigned sqr = isqrt( i );
    


  • dein code braucht bei mir bis 100mio 4 minuten und 24 sec

    dein code ist langsamer weil du die stl verwendest.Wenn ich vectoren verwenden würde wäre mein code auch langsamer.

    viande hat den bis jetzt wohl besten algorithmus wobei ich da noch etwas optimierungspotential in der implemetierung sehe.
    z.B. im Schleifen kopf so wenig wie möglich rechnen die abbruch bedingung vorher ausrechnen da sie sonst bei jedem durchlauf erneut berechnet wird.



  • Mh nagut, kann mir aber vorstellen, dass der Compiler sowas selber macht. Hatte das irgendwie übersehen 😉
    Achso, bis 100 Mio brauch das Ding 5 Sek und 1 Milliarde passt leider nicht in meinen Ram (hab nur 768 MB :(). Wollte das Ding auch mal so umbauen, dass man eine Zahl pro Bit hat und nicht pro byte, aber war irgendwie nicht dazu gekommen. Vielleicht setz ich mich da nochmal ran 🙂



  • Der erste Ansatz, funktioniert zwar und verbraucht wesentlich weniger Speicher, aber ist auch ziemlich lahm und hässlich 😞 Braucht etwa 20 Sek. auf meinem gammel Rechner hier auf der Arbeit bis 10 Mio.
    Zum Vergleich: der andere braucht 1,6 Sek.

    PS: Wer Angst vor Augenkrebs hat, sollte sich das nicht angucken 😉
    PPS: RDTSC holt CPU-Ticks. Ist eigentlich von INTEL, aber AMD hat das inzwischen auch. Ist aber für den Algo irrelevant.

    #include <sstream>
    using std::istringstream;
    #include <iostream>
    using std::cout;
    using std::endl;
    #include <cmath>
    #include <ctime>
    using std::clock;
    
    #define RDTSC(llptr) { \
    __asm__ __volatile__ ( \
    "rdtsc" \
    : "=A" (llptr) \
    ); }
    
    long long countPrime(long long);
    
    namespace bit
    {
        char const bit[] = {
        0x80, //10000000
        0x40, //01000000
        0x20, //00100000
        0x10, //00010000
        0x8,  //00001000
        0x4,  //00000100
        0x2,  //00000010
        0x1   //00000001
        };
        const char all = 0xFF;//11111111
        const char odd = 0x55; //01010101 -> for example the second byte:
                                 //8, 9, 10, 11, 12, 13, 14, 15 we set only the odd numbers
    }
    
    namespace revBit
    {
        char const bit[] = {
        0x7F, //01111111
        0xBF, //10111111
        0xDF, //11011111
        0xEF, //11101111
        0xF7, //11110111
        0xFB, //11111011
        0xFD, //11111101
        0xFE  //11111110
        };
        const char all = 0x0; //00000000
    }
    
    int main(int argc, char *argv[])
    {
        if(argc < 2)
            return 1;
        istringstream isst;
        long long border, res;
        isst.str(argv[1]);
        isst >> border;
    
        long start = clock();
        res = countPrime(border);
        long end = clock();
        cout << res << " Primzahlen bis " << border << " in " << ((end-start)/static_cast<double>(CLOCKS_PER_SEC))*1000 << " ms at " << (sizeof(char)*std::floor(border/8.0f))/(1024.0f*1024.0f) << " MB ";
    
    	return 0;
    }
    
    long long countPrime(long long border)
    {
        //unsigned long long start, end;
        //RDTSC(start);
    
        long long size = static_cast<long long>(std::ceil(border/8.0f));
        char *mat = new char[size];
        for(long long i=0;i<size;i++)
            mat[i] = bit::odd;
        mat[0] &= revBit::bit[1]; //unset 1
        for(int i=8-(size*8-border)+1;i<8;i++)//unset the unused bits at the end
            mat[size-1] &= revBit::bit[i];
    
        long double s = std::sqrt(border);
        long long sSize = static_cast<long long>(std::ceil(s/8));
        long long p = 3;
        long long sVal = p * 3;
        long long addVal = p * 2;
        long long primes = border/2; //we exclude the even numbers
        bool foundNewP = false;
        long long posArr = 0, posBit = 4; //position of 3 + 1 (next element)
    
        //RDTSC(end);
        //cout << "clock cycles (init): " << end - start << endl;
    
        while(p < s)
        {
            for(long long i=0;i<size;i++)
            {
                //RDTSC(start);
                if(sVal > 7)
                {
                    sVal = sVal - 8;
                    continue;
                }
                //RDTSC(end);
                //cout << "clock cycles (eleminate): " << end - start << endl;
                for(;sVal<8;sVal+=addVal)
                {
                    if(mat[i] & bit::bit[sVal])
                    {
                        mat[i] &= revBit::bit[sVal];
                        primes--;
                    }
                }
                if(sVal > 7)
                    sVal = sVal - 8;
            }
    
            //RDTSC(start);
            foundNewP = false;
            for(;posArr<sSize;posArr++)
            {
                for(;posBit<8;posBit++)
                {
                    if(mat[posArr] & bit::bit[posBit])
                    {
                        p = posArr * 8 + posBit;
                        sVal = p * 3;
                        addVal = p * 2;
                        foundNewP = true;
    
                        if(posBit >= 7)
                        {
                            posArr++;
                            posBit = 0;
                        }
                        else
                            posBit++;
                        break;
                    }
                }
                if(foundNewP)
                    break;
                if(posBit >= 7)
                    posBit = 0;
            }
            //RDTSC(end);
            //cout << "clock cycles (search): " << end - start << endl;
            if(!foundNewP)
                break;
        }
    
        delete[] mat;
    
        return primes;
    }
    


  • Ich weiß, dass Doppelposting nervt, aber deswegen habe ich jetzt mal eingeloggt :p
    Langsam fängt das an Spaß zu machen. Also, ich habe das mit den einzelnen Bits nochmal vernünftig gemacht. Das läuft etwa doppelt so schnell wie mein Ausgangsversuch und verbraucht dabei auch noch viel weniger Speicher. Ich habe auch mal ein paar Kommentare reingeworfen.

    #include <climits>
    #include <cmath>
    #include <ctime>
    
    #include <sstream>
    #include <iostream>
    
    using std::istringstream;
    using std::cout;
    using std::endl;
    using std::clock;
    
    //Bits pro Wort
    #define BITS_PER_WORD (CHAR_BIT*sizeof(unsigned long))
    //Anzahl der benötigten Wörter, um die Bits aufzunehmen
    #define WORDS(n) \
     ((n) < 1 ? 1 : ((n) + BITS_PER_WORD - 1)/BITS_PER_WORD)
    
    typedef unsigned long word;
    typedef unsigned long long ull;
    
    static const word WORD_MAX = ULONG_MAX;
    
    ull countPrime(ull);
    
    namespace bit
    {
        word bit[BITS_PER_WORD]; //Bits an einer bestimmten Position
        word revBit[BITS_PER_WORD]; //alle Bits, bis auf das einer bestimmten Position
                                    //(Umkehrung von bit::bit)
        const word odd = 0xAAAAAAAA; //setzt alle ungeraden Bits für ein 32 Bit Int
    
        void init() //initialisierung der Bits(keine Lust die ganzen Matrizen von Handzu machen :p)
                    //ich hoffe einfach mal, dass der Compiler sieht, dass das alles Konstanten sind :o
        {
            for(unsigned int i=0;i<BITS_PER_WORD;i++)
                bit[i] = static_cast<word>(1) << (i % BITS_PER_WORD);
            for(unsigned int i=0;i<BITS_PER_WORD;i++)
                revBit[i] = ~(static_cast<word>(1) << (i % BITS_PER_WORD));
        }
    
        bool test_bit(word* mat, ull _pos) //Bit in _pos testen
        {
            return mat[_pos/BITS_PER_WORD] & bit[_pos % BITS_PER_WORD];
        }
    
        void unset_bit(word* mat, ull _pos) //Bit in _pos auf 0 setzen
        {
            mat[_pos/BITS_PER_WORD] &= revBit[_pos % BITS_PER_WORD];
        }
    
        void set_odd_bits(word* mat, ull _size) //alle ungeraden Bits setzen
        {
            for(unsigned long long i=0;i<_size;i++)
            {
                mat[i] = odd;
            }
        }
    }
    
    int main(int argc, char *argv[])
    {
        if(argc < 2)
            return 1;
        istringstream isst;
        ull border, res;
        isst.str(argv[1]);
        isst >> border;
    
        bit::init();
    
        long start = clock();
        res = countPrime(border);
        long end = clock();
    
        cout << res << " Primzahlen bis " << border
        << " in " << ((end-start)/static_cast<double>(CLOCKS_PER_SEC))*1000
        << " ms bei " << (sizeof(word)*WORDS(border))/(1024.0f*1024.0f) << " MB ";
    
        return 0;
    }
    
    ull countPrime(ull border)
    {
        size_t size = WORDS(border+1);
        word *mat = new word[size];
    
        bit::set_odd_bits(mat, size);
    
        //0 und 1 sind keine Primzahlen :)
        bit::unset_bit(mat, 0);
        bit::unset_bit(mat, 1);
    
        ull s = static_cast<unsigned long long> (std::sqrt(border)); //obere Schranke
    
        ull p = 3, t; //da die 2 bereits ausgeschlossen ist, fangen wir bei 3 an
    
        ull primeCount = border/2; //die 2 ist weg, also auch die Hälfte der möglichen Primzahlen
    
        bool foundNewP = false; //wird gesetzt, wenn eine neue Primzahl gefunden wird
    
        while(p <= s) //es wird nur bis zur oberen Schranke iteriert
        {
            t = p * 2; //die Addition von 2 ungeraden Zahlen ergibt wiederum eine gerade
                        //die brauchen wir aber nicht testen
    
            //@i = p * 3 : p selber ist ist prim, wird daher nicht ausgeschlossen
            //P * 2 ist gerade, somit fangen wir erst bei p * 3 an
            for(ull i=p*3;i<=border;i+=t)
            {
                if(bit::test_bit(mat, i))
                {
                    primeCount--;
                    bit::unset_bit(mat, i);
                }
            }
    
            foundNewP = false;
            for(unsigned long long i=p+1;i<=s;i++)
            {
                if(bit::test_bit(mat, i))
                {
                    p = i;
                    foundNewP = true;
                    break;
                }
            }
    
            if(!foundNewP)
                break;
        }
    
        delete[] mat;
    
        return primeCount;
    }
    

Anmelden zum Antworten