Mein Primzahlprogramm

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

wurzel sollte so aussehen:

const unsigned sqr = isqrt( i );

walljumper

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

viande

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

Ich hab im Netz mal einen Prime-Sieve basierenden Code gefunden:

Zeiten auf meinem Arbeitsrechner, AMD Athlon XP 3200+ (ist noch ein alter 32bitter), 1 GB RAM:
1 - 1.000.000: 0,219 Sek
1 - 10.000.000: 0,296 Sek
1 - 100.000.000: 1,156 Sek
1 - 1.000.000.000: 9,875 Sek

Der Code ist von Kim Walisch, ich finde im Moment den Link aber nicht mehr. Was schnelleres habe ich damals nicht gefunden...

Auf planet3dnow hat einer bis 4mrd in 1,5sec geschafft

Aber bestimmt nicht auf so einer lahmen Kröte wie meiner...

Ich hab den Link wiedergefunden: http://www.primzahlen.de/files/referent/kw/index.htm

walljumper

ne war ein Core2duo @3,2ghz trotzdem ist sein alg um einiges besser

Was ich erstaunlich find ist das mein code bis 1mio wesentlich schneller ist als deiner aber bei 10mio ist meiner viel langsamer.

hab jetzt bis 1mio 0,032sec
bis 10 mio ~2,2 sec

mein code wird nur langsam etwas unübersichtlich ^^ sollte da vieleicht mal aufräumen.

//---------------------------------------------------------------------------
#include<iostream>
#include<time.h>
#include<fastmath.h>
#pragma hdrstop
using namespace std;
//---------------------------------------------------------------------------
typedef unsigned int Int;
const Int lim=1000000;

int main()
{
double time1,time2,time3,time4;
time1=-clock();
// TEST ALGORITHMUS BEGIN
{
 bool *test = new bool[lim];
 for(Int i=0;i<lim;++i)
	test[i]=false;
 Int *Primes = new Int[lim];	

 Int firstborder = sqrt(lim);
 Int end;
 Int Zahl=1;
 Int stop;
 Int Prim=0;
 Int P2;
 Int tmp;
 Int start;
 bool isPrime;

 Primes[0]=2;

 while(Zahl<=firstborder)
	{
	Zahl+=2;
	if(test[Zahl]) continue;
	isPrime=true;
	stop = sqrt((float)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;
		for(Int i=2;i<end;++i)
			test[i*Zahl]=true;

		}
	}
	for(Int i=firstborder;i<lim;++i)
		{
			test[i]=true;
		}
 P2=Prim;
 time3=-clock();
 for(Int i=0;i<Prim;++i)
	{
		end=lim/Primes[i];
		start=firstborder/Primes[i];
		for(Int n=start;n<=end;++n)
			{
			tmp=n*Primes[i];
			test[tmp]=false;
			}
	}
	time4=clock();
 for(Int i=firstborder;i<lim;++i)
		{
			if(test[i]) ++P2;
		}
cout<<P2;
}
// TEST ALGORITHMUS END
time2 = clock();
cout<<endl<<(time1+time2)/CLOCKS_PER_SEC<<endl<<(time3+time4)/CLOCKS_PER_SEC;
getchar();
	return 0;
}

Gregor

Hmmm... in einem Primzahlenthread darf mein altes Primzahlzählprogramm nicht fehlen. ...in Java:

public class TestPrime
{
   public TestPrime ()
   {
   }

   public static void main (String[] args)
   {
      int end = 1000000000;
      long time = System.currentTimeMillis ();
      int [] sieve = new int [(end >> 5) + 1];
      int x, y, i;
      int primes = 1;
      int sqrt = (int)Math.sqrt((double)end);
      x = 3;
      while (x <= sqrt)
      {
         if ((sieve[x >> 5] & (0x1 << (x & 0x1f))) == 0)
         {
            y = x * x;
            i = x << 1;
            ++primes;
            while (end > y)
            {
               sieve[y >> 5] |= (0x1 << (y & 0x1f));
               y += i;
            }
         }
         x += 2;
      }
      while (x <= end)
      {
         if ((sieve[x >> 5] & (0x1 << (x & 0x1f))) == 0) ++primes;
         x += 2;
      }
      System.out.println ("Zwischen 0 und " + end +
                          " liegen " + primes + " Primzahlen.");
      System.out.println ("GesamtZeit : " +
                          (System.currentTimeMillis () - time) +
                           " Millisekunden");
   }
}

gregor@linux:~/JavaProjects/Test/TestPrime> /opt/jdk1.6.0/bin/java -Xmx192m TestPrime
Zwischen 0 und 1000000000 liegen 50847534 Primzahlen.
GesamtZeit : 41291 Millisekunden

Auf einem Pentium M mit 1,86GHz, 512MB RAM.

BTW: Bei den neuen Mehrkernprozessoren könnte man sich ja mal überlegen, ob man da irgendwie einen parallel arbeitenden Algorithmus schreiben kann.

EDIT: volkard hatte damals noch etwas deutlich schnelleres anzubieten. Vielleicht postet er seinen Code ja auch nochmal.