Sehr schneller Primzahl-Finder

[negaH]

@Phantom1:

so es is Wochenende und ich habe mir mal deinen Algo genauer angeschaut.
Vorweg: unsere beiden Siebe unterscheiden sich an wesentlichen Punkten, sind also nicht identisch, wenn auch die mathematische Grundlage identisch scheint.

Bevor wir aber weiter versuchen deinen Algo zu optimieren, müssen wir ihr erstmal korrekt lauffähig bekommen.
D.h. als erstes habe ich überprüft ob dein Algo korrekt arbeitet. Leider ist dies nicht der Fall.

Mein Testcode dazu ist:

Delphi-Quellcode:

function SavePrimes(MaxPrime: Cardinal): Cardinal;
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;
  Count,FailsIndex,FailsNonPrime,P,Q,Index, i, j, k, PrimeLen, PrimeBits, Num, Num2, m, mbit, s: Cardinal;
  QisPrime: Boolean;
  f: TextFile;
begin
  Index := 1;
  FailsIndex := 0;
  FailsNonPrime := 0;
  Count := 0;

  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;


     for I := 0 to PrimeLen-1 do
       for J := 0 to 7 do
       begin
         Q := k * PrimeBits + I * 30 + Stempel[J];
         if Q > MaxPrime then Break;
         if not ((I = 0) and (J = 0) and (K = 0)) and (Primes[I] and (1 shl J) = 0) then
         begin
           P := Prime.Primes[Index];
           if Q <> P then
           begin
             QisPrime := Prime.IsPrime(Q);
             Inc(FailsIndex);
             if not QisPrime then Inc(FailsNonPrime);
             WriteLn('Fails: ', FailsIndex:5, '/', FailsNonPrime:5,  ', Index: ', Index:12, ', Q: ', Q:12, ', ',QisPrime:6, ', P: ', P:12, ', ', Prime.IsPrime(P):6);
             while (Q > P) and (Index < Prime.Primes.MaxIndex) do
             begin
               Inc(Index);
               P := Prime.Primes[Index];
               if P <> Q then
               begin
                 Inc(FailsIndex);
                 WriteLn('Fails: ', FailsIndex:5, '/', FailsNonPrime:5,  ', Index: ', Index:12, ', Q: ', Q:12, ', ',QisPrime:6, ', P: ', P:12, ', ', Prime.IsPrime(P):6);
               end;
             end;
             WriteLn;
             if QIsPrime then Inc(Index);
           end else Inc(Index);
           Inc(Count);
         end;
       end;
  end;
  Result := count;
end;

Wie du siehtst habe ich das ganze Dateihandling rausgenommen, da es irrelevant ist. Dafür habe ich aber die Anzahl der Primzahlen integriert und später dann mein eigenes Sieb parallel arbeitend integriert, um die Resultate direkt vergleichen zu können.

Ein Ausschnitt des obigen Algo. in meiner Console sieht dann so aus:

Code:

Fails:     1/    0, Index:            1, Q:            7,   TRUE, P:            2,   TRUE
Fails:     2/    0, Index:            2, Q:            7,   TRUE, P:            3,   TRUE
Fails:     3/    0, Index:            3, Q:            7,   TRUE, P:            5,   TRUE

Fails:     4/    1, Index:    203233127, Q:   4293918781,  FALSE, P:   4293918787,   TRUE
Fails:     5/    2, Index:    203233128, Q:   4293918791,  FALSE, P:   4293918793,   TRUE
Fails:     6/    3, Index:    203233131, Q:   4293918877,  FALSE, P:   4293918883,   TRUE
Fails:     7/    4, Index:    203233132, Q:   4293918887,  FALSE, P:   4293918907,   TRUE
Fails:     8/    5, Index:    203233133, Q:   4293918913,  FALSE, P:   4293918953,   TRUE
Fails:     9/    6, Index:    203233133, Q:   4293918947,  FALSE, P:   4293918953,   TRUE

....

Fails: 38283/38280, Index:    203280214, Q:   4294967101,  FALSE, P:   4294967111,   TRUE
Fails: 38284/38281, Index:    203280214, Q:   4294967107,  FALSE, P:   4294967111,   TRUE
Fails: 38285/38282, Index:    203280219, Q:   4294967213,  FALSE, P:   4294967231,   TRUE
Fails: 38286/38283, Index:    203280220, Q:   4294967273,  FALSE, P:   4294967279,   TRUE

....

Fails: 38287/38284, Index:    203280222, Q:           15,  FALSE, P:            0,  FALSE
Fails: 38288/38285, Index:    203280222, Q:           21,  FALSE, P:            0,  FALSE
Fails: 38289/38285, Index:    203280222, Q:           37,   TRUE, P:            0,  FALSE
Fails: 38290/38285, Index:    203280223, Q:           61,   TRUE, P:            0,  FALSE
Fails: 38291/38286, Index:    203280224, Q:           75,  FALSE, P:            0,  FALSE
Fails: 38292/38287, Index:    203280224, Q:           81,  FALSE, P:            0,  FALSE

....

