Ein Algorithmus zur Primzahlenenumeration

Gammatester · Beiträge: 328 Erhaltene Danke: 101

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 von D.H. Lehmer aus dem Jahre 1959 (siehe auch Wiki, 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:

markieren

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

Einloggen, um Attachments anzusehen!

Horst_H · Verfasst: Mo 23.07.12 11:14

Hallo,

beeindruckend schneller!
Aber AMD dividiert nicht gern, besonders auffällig bei dem Originalcode von Mathematiker.

markieren

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.

markieren

Delphi-Quelltext

1:
2:
3:
4:
5:
6:
7:
8:
9:

const
  BaseA = 6;
  ProdA = 30030;
  PrA_2 = 15015;
  PrA_4 = 7508;
  PhiPA = 5760;

  {compressed table of Phi(x mod 30030, 6), TabPhiA[i] = lphi(2*i,5)}
  TabPhiA: array[0..PrA_4] of WORD = (...

Gruß Horst

Gammatester · Beiträge: 328 Erhaltene Danke: 101

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):

markieren

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 · Verfasst: 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 · Verfasst: Mi 24.04.13 02:08

Hallo,

ich dachte es passt hier besser als im Beitrag riesieges Boolean Feld...
Also die Erweiterung von user profile icon

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 :-)

markieren

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:
primes.utm.edu/howmany.shtml#pi_def
Gruß Horst
EDIT
1e13 wird falsch berechnet :-(

ich suche noch..

markieren

Quelltext

1:
2:

richtig: 346,065,536,839
         368,028,891,082

Einloggen, um Attachments anzusehen!

Gammatester · Beiträge: 328 Erhaltene Danke: 101

Horst_H hat folgendes geschrieben :

1e13 wird falsch berechnet :-(

ich suche noch..

markieren

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)

markieren

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

markieren

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

Einloggen, um Attachments anzusehen!

Horst_H · Verfasst: 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:

markieren

Delphi-Quelltext

1:
2:
3:

        x1 := x div ProdA;
        a := -x1 *ProdA+x; // a := X MOD ProdA < 2310;
        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 · Beiträge: 328 Erhaltene Danke: 101

Horst_H hat folgendes geschrieben :

markieren

Delphi-Quelltext

1:	a := -x1 *ProdA+x; // a := X MOD ProdA < 2310;

Mist ....

Kein Problem, diesmal war ich halt ein paar ms schneller

Horst_H · Verfasst: 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.

markieren

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

Einloggen, um Attachments anzusehen!

Gammatester · Beiträge: 328 Erhaltene Danke: 101

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):

markieren

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:

markieren

Delphi-Quelltext

1:	const NumPrimes = 5773445; {trunc(MAXPRIME/(ln(MAXPRIME)-1.1));}

Gruß Gammatester

Horst_H · Verfasst: 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 · Beiträge: 328 Erhaltene Danke: 101

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 (und nicht nur seine von mir 'optimierte' Demo-Version); beschrieben in dem von phirec referenzierten Dokument 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 · Verfasst: 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 · Beiträge: 600 Erhaltene Danke: 23 Delphi 7 PE

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:
www.entwickler-ecke....der=asc&start=30 (erster post auf der seite)

_________________
Sucht "neueres" Delphi

Wer nicht brauch was er hat, brauch auch nicht was er nicht hat!

	Mitgliederliste
	Gruppen
	Das Team
	Richtlinien
	Synonyme