Sehr schneller Primzahl-Finder

[GTA-Place]

Das DF (besonders Phantom1) hat mal weiter probiert und Phantom1 hat diese Funktion geschrieben:

Delphi-Quellcode:

procedure SavePrimes(MaxPrime: Cardinal; const FileName: string = '');  
const  
  CACHE = 64*1024;  
  STEMPEL: array[0..7] of Byte = (1, 7, 11, 13, 17, 19, 23, 29);  
  MODS: array[0..29] of Byte = (0, 1, 0, 0, 0, 0, 0, 2, 0, 0, 0, 4, 0, 8, 0, 0,  
                                0, 16, 0, 32, 0, 0, 0, 64, 0, 0, 0, 0, 0, 128);  
var  
  Primes, PrimesLUT: array of Byte;  
  i, j, k, PrimeLen, PrimeBits, Num, Num2, m, mbit, s: Cardinal;  
  f: TextFile;  
begin  
  if FileName<>'' then begin  
    AssignFile(f, FileName);  
    ReWrite(f);  
    WriteLn(f, '2'+#10#13+'3'+#10#13+'5');  
  end;  
 
  SetLength(PrimesLUT, Trunc(Sqrt(MaxPrime)/30)); // max 2184 Byte für 2^32 ;-)  
  PrimesLUT[0]:=$01;  
  PrimeLen:=Length(PrimesLUT);  
  PrimeBits:=PrimeLen*30;  
  for i:=0 to Trunc(Sqrt(PrimeBits)/30) do  
    for j:=0 to 7 do  
      if PrimesLUT[i] and (1 shl j)=0 then begin  
        s:=STEMPEL[j];  
        Num:=i*30+s;  
        Num2:=Num*Num;  
        mbit:=Num2 mod 30;  
        m:=(Num2-mbit) div 30;  
        while m<PrimeLen do begin  
          PrimesLUT[m]:=PrimesLUT[m] or MODS[mbit];  
          Inc(m, i);  
          Inc(mbit, s);  
          if mbit>29 then begin  
            Dec(mbit, 30);  
            Inc(m);  
          end;  
        end;  
      end;  

  SetLength(Primes, CACHE);  
  PrimeLen:=Length(Primes);  
  PrimeBits:=PrimeLen*30;  
  for k:=0 to MaxPrime div PrimeBits do begin  
    FillChar(Primes[0], PrimeLen, 0);  
    for i:=0 to Trunc(Sqrt((k+1)*PrimeBits)/30) do  
      for j:=0 to 7 do  
        if PrimesLUT[i] and (1 shl j)=0 then begin  
          s:=STEMPEL[j];  
          Num:=i*30+s;  
          if k=0 then  
            Num2:=Num*Num  
          else  
            Num2:=Trunc(k*PrimeBits/Num)*Num+Num;  
          mbit:=Num2 mod 30;  
          m:=(Num2-mbit) div 30-k*PrimeLen;  
          while m<PrimeLen do begin  
            primes[m]:=Primes[m] or MODS[mbit];  
            Inc(m, i);  
            Inc(mbit, s);  
            if mbit>29 then begin  
              Dec(mbit, 30);  
              Inc(m);  
            end;  
          end;  
        end;  
    if FileName<>'' then  
      for i:=0 To PrimeLen-1 do  
        for j:=0 to 7 do begin  
          if k*PrimeBits+i*30+STEMPEL[j]>MaxPrime then  
            Break;  
          if not ((i=0) and (j=0) and (k=0)) and (Primes[i] and (1 shl j)=0) then  
            WriteLn(f, IntToStr(k*PrimeBits+i*30+STEMPEL[j]));  
        end;  
  end;  
  if FileName<>'' then  
    CloseFile(f);  
end;

Um 500.000.000 Zahlen zu überprüfen, braucht diese Funktion nur ca. 4 Sekunden.


EDIT: Code wurde von Phantom1 verbessert: 500 mio Zahlen in ca. 2.65 Sekunden.

[negaH]

Toll, das ist in fact ein 8/30 Comb Sieve, so wie es in meinem Posting im Attachment zu finden ist, sehr wahrscheinlich sogar von meinen Source abgekupfert.

Zitat:

Delphi-Quellcode:

 SetLength(PrimesLUT, Trunc(Sqrt(MaxPrime)/30)); // max 2184 Byte für 2^32 ;-)

Kannst du mir mathematisch erklären werum im Remark "max 2184 Bytes für 2^32" steht ? bzw. wie sich diese Schranke auf Grund des Verfahrens ergibt ?


Der Unterschied ist halt nur in der Performance und Flexibilität zu suchen.

Zitat:

EDIT: Code wurde von Phantom1 verbessert: 500 mio Zahlen in ca. 2.65 Sekunden.

Auf welchem Rechner ermittelt ?

Mein obiges Sieb schafft die 500Mio in 2.7 Sekunden auf einem P4 1.5GHz.
Falls die 2.65 Sekunden auch auf einem P4 1.5GHz ermittelt wurden so wundert mich dies schon, da mein Sieb mit hochoptimierten Assembler Routinen arbeitet die im Grunde wesentlich schneller als reiner PASCAL Source sein sollten.

Gruß Hagen

[Phantom1]

Zitat:

Zitat von negaH

Toll, das ist in fact ein 8/30 Comb Sieve, so wie es in meinem Posting im Attachment zu finden ist, sehr wahrscheinlich sogar von meinen Source abgekupfert.

Ich hab mir dein Code bis gerade eben noch nicht angeschaut. Sieht ziemlich komplex aus, schade das ich kein Assembler kann.

Zitat:

Zitat von negaH

Kannst du mir mathematisch erklären werum im Remark "max 2184 Bytes für 2^32" steht ? bzw. wie sich diese Schranke auf Grund des Verfahrens ergibt ?

Mein Programm sucht nur primzahlen bis 2^32. Da man nur die Primzahlen bis zur Wurzel(2^32) braucht um die restlichen wegzustreichen zu können, wären das 65536. Da ich ein 30er Stempel verwende ->65536/30 = 2184 Bytes nötig um alle Primzahlen bis 65536 (als bits) zu speichern. Dieses Array (PrimeLUT) benutze ich dann für die spätere Berechnung der einzelnen Array-Abschnitte und spare mir dadurch viele unnötige Streichungen.

Zitat:

Zitat von negaH

Auf welchem Rechner ermittelt ?
Mein obiges Sieb schafft die 500Mio in 2.7 Sekunden auf einem P4 1.5GHz.
Falls die 2.65 Sekunden auch auf einem P4 1.5GHz ermittelt wurden so wundert mich dies schon, da mein Sieb mit hochoptimierten Assembler Routinen arbeitet die im Grunde wesentlich schneller als reiner PASCAL Source sein sollten.

Mein Rechner ist ein Athlon XP-M@2000MHz. Gemessen wurde mit RealTime-Priority.

mfg
Phantom1

[PierreB]

Zitat:

Gemessen wurde mit RealTime-Priority.

Aähm nur mal kurz ne ganz kleine Frage, ich steh grad aufm Schlauch: Wie kann ich denn die Zeit die ein Algorithmus braucht messen ?? :gruebel:

[negaH]

Zitat:

Ich hab mir dein Code bis gerade eben noch nicht angeschaut. Sieht ziemlich komplex aus, schade das ich kein Assembler kann.

Alle relevanten Assemblerparts sind auch auskommentiert in reinem PASCAL in diesem Source enthalten.
Verzichtet man in meiner Implementation auf alle zusätzlichen Features, wie:
- Index Berechnungen von Primzahlen
- speichern, laden, erzeugen einer Cache Datei
- Anpassungen des Sources für D3 etc.
- direkte Berechnung der Residue Tabelle zu einem beliebigen Startwert

und packt dies alles in eine einzigste Funktion, so dürfte ein Source rauskommen der zu deinem Source fast 1 zu 1 identisch ist.

Es verwundert mich eben schon weil ich damals definitiv im Netz keinen einzigsten Source finden konnte der dieses 8/30Comb Sieb implemntierte. Meine Referenz war ein PostScript von D.J.Bernstein, eine mathematiche Abhandlung über dieses Sieb.

Gruß Hagen

[Kara]

<OT>Wisst ihr, was mir zu diesem Thread einfällt?
"Meiner ist läää-schneller" :stupid: </OT>

[negaH]

@Kara:

die Schnelligkeit einer Implementierung hängt von verschiedenen Faktoren ab:

1.) der mathematische Algorithmus bestimmt die grundsätzliche Laufzeitkomplexität der Implementation
2.) die programmiertechnische Umsetzung der Datenstrukturen und Algorithmen bestimmt die Gesamtlaufzeit