Fails: 113112/105327, Index:    203288004, Q:       917403,  FALSE, P:            0,  FALSE
Fails: 113113/105328, Index:    203288004, Q:       917425,  FALSE, P:            0,  FALSE
Fails: 113114/105329, Index:    203288004, Q:       917433,  FALSE, P:            0,  FALSE
Fails: 113115/105330, Index:    203288004, Q:       917463,  FALSE, P:            0,  FALSE
Fails: 113116/105331, Index:    203288004, Q:       917467,  FALSE, P:            0,  FALSE
Fails: 113117/105332, Index:    203288004, Q:       917497,  FALSE, P:            0,  FALSE

Aufruf: mit SavePrimes($FFFFFFFB) -> $FFFFFFFB = 4294967291 ist höchste Primzahl < 2^32 mit Primzahlindex 203280221.

Fails: gibt ab wieviele Fehlresultat der Algo macht -> Anzahl Indexfehler/Anzahl fehlerhafter Primzahlen (sind also zusammengesetzte Zahlen die dein Algo als Primzahlen ausgibt).
Index: der Index der Primzahl in der Primzahltabelle.
Q: die Zahl die dein Algo als Primzahl ausgibt, danach FALSE/TRUE je nachdem ob Q tatsächlich eine Primzahl ist.
P: die Primzahl die mein Algo. zum Index berechnet, FALSE/TRUE jenachdem ob P eine Primzahl ist


Die ersten 3 Zeilen können wir ignorieren, da dein Algo keine direkte Berechnung zu den ersten 3 Primzahlen bietet.

Ab Primzahlindex 203233127 erzeugt dein Algo nun Zahlen die keine Primzahlen sind, er arbeitet ab da falsch.
Das geht bis Primzahlindex 203280220 was exakt der Moment ist wo der Algo im Grunde terminieren müsste.

Bis zu diesem Bereich hat dein Algo also 38287 Zahlen als Primzahlen erzeugt die aber zusammengesetzte Zahlen sind.
Der letzte Block ab dem P= 0 ist, hätte dein Algo bei der Zeile

Delphi-Quellcode:

     for I := 0 to PrimeLen-1 do
       for J := 0 to 7 do
       begin
         Q := k * PrimeBits + I * 30 + Stempel[J];
         if Q > MaxPrime then Break;                    // <---

schon längst terminieren müssen. Er rechnet aber weiter da wie man sieht Q nun < $FFFFFFFB ist, sprich ein Integer Überlauf dazu führt das der Algo noch mehr falsche Zahlen berechnet. Da er aber nicht terminiert beendet sich der Algo erst mit der äußersten Schleife k und erzeugt somit 113117 falsche Antworten.


Nun zur Performance:

ich habe mit nachfolgendem Source getestet. Auch hier wieder die Dateioperationen und Stringkonvertierungen entfernt, dafür aber unterschieden ob man nur die Anzahl der Primzahlen oder auch die Primzahlen wertmäßig exakt berechnen möchte. Diese Unterscheidung ist wichtig beim direkten Vergleich der unterschiedlichen Verfahren.

Delphi-Quellcode:

function SavePrimes1(MaxPrime: Cardinal; CalcPrimes: Boolean): Cardinal;
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;
  Count, i, j, k, PrimeLen, PrimeBits, Num, Num2, m, mbit, s: Cardinal;
begin

  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 CalcPrimes 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
             Inc(Count);
         end;
  end;
  Result := count;
end;

procedure TestPrimes;
var
  I,P: Cardinal;
begin
  StartTimer;
  Primes.IndexOf($7FFFFFFF);
  Writeln(Secs:10:2, ' sec');

  StartTimer;
  SavePrimes1($7FFFFFFF, False);
  Writeln(Secs:10:2, ' sec');

  I := 1;
  StartTimer;
  repeat
    P := Primes[I];
    Inc(I);
  until P >= $7FFFFFFF;
  Writeln(Secs:10:2, ' sec');


  StartTimer;
  SavePrimes1($7FFFFFFF, True);
  Writeln(Secs:10:2, ' sec');
end;

Laufzeiten:

Code:

Primes.IndexOf($7FFFFFFF)         =  12.74 sec
SavePrimes($7FFFFFFF, False)      =  51.38 sec

loop Primes[i] until >= $7FFFFFFF =  30.64 sec
SavePrimes($7FFFFFFF, True)       =  67.98 sec

auf einem P4 1.5GHz 512Mb und Delphi 5.

Die gravierenden Performanceunterschiede zum AMD sind für mich nur durch die Anwendung der Opcodes BTS/BTC in meinen Bitarray[] Funktionen erklärbar. AMD Prozessoren sind mit diesen Operationen wesentlich langsammer als Pentium CPU's. Für AMD's müsste man diese Opcodes durch indizierte Array[] Zugriff + AND/OR Bitmasken per SHL/SHR ersetzen.

Alledings, als erstes muß man den Algorithmus sauber lauffähig bekommen, denn Speed nutzt nur dann etwas wenn die berechneten Resultate auch wirklich mathematisch korrekt sind.

Gruß Hagen