Entwickler-Ecke
Algorithmen, Optimierung und Assembler - Ein Algorithmus zur Primzahlenenumeration
Delphi-Laie - So 17.06.12 17:02
Titel: Ein Algorithmus zur Primzahlenenumeration
Eine (einfache oder komplexere) Formel zur Erzeugung der Primzahlen, optimalerweise der x. Primzahl, soll es ja bis heute nicht geben. Wohl gibt es aber inzwischen von Gerald Bühler
einen siebfreien (i.Ggs.z. Eratosthenes, Atkin) Algorithmus [
http://www.geraldbuehler.de/primzahlen/] zur Generierung derselben, der allerdings mit einem Hiflsarray arbeitet und bei größeren Zahlen auch spürbar langsamer wird. Das geht aber vom algorithmischen Niveau her immer noch um Größenordnungen über die triviale und von der Laufzeit her inakzeptable Schulmethode heraus, über Probedivisionen aller ungeraden Zahlen (und ggf. Einbezug noch weiterer Teilbarkeitsregeln) bis hin zu den Wurzeln der Zahlen zu operieren und so die Primeigenschaft nach der Versuch-Und-Irrtum-Methode festzustellen. Das wird bei wachsenden Eingabegrößen / -zahlen recht schnell ziemlich zäh.
Interessanterweise ist jener Algorithmus erweiterbar, wie ich hoffte, probierte und herausfand, sodaß man sich nicht vorher festlegen muß, bis zu welcher natürlichen Zahl n (oder alternativ bis zur Primzahl bis zur Nummer soundso) man enumerieren möchte (im Gegensatz zu den Sieben). Es sind also beliebig viele Primzahlen aufsteigend lückenlos generierbar bzw. enumerierbar.
Im Anhang befinden sich 3 Demonstrationsprogramme dazu, die ich auf die Schnelle erstellte. Primgen 1 ist das nach Delphi übersetzte Original, während Primgen 2 und 3 die nach oben offene Enumeration demonstrieren.
Edit: Noch eine kleine Unsauberkeit korrigiert.
Mathematiker - So 17.06.12 19:18
Hallo Delphi-Laie,
die drei Programme sind interessant, allerdings beim Zählen der Primzahlen nicht wirklich schnell.
Wesentlich schneller ist folgende rekursive Methode, die vor allem auf größere Datenmengen verzichtet.
Delphi-Quelltext
1:
2:
3:
4:
5:
6:
7:
8:
9:
10:
11:
12:
13:
14:
15:
16:
17:
18:
19:
20:
21:
22:
23:
24:
25:
26:
27:
28:
29:
30:
31:
32:
33:
34:
35:
36:
37:
38:
39:
40:
41:
42:
43:
44:
45:
46:
47:
48:
49:
50:
51:
52:
53:
54:
55:
56:
57:
58:
59:
60:
61:
62:
63:
64:
65:
66:
67:
68:
69:
70:
71:
| Function Primzahlanzahl( N:integer ):integer;
var r, p:integer;
function isprime(n:integer):Boolean;
var i,w:integer;
teiler:boolean;
begin
if n = 1 then begin
isprime:=false; exit;
end;
if n in [2,3,5,7,11,13] then begin
isprime:=true; exit;
end;
if odd(n) then begin
if n mod 3 = 0 then begin
isprime:=false; exit;
end;
if n mod 5 = 0 then begin
isprime:=false; exit;
end;
if n mod 7 = 0 then begin
isprime:=false; exit;
end;
if n mod 11 = 0 then begin
isprime:=false; exit;
end;
i:=13;
w:=trunc(sqrt(n)+1);
if w<13 then
begin
isprime:=true;
exit;
end;
repeat
teiler:=(n mod i = 0);
inc(i,2);
until teiler or (i>w);
isprime:=not teiler;
end
else isprime:=false;
end;
Function Entfernen( n,p: integer):integer;
var nklein, p2, d: integer;
begin
nklein := n div p;
if nklein<p then
begin
Entfernen:=1;
exit
end;
d := nklein;
for p2:=2 to p-1 do
if IsPrime(p2) then
d := d - Entfernen(nklein,p2);
Entfernen:=d
end;
begin
r := N-1;
for p:=2 to n do
if p*p>n then
begin
Primzahlanzahl := r;
exit
end
else
if IsPrime(p) then
begin
r := r - Entfernen(n,p);
r := r+1;
end;
Primzahlanzahl := r;
end; |
Bei einem Test auf meinem Rechner benötigte dein Programm bis 10 Millionen 3,74 s, meine Routine 0,11 s. Beim Zählen bis 100 Millionen wurde es noch dramatischer 75,74 s zu 0,89 s.
Zu meinem Text: Die scheinbar umständliche Funktion isprime sorgt dafür, dass möglichst schnell kleine Primteiler ausgesondert werden. Wahrscheinlich kann man isprime durch eine effizientere Funktion ersetzen.
Beste Grüße
Mathematiker
Horst_H - So 17.06.12 21:53
Hallo,
es kann doch kaum einer gegen die Version des Kim Walisch
http://code.google.com/p/primesieve/ anstinken, der jüst eine neue Version hat 3.7 ;-)
| Zitat: |
Bis 1e9 in 0,09 Sekunden für einen
Intel Core i7-920 4 x (2.66GHz, 32K L1 Data Cache) |
Und soviel Platz braucht die auch nicht, sonst wäre sie nicht so schnell.
Ich die Version 1 jetzt mal bis 1e8 rechnen lassen und die belegt dabei maximal knapp über 50 Mb , das ist wohl nicht sehr sparsam.
Gelegenheit mal die PowerShell und Measure-Command zu testen.
Quelltext
1:
2:
3:
4:
5:
6:
7:
8:
9:
10:
11:
12:
13:
| PS ...\Primgen 1> Measure-Command {start-process Project1 -Wait}
Days : 0
Hours : 0
Minutes : 0
Seconds : 45
Milliseconds : 590
Ticks : 455905584
TotalDays : 0,0005276685
TotalHours : 0,012664044
TotalMinutes : 0,75984264
TotalSeconds : 45,5905584
TotalMilliseconds : 45590,5584 |
Ich habe 1e8 per Copy und paste eingefügt und mich auch sonst beeilt.
Also 40 Sekunden werden es gewiß gewesen sein. ;-)
Das ist jetzt nicht wirklich überzeugend, wozu wurde der Algorithmus notariell bestätigt?
Gruß Horst
Selbst mein Uralt Programm ist wesentlich schneller.
0,3 Sekunden bis 1e8 bei 2,6 Mb davon 835kb privater Speicher.
SIeve of Atkin in seiner völlig unoptimierten, ursprünglichen Form, die mir Fiete zeigte, braucht über 5 Sekunden bei 100 Mb Speicherbedarf.
alzaimar hatte dieses hochoptimiert, nach einer Vorlage in C, zu Delphi übertragen
http://www.entwickler-ecke.de/topic_Optimierung+einer+Primzahlfunktion_33013,232.html
-> PrimeGenDelphi.zip
http://www.entwickler-ecke.de/download.php?id=1364
Delphi-Laie - So 17.06.12 22:54
Danke für Eure Reaktionen!
Daß jener Algorithmus nicht die beste Komplexität hat, ging aus meinem ersten Beitrage ja schon hervor. Dafür ist er insofern eleganter, als daß er sich von der naheliegenden und seit der Vorantike bekannten Versuch-Und-Irrtum-Methode der Probedivisionen entfernt / unterscheidet. Die Siebungen, wenn es einem nur um Geschwindigkeit geht, sind ja noch viel schneller, doch speicherverschwendend und eben leider nicht enumerierend: Man muß sich vorher festlegen bzw. beschränken, und danach ist Schluß bzw. man muß mit einem größeren Startwert neu beginnen, während eine solche Enumeration eben beliebig wiederholt bzw. fortgesetzt werden kann. Bei großen Startmengen wird der Speicherverbrauch der Startmengen für die Siebungen sogar für GByte-Computer deutlich bis heftig (die speichersparenden Datenkomprimierungen im "Marktführer" jetzt mal außer Betracht gelassen, die Tendenz ist dennoch gegeben).
Ohne etwas in den Himmel zu loben, was ohnehin im wesentlichen nicht von mir stammt: Konkret die x. Primzahl (ohne über den Umweg der Enumeration) zu generieren, ist auch dieser Algorithmus leider nicht imstande. Für Primzahlzerlegungen und -faktorisierungen hingegen liefert er in Einzelschüssen jedoch genau die saubere Startmenge, die man dafür benötigt.
Horst_H - So 17.06.12 23:18
Hallo,
nichtsdestotrotz finde ich bei einer Speicherbelegung von 50 Mb um alle Zahlen bis 1e8 mit Primgen 1 zu bestimmen, völlig gegensätzlich zu der Aussage von einem platzsparendem Algorithmus.
Gruß Horst
Delphi-Laie - Mo 18.06.12 10:05
 Horst_H hat folgendes geschrieben  : |