Man kann sich immer hinstellen und die Arbeit anderer diskreditieren ohne es real besser machen zu können oder eben nicht im entferntesten verstehen zu können welche Arbeit hinter einer solchen Umsetzung tatsächlich steckt.

Der schlußendlich Vergleich der Effizienz einer jeweiligen Implementierung hat also rein garnichts "Penis-Längen" zu tuen, sondern mit dem Ehrgeiz eines Programmierers seine Arbeit so gut wie nur möglich zu machen. Das sich daraus direkt auch der Stolz auf die eigene Arbeit ergibt, die man natürlicherweise auch mit der Arbeit anderer vergleichen will, ist doch wohl ein rein menchliches Bedürfnis.

Wenn man nun bedenkt das die Primzahlfindung als Operation wiederum die Laufzeiten der darauf basierenden Algorithmen beeinflusst so ist es ein wichtiges Bedürfnis diese Primzahlberechnungen so effizient wie nur möglich hinzubekommen.

Sich also, wenn auch offtopic, sich über solche Vergleiche lustig zu machen, deutet nur daraufhin das man das grundsätzliche Problem überhaupt nicht verstehen konnte. (schade eigentlich da du anscheinend nicht bereit bist an deinen Schwächen zu arbeiten).

Denn ich behaupte mal, das jeder Poster in diesem Thread, mit seiner eigenen Implementierung sich viel Mühe und Arbeit gemacht hat sie
1.) mathematisch zu begreifen
2.) sich im WEB nach relevanten Informationen auf Suche zu begeben
3.) lange an der Implementierungen rumgefeilt hat um sie so effizient wie möglich hinzubekommen

