|
Registriert seit: 25. Jun 2003
Ort: Thüringen
2.950 Beiträge
|
Re: Sehr schneller Primzahl-Finder
26. Aug 2005, 22:12
@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
|