| nichtsdestotrotz finde ich bei einer Speicherbelegung von 50 Mb um alle Zahlen bis 1e8 mit Primgen 1 zu bestimmen, völlig gegensätzlich zu der Aussage von einem platzsparendem Algorithmus. |
Meinst Du das beispielhafte C++-Programm? Zugegebenermaßen, als ich Eratosthenes las, war mein Interesse schon wieder erloschen.
Die Sache interessierte mich nun aber doch weiter: Also, die Summe der Länge beider Arrays (des eigentlichen, in dem die Primzahlen abgelegt werden und die des Hilfsarrays) beträgt nur einen Bruchteil dessen, als wenn man ein Array einfach mit allen natürlichen Zahlen füllen würde und Erathostenes darauf losließe. Gut, bei Eratosthenes muß man nur einen Wert, nämlich boolean, ablegen, das reicht ja als Information über die Primeigenschaft. Ob das von Delphi speichersparend umgesetzt wird, weiß ich nicht, vermute aber, daß für jeden Wert des Types "boolean" analog zu "bool" letztlich doch eine integre Variable "verschwendet" wird. Bitweiser Zugriff auf Integertypen wäre ja auch langsamer.
Dennoch ist die Tendenz eindeutig: Der Speicherbedarf (zumindest i.S. der Summe der Längen beider Arrays) bleibt gegenüber der Länge des Eratosthenes-Arrays immer mehr zurück, je mehr man enumeriert. Er beträgt nur einen (eben immer kleiner werdenden) Bruchteil.
Mathematiker - Mo 18.06.12 11:46
Hallo Hort_H,
| Horst_H hat folgendes geschrieben : |
| alzaimar hatte dieses hochoptimiert, nach einer Vorlage in C, zu Delphi übertragen |
Ich habe mir Deinen Hinweis auf PrimeGen angesehen und kann nur sagen: faszinierend.
Aber! Irgendwie stelle ich mich zu dämlich an. :autsch:
Ich schaffe es weder die Anzahl der Primzahlen bis zu einer genauen(!) Grenze, sagen wir 5 Millionen, berechnen zu lassen noch irgendeine Primzahl (zum Beispiel die letzte berechnete) anzeigen zu lassen.
Aus dem Quelltext werde ich einfach nicht schlau. Die Prozedur berechnet offensichtlich die Primzahlen:
Delphi-Quelltext
1:
2:
3:
4:
5:
6:
7:
8:
9:
10:
11:
12:
13:
14:
15:
16:
17:
18:
19:
20:
21:
22:
23:
24:
25:
26:
27:
28:
| Procedure DoIt;
Var
i,j : Integer;
Buf : pABuffer;
pq, Nexti : pCardinal;
k,q,kPos,bits,data : Cardinal;
Begin
fillChar (a,sizeOf (a), 0);
For i:=0 to 7 do begin
Nexti := @next[i];
Buf := @a[i];
pq := @qtab[0];
For j:=0 to 3508 do Begin
k := nexti^;
if (k < ccBlockLen) Then Begin
q := pq^;
Repeat
kPos := k shr 5;
buf^[kPos] := buf^[kPos] or (1 shl (k and 31));
inc (k,q);
Until k >= ccBlockLen;
End;
nexti^ := k - ccBlockLen; inc (pq);
inc (nexti);
End;
End;
End; |
Aber, wo bitte stehen diese dann. Die vielen Zeiger und Adresszuweisungen sind für mich kryptisch.
Kennst Du ein Beispiel, wie man PrimeGen vernünftig einsetzen kann?
Danke.
Mathematiker
Horst_H - Mo 18.06.12 12:54
Hallo,
@Delphi-Laie:
Wie kommt es denn nun, das Dein PrimeGen 1 maximal 50 Mb=5e7 Byte belegt um die Primzahlen bis 1e8 zu zählen.Das ist so groß wie das Sieb des Eratosthenes mit ungeraden Zahlen.Sind das nur Datenleichen auf dem Weg nach oben?
@Mathematiker:
alzaimar hat damals, 2005, das Programm primeGen
http://cr.yp.to/primegen.html nach Delphi übertragen, das ist nicht auf seinem Mist gewachsen, was er auch nie behauptet hat.Ich habe doch in dem thread auch gefragt, wie man den jetzt zu einer bestimmten Zahl zählen muss ;-)
Mal schauen, ob ich drauf komme.
Gruß Horst
Gammatester - Mo 18.06.12 12:54
Einen Algorithmus der aller Primzahlen für bis zu einer bestimment großen Zahl bestimmt braucht man wohl eher nur für theroetische Zwecke, mal abgesehen vom Speicher (Problem 50MB siehe oben), der praktisch linear anwächst (ca n / ln(n)). Zwei Aufgaben, die in der Praxis vorkommen, sind u.a.
1. Ein Primzahlgenerator bzw. -Enumerator, der fortlaufende Primzahlen erzeugt ohne diese zu speichern.
2. Ein segmentiertes Sieb, daß zB die auf eine Zahl folgenden 1000 oder 10000 Primzahlen erzeugt und speichert (zB für Primzahlen im RSA-Algorithmus)
Der Generator aus 1. könnte theoretisch ohne Speicher auskommen, ist dann aber langsam. Praktisch wird er auch eine Art segmetiertes Sieb genutzen. Für das Problem 5000000-te Primzahl liefert mein Generator aus
MPArith [
http://www.wolfgang-ehrhardt.de/misc_de.html#mparith] für 32-Bit-Integer
Quelltext
1:
2:
3:
| X:\TMP>T_PS.EXE
Zeit [s]: 0.214
5000000-te Primzahl = 86028121 |
was man leicht mit Wolfram Alpha verifizieren kann. Der Quellcode dazu:
Delphi-Quelltext
1:
2:
3:
4:
5:
6:
7:
8:
9:
10:
11:
12:
13:
| uses
mp_prime;
var
sieve: TSieve;
np,n,p: longint;
begin
prime_sieve_init(sieve,2);
np := 5000000;
for n:=1 to np do begin
p := prime_sieve_next(sieve);
end;
writeln(np,'-te Primzahl = ',p);
end. |
Delphi-Laie - Mo 18.06.12 13:38
| Horst_H hat folgendes geschrieben : |
@Delphi-Laie:
Wie kommt es denn nun, das Dein PrimeGen 1 maximal 50 Mb=5e7 Byte belegt um die Primzahlen bis 1e8 zu zählen.Das ist so groß wie das Sieb des Eratosthenes mit ungeraden Zahlen.Sind das nur Datenleichen auf dem Weg nach oben? |
Also, "mein" Primgen ist es genaugenommen nicht. Ich übertrug im wesentlichen nur den Algorithmus nach Delphi. Auch kann ich weder zum Speicherverbrauch noch über irgendwelche "Datenleichen" i.S. einer Begründung Auskunft geben. Ich bin schließlich nicht das Programm selbst. Bei mir benötigt das Programm bei n= ca. 5 Mio. nur knapp 4 MByte im Speicher - die Tendenz ist natürlich allmählich steigend. Unelegant ist in meinen Augen höchstens das ständige Anpassen der Arraylängen, doch auch nur, seitdem ich weiß, daß das die Speicherfragmentierung forciert (was aber nach Programmende wieder behoben sein dürfte).
Was ich jedoch kann, ist, die Länge der beiden Arrays (konkret: deren Summe) mit der natürlichen Zahl n, bei der der Algorithmus gerade werkelt, vergleichen und ins Verhältnis setzen. Bei größeren Werten von n, im Bereich von einigen Millionen, sinkt er auf unter 1/7 der Zahl n ab. n wäre die Länge des Eratosthenes-Arrays.
Ergänzungen:
1. Bei 100 Mio. steigt der Speicherbedarf "meines" obigen Programmes tatsächlich bis auf ca. 50 MByte an.
2. Ein bis 100 Mio. dimensionertes Cardinalarray erfodert jedoch knapp 400 MByte!
3. Ein bis 100 Mio. dimensionertes Wordarray erfodert knapp 200 MByte.
4. Ein bis 100 Mio. dimensionertes Bytearray erfodert - nunmehr nicht überraschend - ca. 100 MByte.
Also werden Integervariablen anscheinend speichersparsam abgelegt.
5. Ein bis 100 Mio. dimensionertes Booelanarray erfodert ebenfalls ca. 100 MByte. Anscheinend wird tatsächlich pro booleschen Wert ein Byte benötig.
Mithin ist "mein" obiges Programm eben doch speichersparsamer. Allerdings war mein simpler Vergleich der Arraylängen zu oberflächlich, da ein Sieb nur mit Boolean gefüllt werden muß: Es bleiben also 50 MByte zu 100 MByte, also etwa die Hälfte.
Und noch eine Ergänzung: Wenn man Eratosthenes nur mit ungeraden Zahlen füllt, ja, dann ist tatsächlich Gleichstand erreicht!
Delphi-Laie - So 24.06.12 18:28
Zurecht wurde hier kritisiert, daß der Speicherverbrauch dieses hier vorgestellten Algorithmus' einem Eratoshenes-Sieb, das nur ungerade Zahlen speichert, gleicht.
Ich nahm mich der Sache an. Schon in der ersten Version fiel mir auf, daß das Hilfsarray überwiegend (in der weit überwiegend höheren und damit größeren "Hälfte") aus Nullen besteht. Tatsächlich wird es an falscher Stelle des Algorithmus' falsch dimensioniert, konkret: natürlich überdimensioniert. Jetzt müßte es jedoch stimmen. Ich lud bei den drei Demonstrationsprogrammen die jeweils verbesserte Version oben hoch.
Nach wie vor bin ich der Meinung, daß dieser Algorithmus ein höheres (algorithmisches) Nivaus als jedes noch so ausgeklügeltes Sieb darstellt, weil er enumeriert und damit - ohne Neuanfang - "nach oben offen" ist.
Viele Grüße
Delphi-Laie
Edit: Nach dieser Korrektur wächst die Länge des Hilfsarrays anscheinend nur noch logarithmisch mit der Anzahl der ermittelten Primzahlen, fällt also bei wachsender Primzahlmenge bwz. wachsendem Primzahlarray immer weniger bis kaum noch ins Gewicht.
Gammatester - Mo 09.07.12 10:11
| Mathematiker hat folgendes geschrieben : |
... Bei einem Test auf meinem Rechner benötigte dein Programm bis 10 Millionen 3,74 s, meine Routine 0,11 s. Beim Zählen bis 100 Millionen wurde es noch dramatischer 75,74 s zu 0,89 s. Zu meinem Text: Die scheinbar umständliche Funktion isprime sorgt dafür, dass möglichst schnell kleine Primteiler ausgesondert werden. Wahrscheinlich kann man isprime durch eine effizientere Funktion ersetzen.
|
Hallo Mathematiker,
entschuldige, daß ich nochmal auf diesen Beitrag zurückkomme. Ich habe mal Deine Prozedur ein wenig optimiert, Zeiten auf meinem Rechner:
Quelltext
1:
2:
| bis 100 Millionen: Dein Original 0.318s, meine V3: 0.048
bis 1000 Millionen: Dein Original 3.437s, meine V3: 0.422 |
Also ein relativ schneller Code für die Funktion PrimePi(n).
Frage: Hast Du den Entfernen-Algorithmus selbstentwickelt oder gibt es dazu Quellen, was genau macht Entfernen(n,p)?
Gruß Gammatester
Delphi-Laie - Mo 09.07.12 10:16
Ich beschäftigte mich mit dem Algorihtmus noch einmal näher.
Also, das "Rumgemache" mit der Quadratwurzel kann ebensogut entfallen. Und dann lassen sich die inkrementierte Variablen (c und i) auch um 2 inkrementieren, um gleich alle geraden Zahlen auszuschließen. Die 2 als erste bzw. kleinste und einzige gerade Primzahl muß dann besonders berücksichtigt werden.
Immerhin, und das ist ja das Ziel des originalen Programmautors, entfallen die Probedivisionen. Inwieweit das wirklich einen Geschwindigkeitsvorteil bringt, ist mir nicht bekannt. Div und mod bei den Integern sind ja auch recht schnell. Die Geschwindigkeit dieses Algorithmus' mit der eines Siebes zu vergleichen, bringt ihn unweigerlich ins Hintertreffen. Der Vorteil liegt gegenüber dem der Probedivisionen und hat ihm gemein, eine nach oben offene Enumerierbarkeit zu ermöglichen, die bei Sieben nur "pseuso" möglich ist, indem die "voreingestellte" Siebgröße erhöht wird und dann erneut gesiebt wird. Außerdem entfällt der hohe Speicherbedarf, den die Siebe mit sich bringen.
Mathematiker - Mo 09.07.12 10:47
Hallo Gammatester,
| Gammatester hat folgendes geschrieben : |
Ich habe mal Deine Prozedur ein wenig optimiert, Zeiten auf meinem Rechner:
Quelltext 1:
2:
| bis 100 Millionen: Dein Original 0.318s, meine V3: 0.048
bis 1000 Millionen: Dein Original 3.437s, meine V3: 0.422 | Also ein relativ schneller Code für die Funktion PrimePi(n). |
Klingt sehr gut. Es wäre nett, wenn Du die Änderungen veröffentlichen würdest.
| Gammatester hat folgendes geschrieben : |
| Hast Du den Entfernen-Algorithmus selbstentwickelt oder gibt es dazu Quellen, was genau macht Entfernen(n,p)? |
Nein, der ist nicht von mir. Vor einigen Jahren, etwa 1995-2000, habe ich den Algorithmus in einer Zeitschrift (oder im Netz?) gesehen. Ich weiß nicht mehr wo und kann Dir auch nicht genau sagen, von wem er stammt. Ich glaube mich zu errinnern, dass er von einem der Borweins war. :nixweiss: Auf jeden Fall war er in C geschrieben.
Verstanden habe ich Ihn auch nicht so richtig. Aber er funktioniert erstaunlicher Weise genau.
Beste Grüße
Mathematiker
Gammatester - Mo 09.07.12 10:50
@Delphi-Laie,
schon klar, meine Frage bezog sich wirklich nur konkret auf des Mathematikers Funktion Entfernen. Hier meine Implementation der PrimePi-Funktion (setzt MPArith voraus):
Delphi-Quelltext
1:
2:
3:
4:
5:
6:
7:
8:
9:
10:
11:
12:
13:
14:
15:
16:
17:
18:
19:
20:
21:
22:
23:
24:
25:
26:
27:
28:
29:
30:
31:
32:
33:
34:
35:
36:
37:
38:
39:
40:
41:
42:
43:
44:
45:
46:
47:
48:
49:
50:
51:
52:
53:
54:
55:
56:
57:
58:
59:
60:
61:
62:
63:
64:
65:
66:
67:
68:
69:
70:
71:
72:
73:
74:
75:
| program t_psv3;
{$i STD.INC}
{$i mp_conf.inc}
{$x+}
{$i+}
{$ifdef APPCONS}
{$apptype console}
{$endif}
{$ifndef FPC}
{$N+}
{$endif}
uses
HRTimer,
mp_prime;
function remove(n,p: longint):longint;
var
d: longint;
k: integer;
w: word;
begin
n := n div p;
if n<p then remove := 1
else begin
d := n;
for k:=1 to NumPrimes16 do begin
w := Primes16[k];
if w+1 > p then break;
d := d - remove(n,w);
end;
remove := d;
end;
end;
function primepi(n: longint): longint;
var
r: longint;
k: integer;
p: word;
w: word;
begin
r := n-1;
if n > 2 then begin
w := trunc(sqrt(n));
for k:=1 to NumPrimes16 do begin
p := Primes16[k];
if p > w then break;
r := succ(r) - remove(n,p);
end;
end;
primepi := r;
end;
var
np,n: longint;
HR: THRTimer;
begin
writeln('-----------------------');
np := 2100000000;
StartTimer(HR);
n := primepi(np);
writeln('Zeit [s]: ',ReadSeconds(HR):1:3);
writeln('V3 PrimePi(',np,') = ',n);
end. |
Die Geschwindigkeit [PrimePi(2100000000) = 102886526 in 0.855s] gefällt mir so gut, daß ich beabsichtigte, sie in die nächste MPArith aufzunehmen, selbstverständlich mit Referenz an Mathematiker und Entwicklerecke wie im Source.
Dazu wäre es natürlich angebracht, genau zu beschreiben, was Entfernen/Remove macht.
Gruß Gammatester
Horst_H - Mo 09.07.12 11:10
Hallo,
was remove macht, ist doch klar.
n div p ->Anzahl aller Vielfachen von p (Faltoren k = 1,2,3..n div p)
Das sind die theoretischen Ausstreichungen durch p und davon zieht ( remove ) man die ab, die durch vorangegangenen Primzahlen schon gestrichen wurden.
Gruß Horst
Mathematiker - Mo 09.07.12 12:05
Hallo,
ich habe die Quelle gefunden.
Irgendwie werde ich alt, dass ich mir nicht mehr merken kann, von wem etwas war. Vielleicht ist es auch BSE? :lol:
Beste Grüße
Mathematiker
Nachtrag: Original-Datei entfernt!
Gammatester - Mo 09.07.12 12:28
| Horst_H hat folgendes geschrieben : |
..was remove macht, ist doch klar.
n div p ->Anzahl aller Vielfachen von p (Faltoren k = 1,2,3..n div p)
Das sind die theoretischen Ausstreichungen durch p und davon zieht ( remove ) man die ab, die durch vorangegangenen Primzahlen schon gestrichen wurden. |
Soweit klar. Nur das 'vorangegangen' durch die Implementation festgelegt wurde, man könnte die Schleife in PrimePi auch rückwärts laufen lassen. Hältst Du diese Charakterisierung für richtig?
"Remove(n,p) mit p prim <= n zählt die Vielfachen k*p <= n, die nicht Vielfache der Primzahlen p_j < p sind."
Gruß Gammatester
Horst_H - Mo 09.07.12 15:14
Hallo,
| Zitat: |
Hältst Du diese Charakterisierung für richtig?
"Remove(n,p) mit p prim <= n zählt die Anzahl der Vielfachen k*p <= n, die nicht Vielfache der Primzahlen p_j < p sind." |
Im Endeffekt ist es genau das.
Man kann es auch anschaulich so formulieren,
Primzahl P sagt ich streiche nsmall (n DIV P ) Zahlen und fragt jetzt die alle kleineren, von klein zu groß hin, wieviele sie davon schon gestrichen haben und zieht diese von seiner Zahl ab.
Und die fragen ihrerseits genauso.
remove selber gibt durch Rekursion die Anzahl, der nur durch die Primzahl p_j zuerst ausgestrichen Zahlen an ( also ist p_j ist der kleinste Faktor der gestrichenen Zahl ).
In pi.pas ist das ja als dismissed bezeichnet, wenn google translate recht hat, soll man diese Zahlen dann als die durch p_j verworfen ( nicht prim ) betrachten.
Ich dachte zuerst, wie kann das je fertig werden wenn Anzahl(Primzahlen bis p)*Anzahl (Primzahlen bis p)/2 Aufrufe getätigt werden, aber das sind für Dword auch nur 6543*6543/2 ~ 21 Milliionen , zudem landet man ja schnell bei extrem kleinen Zahlen nsmall, sodass die Rekursionstiefe nicht 6542 ist.
Da 2 die kleinste Primzahl, löst aber keine Rekursion mehr aus, also ist und 65521? ( ~3^10) für DWord die Größte, ist maximale die Rekursionstiefe 10 oder 11.
Delphi-Quelltext
1:
2:
3:
4:
5:
6:
7:
8:
9:
| nsmall := n div p;
if nsmall<p then begin
DismissedBy :=1;
exit
end;
d := nsmall;
for p2:=2 to p-1 do
if IsPrime(p2) then
d := d - DismissedBy( nsmall, p2 ); |
Gruß Horst
Gammatester - Mo 23.07.12 09:37
Nach zwei Wochenenden habe ich die Rechenzeit auf weniger als 1/10 der schon 'guten' Zeit gedrückt. Allerdings basieren die neuen Routinen nicht mehr auf dem Mathematiker-Code (bzw. dem von Albert van der Horst), sondern auf den Algorithmen von Legendre/Meissel/Lehmer beschrieben in dem lesenswerten
Originalartikel [
http://projecteuclid.org/euclid.ijm/1255455259] von D.H. Lehmer aus dem Jahre 1959 (siehe auch
Wiki [
http://en.wikipedia.org/wiki/Prime-counting_function#Algorithms_for_evaluating_.CF.80.28x.29], gibt's warum auch immer nicht auf Deutsch).
Da MPArith V1.21 noch einige Zeit in Arbeit ist, habe ich im Anhang schon mal die neuen Funktionen im Quellcode für MPArith V1.20 und eine Delphi6-EXE eingepackt. Auf meinem Rechner hat man diese Zeiten:
Quelltext
1:
2:
3:
4:
5:
6:
7:
8:
9:
| PrimePi mit MPArith V1.20.24 (31/32 bit) (c) 2012 Mathematiker/WE(Gammatester)
CPUFrequenz: 2.826 GHz
Mathemetiker Originalcode [s]: 6.805 PrimePi(2000000000) = 98222287
Mathemetiker mit isprime32 [s]: 3.465 PrimePi(2000000000) = 98222287
Mathemetiker / WE Primes16 [s]: 0.937 PrimePi(2000000000) = 98222287
v.d. Horst / WE optimiert [s]: 0.728 PrimePi(2000000000) = 98222287
Legendre: WE mit lsumphi32 [s]: 0.371 PrimePi(2000000000) = 98222287
Meissel: WE mit lsumphi32 [s]: 0.099 PrimePi(2000000000) = 98222287
Lehmer: WE mit lsumphi32 [s]: 0.051 PrimePi(2000000000) = 98222287 |
Gruß Gammatester
Horst_H - Mo 23.07.12 11:14
Hallo,
beeindruckend schneller!
Aber AMD dividiert nicht gern, besonders auffällig bei dem Originalcode von Mathematiker.
Quelltext
1:
2:
3:
4:
5:
6:
7:
8:
9:
10:
11:
| PrimePi mit MPArith V1.20.24 (31/32 bit) (c) 2012 Mathematiker/WE(Gammatester)
CPUFrequenz: 3.191 GHz
Mathemetiker Originalcode [s]: 14.937 PrimePi(2000000000) = 98222287
Mathemetiker mit isprime32 [s]: 3.289 PrimePi(2000000000) = 98222287
Mathemetiker / WE Primes16 [s]: 1.178 PrimePi(2000000000) = 98222287
v.d. Horst / WE optimiert [s]: 0.902 PrimePi(2000000000) = 98222287
Legendre: WE mit lsumphi32 [s]: 0.754 PrimePi(2000000000) = 98222287
Meissel: WE mit lsumphi32 [s]: 0.168 PrimePi(2000000000) = 98222287
Lehmer: WE mit lsumphi32 [s]: 0.060 PrimePi(2000000000) = 98222287
Weiter mit Enter |
Vielleicht kann man das Prinzip von TabPhiA statt 5 = 2*3*5*7*11=2310 auf 6 == 30030 erweitern.Hauptsache es bleibt im Level_1 Cache der CPU mit 7508*2= 15016 Byte.
Delphi-Quelltext
1:
2:
3:
4:
5:
6:
7:
8:
9:
| const
BaseA = 6;
ProdA = 30030;
PrA_2 = 15015;
PrA_4 = 7508;
PhiPA = 5760;
TabPhiA: array[0..PrA_4] of WORD = (... |
Gruß Horst
Gammatester - Mo 23.07.12 11:41
| Horst_H hat folgendes geschrieben : |
Vielleicht kann man das Prinzip von TabPhiA statt 5 = 2*3*5*7*11=2310 auf 6 == 30030 erweitern.Hauptsache es bleibt im Level_1 Cache der CPU mit 7508*2= 15016 Byte.
|
Wahrscheinlich bringt das nur Wesentliches für die Algorithmen Legendre und Meissel. Für Lehmer wird maximal Phi(x, x^0.25) aufgerufen, da ist der Gewinn für 32-Bit-Integer kaum/nicht merkbar. Für die MPArith-Version habe ich im Gegenteil sogar (hauptsächlich für 16-Bit) die Tabelle für BaseA=4 implemetiert, d.h. BaseA=4/5 wird über $ifdef bestimmt. Hier die Zeiten (erst BaseA=4, dann BaseA=5):
Quelltext
1:
2:
3:
4:
5:
6:
7:
8:
9:
10:
11:
12:
13:
| X:\TMP>t_primpi_a4.exe
-----------------------
CPUFrequenz: 2.826 GHz
Legendre: WE mit lsumphi32 [s]: 0.409 PrimePi(2000000000) = 98222287
Meissel: WE mit lsumphi32 [s]: 0.125 PrimePi(2000000000) = 98222287
Lehmer: WE mit lsumphi32 [s]: 0.055 PrimePi(2000000000) = 98222287
X:\TMP>T_PRIMPI_A5.EXE
-----------------------
CPUFrequenz: 2.826 GHz
Legendre: WE mit lsumphi32 [s]: 0.362 PrimePi(2000000000) = 98222287
Meissel: WE mit lsumphi32 [s]: 0.099 PrimePi(2000000000) = 98222287
Lehmer: WE mit lsumphi32 [s]: 0.051 PrimePi(2000000000) = 98222287 |
Gruß Gammatester
Horst_H - Mo 23.07.12 11:55
Hallo,
absolut bringt es wirklich wenig, aber relativ bringt der Übergang von 4 auf 5 8..12% Beschleunigung.Das kann man ja für INt64/Qword im Hinterkopf behalten.
Gruß Horst
Horst_H - Mi 24.04.13 02:08
Hallo,
ich dachte es passt hier besser als im Beitrag riesieges Boolean Feld...
Also die Erweiterung von
Gammatester in größere Zahlenbereiche.
Ich hatte die Bestimmung der Abschätzung für pi(x) = x/(ln(x)-1.08...) falsch aus einem anderen Programm übernommen, weil ich dort die Anzahl bis zur Wurzel wissen wollte.
Dieses hatte ich nicht korrigiert :oops:
Durch zählen der Rauf und Runtersuche kam ich drauf, dass ich fast immer aus Regionen von 350 bis hinauf auf 115000 suchen mußte. Das hat natürlich viel Zeit gekostet.
Hier nun die neue Version 4.
Bis 1e11 in 0,9 Sekunden ist annehmbar :-)
Quelltext
1:
2:
3:
4:
| PrimePi mit MPArith V (c) 2012 Mathematiker/WE(Gammatester)
Lehmer: WE mit lsumphi32 [s]: 00:00:00.858
PrimePi(100000000000) = 4118054813 |
Als Referenz:
http://primes.utm.edu/howmany.shtml#pi_def
Gruß Horst
EDIT
1e13 wird falsch berechnet :-( ich suche noch..
Quelltext
1:
2:
| richtig: 346,065,536,839
368,028,891,082 |
Gammatester - Mi 24.04.13 10:48
| Horst_H hat folgendes geschrieben : |
1e13 wird falsch berechnet :-( ich suche noch..
Quelltext 1:
2:
| richtig: 346,065,536,839
368,028,891,082 | |
Man muss aufpassen, daß die Longint/Longword-Bereiche nicht überschritten werden. Die beiden wichtigsten sind:
- in procedure phirec: a := -x1 *ProdA+x;
- in function primepi32: z := x div Primes1E9[i];
Ich habe a und z nun int64 gemacht (was die Geschwindingkeit natürlich nicht gerade steigert). Ich habe das und einige Delphi/FPC-Inkompatibilitäten gefixt. Im Anhang die V5, sie berechnet PrimePi(10^13)
Quelltext
1:
2:
3:
4:
5:
6:
7:
8:
| X:\TMP>prime_piV5.exe
PrimePi mit Lehmer/Meissel (c) 2012-2013 Mathematiker/Gammatester/Horst_H
SiebPrim 99460729
Bis Primzahl 10007 gesucht 100000000
Anzahl Primzahlen 5761455
Lehmer: WE mit lsumphi32 [s]: 00:01:24.546
PrimePi(10000000000000) = 346065536839 |
Gruß Gammatester
Edit: Es hat mir keine Ruhe gelassen, warum a int64 sein soll (dies war der erste Fix-Versuch für 10^13). In der Tat ist an der Stelle a = x mod ProdA, also sehr klein. Habe es also wieder rückgänging gemacht und damit die Geschwindigkeit wieder etwas erhöht. Im Angang als V5a
Quelltext
1:
2:
3:
4:
5:
6:
7:
8:
| X:\TMP>prime_piV5a.exe
PrimePi mit Lehmer/Meissel (c) 2012-2013 Mathematiker/Gammatester/Horst_H
SiebPrim 99460729
Bis Primzahl 10007 gesucht 100000000
Anzahl Primzahlen 5761455
Lehmer: WE mit lsumphi32 [s]: 00:01:07.825
PrimePi(10000000000000) = 346065536839 |
Horst_H - Mi 24.04.13 12:23
Hallo,
Z hatte ich übersehen das es maximal x/ sqrt(sqrt(x)) = x^(3/4) werden kann, weshalb 10^12 funktionierte.
Aber hier
| Zitat: |
| in procedure phirec: a := -x1 *ProdA+x; |
aber dort steht insgesamt:
Delphi-Quelltext
1:
2:
3:
| x1 := x div ProdA;
a := -x1 *ProdA+x; x := x1*PhiPA; |
zudem ist a maximal sqrt(sqrt(x)) ist.
Es spart aber nur wenig Zeit statt 2min15 nun 2min2
Mist ....
Gruß Horst
Gammatester - Mi 24.04.13 12:27
| Horst_H hat folgendes geschrieben : |
Mist ....
|
Kein Problem, diesmal war ich halt ein paar ms schneller :P
Horst_H - Mi 24.04.13 14:29
Hallo,
und es ist wieder bei knapp 2 min bei mir.
Was mich doch erstaunt, dass die Abschätzung für pi(x) = x/(ln(x)-1.08..) so sehr gut funktioniert.
Im Schnitt nur 15 daneben, da versuche ich erst gar nicht den Ansatz einer linearen Suche innerhalb eines Bereiches, denn so oft, wird die Funktion nicht aufgerufen.
Quelltext
1:
2:
| 1e12 9.629.082 Aufrufe Primes1E6Index(x:LongWord):longInt; mit x > $FFFF
1e13 10.430.405 Aufrufe, während phirec 2.708.620.547 200-mal so häufig aufgerufen wird |
Da fällt mir man kann ja phirec in einer 64 und einer 32 Bit Variante probieren.
Und das bringt sogar etwas.Bei 1e13 statt 116 Sekunden nun 83 Sekunden.
Lehmer: WE mit lsumphi32 [s]: 00:01:22.971
PrimePi(10000000000000) = 346065536839
Da hat man nun 3.2 Ghz und immer noch soooooo langsam.
Gruß Horst
Gammatester - Mi 24.04.13 14:57
| Horst_H hat folgendes geschrieben : |
| Da hat man nun 3.2 Ghz und immer noch soooooo langsam. |
Bei mir (Delphi und Core 2 Duo E8300, 2.8GHz) bringt das ganze nicht besonders viel (gerade mal 6s):
Quelltext
1:
2:
3:
4:
5:
6:
7:
| PrimePi mit Lehmer/Meissel (c) 2012-2013 Mathematiker/Gammatester/Horst_H
SiebPrim 99460729
Bis Primzahl 10007 gesucht 100000000
Anzahl Primzahlen 5761455 geschaetzt 5773445
Lehmer: WE mit lsumphi32 [s]: 00:01:01.227
PrimePi(10000000000000) = 346065536839 |
Neben dem schon bekannten "normalen" FPC-Problem Dword -> Longword ist diesmal die Const-Zuweisung an Numprimes problematisch (selbst *1.0 hilft nichts), ich habe sie ersetzt durch:
Delphi-Quelltext
1:
| const NumPrimes = 5773445; |
Gruß Gammatester
Horst_H - Mi 24.04.13 15:47
Hallo,
freepascal macht seine Int64 Berechnungen scheinbar erheblich langsamer mit ernormen Geschiebe auf dem Stack.
Aber es gibt doch sicher noch andere Arten der Beschleunigung.
Ich kann mir gut vorstellen, dass es oft gleiche Kombinationen x,a für phiRec(x,a) gibt.
Vielleicht sollte man das mit einer kleinen Zahl mal testen.
Gruß Horst
Gammatester - Mi 24.04.13 16:13
| Horst_H hat folgendes geschrieben : |
| Aber es gibt doch sicher noch andere Arten der Beschleunigung. Ich kann mir gut vorstellen, dass es oft gleiche Kombinationen x,a für phiRec(x,a) gibt. |
Ja die gibt es: zB den 'richtigen' Algorithmus von
T. Oliveira e Silva [
http://www.ieeta.pt/~tos/primes.html] (und nicht nur seine von mir 'optimierte' Demo-Version); beschrieben in dem von phirec referenzierten Dokument
http://www.ieeta.pt/~tos/bib/5.4.pdf: Dort wird der Aufrufbaum für Phi(x,a) intelligent reduziert.
Aber das grenzt an Arbeit :wink: und sollte nur von wirklich Interessierten implementiert werden.
Gruß Gammatester
Horst_H - Do 25.04.13 08:09
Hallo,
es lohnt wohl nicht den Aufwand, oder erklärt das Missverhältnis Materie zu Antimaterie beim Urknall?
Gruß Horst
IhopeonlyReader - Do 25.04.13 17:32
| Delphi-Laie hat folgendes geschrieben : |
| Ob das von Delphi speichersparend umgesetzt wird, weiß ich nicht, vermute aber, daß für jeden Wert des Types "boolean" analog zu "bool" letztlich doch eine integre Variable "verschwendet" wird. Bitweiser Zugriff auf Integertypen wäre ja auch langsamer. |
Deshalb habe ich die Units uTBitBooleans geschrieben.. hiervon gibt es 2 Varinaten, 1x auf Byte basierend und einmal auf 32 bit basierend... einzeln lesen da ist die Byte-variante schneller, beim schreiben (sowohl ALLE 32/8 Bits auf einmal oder einzeln) ist die 32er Variante schneller... beim erstellen/löschen ebenfalls die 32er variante..
Alelerdings ist die ganze Sache danach um ca. 250% (von 15 auf 38 sek).. zu finden sind die Units hier:
http://www.entwickler-ecke.de/viewtopic.php?t=111430&postorder=asc&start=30 (erster post auf der seite)
Entwickler-Ecke.de based on phpBB
Copyright 2002 - 2011 by Tino Teuber, Copyright 2011 - 2024 by Christian Stelzmann Alle Rechte vorbehalten.
Alle Beiträge stammen von dritten Personen und dürfen geltendes Recht nicht verletzen.
Entwickler-Ecke und die zugehörigen Webseiten distanzieren sich ausdrücklich von Fremdinhalten jeglicher Art!