Es steckt also in jedem Posting hier viel Zeit drinnen.
Statt davon zu lernen kann man sich auch darüber lustig machen, bravo ;)

EDIT:
4.) sie dann hier in der DP der Allgmeinheit zur Verfügung stellt damit Leute wie DU daraus lernen können.

Gruß Hagen

[Kara]

:roll:
Ist ja gut. Du brauchst dich durch meinen kleinen Einwurf nicht angegriffen fühlen.

[negaH]

Warum ist dann überhaupt ein solcher Einwurf nötig gewesen ?

Einfach nur um was zu posten, oder einfach nur um überhaupt was zu sagen, oder einfach aus Frust weil man es nicht besser kann ?

Alles nur Anzeichen dafür das Du unreif bist.

Gruß Hagen

PS: ich hoffe das meine "Gegenkritik" an deinem Verhalten dich nun ermuntert mal 90 Grad den Kopf nach unten zu kippen und zu schauen wie lang er wirklich ist.

[Kara]

Ich hoffe, dass ich nichts sehe... :wink:

Und ich frage mich grad, ob du einen schlechten Tag hast, oder warum meine zwei Zeilen dich zu solchen Antworten treiben.
Aber bevor wir noch länger OT bleiben, ziehe ich mich geschlagen aus diesem Thread zurück und überlasse ihn wieder euren Expertengesprächen.