Autor Beitrag
Gammatester
ontopic starontopic starontopic starontopic starontopic starontopic starontopic starhalf ontopic star
Beiträge: 328
Erhaltene Danke: 101



BeitragVerfasst: Di 18.09.12 12:35 
Das Thema ist eine Art Gegenstück zu Ein Algorithmus zur Primzahlenenumeration von user profile iconDelphi-Laie. In der neuen Version MPArith V1.22.28 habe ich eine einigermaßen schnelle Version der exakten Bestimmung der k-ten Primzahl pk implementiert (function prime32 in mp_prime.pas), die praktisch ohne Speicher auskommt. Die Idee ist folgende: Man bestimme eine Näherung p (ausgehend zB von der Wiki-Näherung),
ausblenden Delphi-Quelltext
1:
2:
3:
l1 := ln(k);
l2 := ln(l1);
p  := trunc(k*(l1+l2-1.0 + (l2-c)/l1 - ((l2-6.0)*l2 + 11.0)/(2.0*l1*l1)));
die immer kleiner ist als die genaue Zahl, berechne PrimePi(p) und zähle dann die verbleibenden Primzahlen ab. Die größte Abweichung hat meine Näherung für k=95550393. Hier ein paar Richtwerte:
ausblenden Quelltext
1:
2:
3:
4:
5:
6:
7:
8:
9:
10:
11:
D:\Xtools\MPArith>t_calc.exe
T_CALC using MPArith V1.22.28 (31/32 bit) [mp_calc]  (c) W.Ehrhardt 2006-2012
[D]:=> prime(95550393)
Result = 1942823447
Time = 84.010 ms
[D]:=> prime(105097565)
Result = 2147483647
Time = 58.589 ms
[D]:=> prime(1000000)
Result = 15485863
Time = 1.035 ms
Man sieht, daß die Korrektur doch ca 50% Rechenzeitaufschlag verbraucht. Gibt es von den üblichen Verdächtigen :wink: irgendwelche Optimierungsvorschläge oder ist der ganze Ansatz fragwürdig? Was ich vermeiden will sind riesige Tabellen.

Gruß Gammatester

Für diesen Beitrag haben gedankt: Mathematiker
Delphi-Laie
ontopic starontopic starontopic starontopic starontopic starontopic starofftopic starofftopic star
Beiträge: 1600
Erhaltene Danke: 232


Delphi 2 - RAD-Studio 10.1 Berlin
BeitragVerfasst: Di 18.09.12 13:19 
user profile iconGammatester hat folgendes geschrieben Zum zitierten Posting springen:
Was ich vermeiden will sind riesige Tabellen.


Dieses Dein Ansinnen ist berechtigt, ja löblich. Grundsätzlich sind rechenzeit- und speicheranspruchsoptimierte (bzw. -minimierte) Programmierung zwei ziemlich konträre Gegenpole. Der von mir vorgestellte (war ja nur eine Entdeckung meinerseits a.a.O.) Alorithmus benötigt allerdings keine "riesige Tabelle" sondern zwei Arrays: Eines, das alle bis dato ermittelten Primzahlen aufsteigend auflistet, ein anderes, das mit der Wurzel der bisher größten ermittelten Primzahl wächst und die Teiler beinhaltet. Dagegen benötigt das zweifellos sehr schnelle Eratosthenessieb ein ziemlich großes Array.

Die Geschwindigkeit Deines Algorithmus' hängt sicher vom Realdatentyp ab, vermutlich verwendest Du nicht real, sondern einen coprozessorunterstützen. Dennoch, dreimalige Logarithmierung und diverse Divisionen "wirken langsam". Schlau, wie genau das Ergebnis ist, werde ich aus Deinem Beitrag nicht. Reicht es, die Fließkommazahl (auf)zurunden?


Zuletzt bearbeitet von Delphi-Laie am Di 18.09.12 16:50, insgesamt 2-mal bearbeitet
Gammatester Threadstarter
ontopic starontopic starontopic starontopic starontopic starontopic starontopic starhalf ontopic star
Beiträge: 328
Erhaltene Danke: 101



BeitragVerfasst: Di 18.09.12 13:53 
user profile iconDelphi-Laie hat folgendes geschrieben Zum zitierten Posting springen:
Die Geschwindigkeit Deines Algorithmus' hängt sicher vom Realdatentyp ab, vermutlich verwendest Du nicht real, sondern einen coprozessorunterstützen. Dennoch, dreimalige Logarithmierung und diverse Divisionen "wirken langsam".
Sicher, aber die Zeit is vernachlässigbar. Richtig 'gerechnet' wird erst, wenn k>6542 ist (sonst wird die vorhandene 16-Bit-Primzahltabelle benutzt) und hier hat man zB prime(6600) = 66103 in 0.037 ms, d.h. die 3-Log-Zeit ist höchstens soviel (vermutlich aber eher im µs-Bereich).
user profile iconDelphi-Laie hat folgendes geschrieben Zum zitierten Posting springen:
Schlau, wie genau ist das Ergebnis ist, werde ich aus Deinem Beitrag nicht. Reicht es, die Fließkommazahl (auf)zurunden?
Das Ergebnis ist immer exakt, die Nährungsformel ist asymptotisch (relativ), die absoluten Abweichungen werden immer größer! Runden reicht da natürlich nicht! Wie gesagt, die größte Abweichung im Longintbereich ist für k=95550393: Der Näherungswert ist hier 1942781249, der richtige Wert 1942823447 (zB via Wolfram Alpha), die Differenz von 42198 wird mit FindNextPrime32-Aufrufen überbrückt.
Delphi-Laie
ontopic starontopic starontopic starontopic starontopic starontopic starofftopic starofftopic star
Beiträge: 1600
Erhaltene Danke: 232


Delphi 2 - RAD-Studio 10.1 Berlin
BeitragVerfasst: Di 18.09.12 14:18 
"Das Ergebnis ist immer exakt" <-> "die Nährungsformel ist asymptotisch (relativ), die absoluten Abweichungen werden immer größer! Runden reicht da natürlich nicht!", tut mir leid, damit kann ich nichts anfangen, das widerspricht sich m.E. diametral. Wahrscheinlich verstehe ich irgendetwas an Deinen Aussagen nicht.

Mit dem obigen Term (p:=...) komme ich auch nicht klar. Dort steht ein c, wo im Original eine 2 erscheint. Aber auch sonst stimmt Dein Term nicht mit Original überein: Es fehlt die Quadrierung von ln ln n (bzw. l2), und mit dem letzten (rechtesten) Term des Originals, der bei Dir gar nicht erscheint, kann ich auch nichts anfangen, weil dort eine Variable o wie aus dem Nichts plötzlich auftaucht.

Kannst Du nicht einfach ein ganz einfaches Programm hier veröffentlichen, daß z.B. die Enumeration oder wahlweise die Generierung der n. (bzw. k.) Primzahl auf Grundlage der Formel demonstriert, bitte? Ich war damit bisher jedenfalls nicht erfolgreich.
Horst_H
ontopic starontopic starontopic starontopic starontopic starontopic starontopic starofftopic star
Beiträge: 1652
Erhaltene Danke: 243

WIN10,PuppyLinux
FreePascal,Lazarus
BeitragVerfasst: Di 18.09.12 14:36 
Hallo,

der Ablauf, der Suche nach der k.ten Primzahl ist jetzt doch die Bestimmung einer Zahl pl , mittels der WIKI-Formel, die garantiert noch nicht die k.te ist und mit der Lehmer/? Algorithmus wird die exakte Primzahlnummer l von pl bestimmt und anschliessend die nächsten k-l Primzahlen bestimmt um zur k.ten zu kommen.

Es geht doch darum abschnittsweise ein "einfache" Näherungsgleichung zu finden, die näher an der Kurve, aber garantiert unterhalb ist, damit die Suchrichtung nur aufwärts ist.

Wie bestimmt man also günstig neue Konstanten 1.0/6.0/11.0 die es besser machen?
Als erstes fällt mir konvexe Hülle ein. Aber ich glaube nicht, damit die Anzahl stark der Punkte erheblich reduzieren zu können.

Gruß Horst
Gammatester Threadstarter
ontopic starontopic starontopic starontopic starontopic starontopic starontopic starhalf ontopic star
Beiträge: 328
Erhaltene Danke: 101



BeitragVerfasst: Di 18.09.12 14:51 
user profile iconDelphi-Laie hat folgendes geschrieben Zum zitierten Posting springen:
"Das Ergebnis ist immer exakt" <-> "die Nährungsformel ist asymptotisch (relativ), die absoluten Abweichungen werden immer größer! Runden reicht da natürlich nicht!", tut mir leid, damit kann ich nichts anfangen, das widerspricht sich m.E. diametral.
Nein überhaupt nicht, das ist doch gerade der Sinn des ganzen, die Näherung ist nicht genau, aber damit liegen man nahe genug am richtigen Ergebnis, analog ist zB die iterative Berechung der Quadratwurzel mit dem Newtonverfahren.

user profile iconDelphi-Laie hat folgendes geschrieben Zum zitierten Posting springen:
Mit dem obigen Term (p:=...) komme ich auch nicht klar. Dort steht ein c, wo im Original eine 2 erscheint.
Richtig, ich habe aber nie gesagt, daß ich genau die Formel benutze sondern "ausgehend zB von ..."

user profile iconDelphi-Laie hat folgendes geschrieben Zum zitierten Posting springen:
Aber auch sonst stimmt Dein Term nicht mit Original überein: Es fehlt die Quadrierung von ln ln n (bzw. l2), und mit dem letzten (rechtesten) Term des Originals, der bei Dir gar nicht erscheint, kann ich auch nichts anfangen, weil dort eine Variable o wie aus dem Nichts plötzlich auftaucht.
Also in meinem Verständnis ist l2*l2 = (ln(ln(k))^2 :?: Außerdem ist das o keine Variable, sondern das allseits bekannte Landau-Symbol!

user profile iconDelphi-Laie hat folgendes geschrieben Zum zitierten Posting springen:
Kannst Du nicht einfach ein ganz einfaches Programm hier veröffentlichen, daß z.B. die Enumeration oder wahlweise die Generierung der n. (bzw. k.) Primzahl auf Grundlage der Formel demonstriert, bitte? Ich war damit bisher jedenfalls nicht erfolgreich.
Ich habe zwar den Link angegeben, aber hier noch mal die Funktion im Quelltext
ausblenden volle Höhe 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:
{---------------------------------------------------------------------------}
function prime32(k: longint): longint;
  {-Return the kth prime if 1 <= k <= 105097565, 0 otherwise}
var
  i,pk: longint;
  c,l1,l2: extended;
  ctx: TPrimeContext;
const
  c1 = 2.006218;
  c2 = 2.004383;
  k1 = 57000000;
begin
  if (k<1or (k>105097565then prime32 := 0
  else if k<=6542 then prime32 := Primes16[k]
  else begin
    if k<k1 then c := c1
    else c := c2;
    l1 := ln(k);
    l2 := ln(l1);
    {Get approximation for pk, compute primepi, and correct}
    pk := trunc(k*(l1+l2-1.0 + (l2-c)/l1 - ((l2-6.0)*l2 + 11.0)/(2.0*l1*l1)));
    {This is the sligtly modified asymptotic expansion for pk (where c = 2), }
    {see e.g. http://functions.wolfram.com/13.03.06.0004.01 or P. Dusart, The}
    {kth prime is greater than k(log k + log log k - 1) for k>=2, Math. Comp.}
    {68, p 411, 1999. The computed pk is not greater than the real pk for all}
    {pk <= MaxLongint; verified by exhaustive computation with the procedure }
    {check_pk_formula4. The maximum deviation of 42198 occurs for k=95550393,}
    {this gives an i value from primepi32 which is 2314 to low.}
    FindFirstPrime32(pk,ctx);
    i := primepi32(ctx.prime);
    while i<k do begin
      FindNextPrime32(ctx);
      inc(i);
    end;
    prime32 := ctx.prime;
  end;
end;
Gammatester Threadstarter
ontopic starontopic starontopic starontopic starontopic starontopic starontopic starhalf ontopic star
Beiträge: 328
Erhaltene Danke: 101



BeitragVerfasst: Di 18.09.12 15:28 
user profile iconHorst_H hat folgendes geschrieben Zum zitierten Posting springen:
Wie bestimmt man also günstig neue Konstanten 1.0/6.0/11.0 die es besser machen? Als erstes fällt mir konvexe Hülle ein. Aber ich glaube nicht, damit die Anzahl stark der Punkte erheblich reduzieren zu können.
Ein Problem bei der Sache ist folgendes: Der Fehlerterm für pk/k ist zwar o(1/ln(ln(k))^2), aber für den Maximalwert k=105097565 hat man immerhin einen geschätzten absoluten Fehlerrest der Größenordnung k/ln(ln(k))^2 = 12358539.42, und im Vergleich dazu scheint mir (mit den gefunden c-Werten) die Abweichung von 42198 relativ klein. Ich habe auch schon Approximationen mit mehr Termen der Entwicklung verwendet, allerdings habe ich keinen besseren Wert als die 42198 erricht. Jedoch könnte man wahrscheinlich die mittlere Abweichung verkleinern.

Danke für den Hinweis konvexe Hülle.

Gruß Gammatester
Mathematiker
ontopic starontopic starontopic starontopic starontopic starontopic starhalf ontopic starofftopic star
Beiträge: 2622
Erhaltene Danke: 1447

Win 7, 8.1, 10
Delphi 5, 7, 10.1
BeitragVerfasst: Di 18.09.12 19:43 
Hallo Gammatester,
user profile iconGammatester hat folgendes geschrieben Zum zitierten Posting springen:
Man sieht, daß die Korrektur doch ca 50% Rechenzeitaufschlag verbraucht. Gibt es von den üblichen Verdächtigen :wink: irgendwelche Optimierungsvorschläge oder ist der ganze Ansatz fragwürdig? Was ich vermeiden will sind riesige Tabellen.

Die Routine ist wirklich sehr schnell. Sehr schön.
Optimierungsvorschläge sehe ich im Moment in prime32 eigentlich nicht. Das ist so optimal.
Da Du nur eine maximale Abweichung = 42198 von der gesuchten Primzahl hast, könnte man es ja auch mit einem "kleinen" Sieb versuchen. Aber Du hast ja geschrieben, dass Du möglichst kaum Speicher belegen möchtest. Damit bringt ein Sieb nichts, wahrscheinlich auch zeitlich nichts.

Ich denke, wenn überhaupt, sind nur Optimierungen in FindNextPrime32 zu suchen. Allerdings gebe ich zu, dass ich trotz Deiner eingefügten Kommentare noch nicht ganz verstanden habe, was bei
ausblenden Delphi-Quelltext
1:
2:
3:
4:
5:
6:
7:
      {move n to next prime residue class mod MP_MOD and index id into diff array}
      for k:=modnpd2(next) to (NPRC_MOD div 2)-1 do begin
        index := NPRC_OddIdx[k];
        {note: loop terminates via break because NPRC_OddIdx[(NPRC_MOD div 2)-1]<>-1}
        if index<>-1 then break;
        inc(next,2);
      end;

geschieht.
Insbesondere den Hintergrund von NPRC_OddIdx und NPRC_Diff verstehe ich noch nicht.

Beste Grüße
Mathematiker

_________________
Töten im Krieg ist nach meiner Auffassung um nichts besser als gewöhnlicher Mord. Albert Einstein
Delphi-Laie
ontopic starontopic starontopic starontopic starontopic starontopic starofftopic starofftopic star
Beiträge: 1600
Erhaltene Danke: 232


Delphi 2 - RAD-Studio 10.1 Berlin
BeitragVerfasst: Di 18.09.12 23:56 
Bei der Enumeration ist der vorgestellte Algorithmus tatsächlich schneller als der von mir publizierte, wie ich nunmehr herausfand.

Jedoch weiß ich leider nicht, inwieweit die 16-Bit-Primzahltabelle dafür ausschlaggebend ist. Vermutlich stark, da bis dahin ja nicht gerechnet wird (wie ich eingangs schon schrieb, entweder man spart beanspruchten Speicherplatz oder Rechenzeit). Leider sehe ich mich nicht imstande, diese Tabelle zu deaktivieren bzw. zu löschen, um diesen hier vorgestellten Algorithmus von Beginn an zum Rechnen zu zwingen. Mit eben diesen ungleichen Waffen ist kein gerechter Vergleich möglich. Mindestens ebenso interessant, wenn nicht sogar interessanter als die Geschwindigkeit ist aber der dahintersteckende Algorithmus, den ich derzeit auch nicht durchschaue.

Ergänzung: Eine Benutzung des Primes16-Arrays ist weder in der FindFirstPrime32- noch in der FindNextPrime32-Prozedur erkennbar.
Gammatester Threadstarter
ontopic starontopic starontopic starontopic starontopic starontopic starontopic starhalf ontopic star
Beiträge: 328
Erhaltene Danke: 101



BeitragVerfasst: Mi 19.09.12 09:36 
user profile iconMathematiker hat folgendes geschrieben Zum zitierten Posting springen:
Insbesondere den Hintergrund von NPRC_OddIdx und NPRC_Diff verstehe ich noch nicht.
Bei Nextprime32 und FindnextPrime32 werden nur geeignete Kandidaten mod 210 getestet (bekannter ist das für mod 30, siehe zB Wiki mod 30). Mit modnpd2(next) wird praktisch n mod 105 gerechnet (ist hier als Funktion, weil in alten Versionen noch die Mod30-Tabellen optional verwendet werden konnten). Dann wird in der Tabelle NPRC_OddIdx der Schrittweiten-Index des ersten möglichen Kandidaten gesucht. Ab da wird dann das Schrittweitenarray zyklisch durchlaufen, also im Prinzip n := n + Diff[i]; i := i+1 mod 48; wobei 48 = EulerPhi(210).

Nextprime32 und FindnextPrime32 unterscheiden sich nur dadurch, daß bei FindnextPrime32 die Initialiliserung nicht gemacht wird (die ist ausgelagert in FindFirstPrime32, der Index in die Schrittweitentabelle und der nächste Kandidat sind im Kontextrecord).

user profile iconDelphi-Laie hat folgendes geschrieben Zum zitierten Posting springen:
Bei der Enumeration ist der vorgestellte Algorithmus tatsächlich schneller als der von mir publizierte, wie ich nunmehr herausfand.

Jedoch weiß ich leider nicht, inwieweit die 16-Bit-Primzahltabelle dafür ausschlaggebend ist. Vermutlich stark, da bis dahin ja nicht gerechnet wird
Nein überhaupt nicht, denn der Bereich ist im Gegenteil ja praktisch = 0, genauer 6542/105097565 = 0.000062.

user profile iconDelphi-Laie hat folgendes geschrieben Zum zitierten Posting springen:
Ergänzung: Eine Benutzung des Primes16-Arrays ist weder in der FindFirstPrime32- noch in der FindNextPrime32-Prozedur erkennbar.
Das ist ja auch der Sinn des ganzen: Wenn man die Primzahleigenschaft einer 32-Bitzahl n durch Primzahl-Testdivisionen bis sqrt(n) bestimmt, ist das zwar theoretisch korrekt aber mit Sicherheit nicht schnell.

Gruß Gammatester

Für diesen Beitrag haben gedankt: Delphi-Laie
Delphi-Laie
ontopic starontopic starontopic starontopic starontopic starontopic starofftopic starofftopic star
Beiträge: 1600
Erhaltene Danke: 232


Delphi 2 - RAD-Studio 10.1 Berlin
BeitragVerfasst: Mi 19.09.12 15:09 
Gut, danke, dann muß ich wohl auch die obige Aussage zum Primes16-Array falsch verstanden haben.

Nach längerer Mühe, gelang es mir, den Algorithmus soweit zu extrahieren, daß ich auf die MPA-Units verzichten konnte. Funktioniert jetzt soweit, wird compiliert und läuft. Jedoch meldet sich der Compiler (Delphi 4) mit einer Meldung, die ich nicht mal als Fehlermeldung zu bezeichnen wage (Anhang). Sowas erlebte ich noch nie. Ist er wohl überfordert. Besonders drollig ist, daß ich NPRC_MOD schon mit ihrem Wert 210 im Quelltext ersetzte. Neuigkeit: Nach dem Neustart der IDE ist diese Meldung verschwunden.

Dieser Algorithmus ist es sicher wert, in separierter Form in die Bibliothek einzufließen, jedoch ist er selbst in dieser extrahierten Form für mich leider nicht nachvollziehbar, und er ist eben leider auch nicht kurz. Es spielen einfach zu viele Unterprogramme und Datenstrukturen mit. Beim Bemühen um ein Verständnis helfen mir auch die Eingangserläuterungen nicht, es ist zum "Mäusemelken"! *seufz*
Einloggen, um Attachments anzusehen!
Gammatester Threadstarter
ontopic starontopic starontopic starontopic starontopic starontopic starontopic starhalf ontopic star
Beiträge: 328
Erhaltene Danke: 101



BeitragVerfasst: Mi 19.09.12 15:46 
user profile iconDelphi-Laie hat folgendes geschrieben Zum zitierten Posting springen:
Beim Bemühen um ein Verständnis helfen mir auch die Eingangserläuterungen nicht, es ist zum "Mäusemelken"! *seufz*
Die Grundidee ist folgende: Wie den letzen Beiträgen Deines Threads gezeigt wurde, kann man die Primzahl-Zählfunktion pi(x), d.h. die Anzahl der Primzahlen <= x berechnen, ohne wirklich alle Primzahlen zu zählen. Das wird hier ausgenutzt:

Es gilt doch pi(p_k) = k. Wenn q eine Primzahl-Näherung für p_k ist mit q <= p_k, dann berechnen wir i = pi(q) <= k. Wenn i=k ist, sind wir fertig, ansonsten bestimmen wir die nächsten k-i Primzahlen.

Gruß Gammatester

Edit-PS: q wird bestimmt als nächste Primzahl größergleich der modifizierten Wiki-Näherung.
Horst_H
ontopic starontopic starontopic starontopic starontopic starontopic starontopic starofftopic star
Beiträge: 1652
Erhaltene Danke: 243

WIN10,PuppyLinux
FreePascal,Lazarus
BeitragVerfasst: Mi 19.09.12 16:57 
Hallo,

der schnelle Primzahltest für DWord wurde zum Beispiel auch von Hagen Reddmann , vor Jahren ( 2000) in
www.swissdelphicente.../showcode.php?id=914 genutzt.

@gammatester, was stört dich an 85 ms? für 42198 Tests. ~ 5e5 Tests/Sekunden

EulerPhi(210=2*3*5*7)ist nicht 48, sondern ((2-1)*(3-1)*(5-1)*(7-1)) = 48= Anzahl der im Intervall 210 zu 2,3,5,7 teilerfremden Zahlen.

Es sind doch alle Primzahlen<65536 schon vorhanden, dann kann man auch sieben ( Ach, schooon wieder..)
Mein segmentiertes Sieb mit 30030 ungeraden Zahlen ( 1..60061) , die teilerfremd zu 2,5,7,11,13 sind, siebt und zählt bis 1e9 in 3,104 Sekunden und bis 1,01e9 ( 1e7 Mio Zahlen mehr) in 32 ms mehr.
Maximal ist die Formel ja nur 42198*210/48 ~185000 davon entfernt.Das sind dann drei Siebungen statt 167 für 1e7.
Zur Initialisierung müßte man aber für maximal 6542 Primzahlen die Position im erstem Sieb (Startzahl = (Zahl DIV 60060) *60060+1) berechnen, bei meinem Rechner etwa 0,2 ms

Das meinte ich übrigens auch für FirstPrime NextPrime allgemein, wenn man viele aufeinanderfolgende Zahlen testen will.
Mit FirstPrime das Sieb initialisieren und mit NextPrime einfach das Sieb beackern.Das macht eben auch bei sehr großen Zahlen Sinn, wenn man schon mit allen Prime16 getestet hat, bevor man MillerRabin drauf loslässt.
Uups, mit 60060 komme ich nur bis 3,607e9 und nicht 4,295e9, aber 30030 Byte passen so schön in den Level-I cache und es war passend für Longint.
Ich probiere das mal umzusetzen.

Gruß Horst
Mathematiker
ontopic starontopic starontopic starontopic starontopic starontopic starhalf ontopic starofftopic star
Beiträge: 2622
Erhaltene Danke: 1447

Win 7, 8.1, 10
Delphi 5, 7, 10.1
BeitragVerfasst: Mi 19.09.12 17:10 
Hallo Horst_H,
user profile iconHorst_H hat folgendes geschrieben Zum zitierten Posting springen:
der schnelle Primzahltest für DWord wurde zum Beispiel auch von Hagen Reddmann , vor Jahren ( 2000) in www.swissdelphicente.../showcode.php?id=914 genutzt.

Dieser Text wird oft genannt, nur leider funktioniert er nicht, da ein Sprung nach @@3 erfolgen soll, die Marke @@3 aber nicht existiert.

Beste Grüße
Mathematiker

_________________
Töten im Krieg ist nach meiner Auffassung um nichts besser als gewöhnlicher Mord. Albert Einstein
Horst_H
ontopic starontopic starontopic starontopic starontopic starontopic starontopic starofftopic star
Beiträge: 1652
Erhaltene Danke: 243

WIN10,PuppyLinux
FreePascal,Lazarus
BeitragVerfasst: Mi 19.09.12 18:22 
Hallo,

das läßt sich sicher rausfinden.
Es geht ja nach primes.utm.edu/prove/prove2_3.html etwa in der Mitte.
Zitat:
If n < 9,080,191 is a both 31 and 73-SPRP, then n is prime.

Dann kann man ja erahnen wo @@3 liegen muss, aber nicht heute ...

Gruß Horst
Delphi-Laie
ontopic starontopic starontopic starontopic starontopic starontopic starofftopic starofftopic star
Beiträge: 1600
Erhaltene Danke: 232


Delphi 2 - RAD-Studio 10.1 Berlin
BeitragVerfasst: Mi 19.09.12 22:12 
user profile iconMathematiker hat folgendes geschrieben Zum zitierten Posting springen:
Hallo Horst_H,
user profile iconHorst_H hat folgendes geschrieben Zum zitierten Posting springen:
der schnelle Primzahltest für DWord wurde zum Beispiel auch von Hagen Reddmann , vor Jahren ( 2000) in www.swissdelphicente.../showcode.php?id=914 genutzt.

Dieser Text wird oft genannt, nur leider funktioniert er nicht, da ein Sprung nach @@3 erfolgen soll, die Marke @@3 aber nicht existiert.


Auch der zweite Algorithmus auf jener Seite ist nach meiner Erfahrung fehlerhaft: Der Startdivisor vor der Schleife muß mit 5 beginnen. 2 und 3 müssen dann als Quotienten explizit vornangestellt werden. Auch wenn ich von Hagen sonst bisher qualitativ hochwertige Quelltexte wahrnahm: Die Bewertung dieser beiden fehlerhaften Algorithmen im mittleren Bereich kommt mir fast schon manipuliert vor.
Gammatester Threadstarter
ontopic starontopic starontopic starontopic starontopic starontopic starontopic starhalf ontopic star
Beiträge: 328
Erhaltene Danke: 101



BeitragVerfasst: Do 20.09.12 09:40 
user profile iconHorst_H hat folgendes geschrieben Zum zitierten Posting springen:
@gammatester, was stört dich an 85 ms?
Mich stören ja eigentlich nicht die 85ms als solche, sondern der Umstand, daß durch die sub-optimale Approximation 35ms verbraten werden, primepi(1942823447) wird in 49.9 ms berechnet.

user profile iconHorst_H hat folgendes geschrieben Zum zitierten Posting springen:
EulerPhi(210=2*3*5*7)ist nicht 48, sondern ((2-1)*(3-1)*(5-1)*(7-1)) = 48= Anzahl der im Intervall 210 zu 2,3,5,7 teilerfremden Zahlen.
:?: Also für mich (und Pari, Maple etc) ist EulerPhi(210)=48. Und für Dich doch irgendwie auch, also was willst Du damit sagen?

user profile iconHorst_H hat folgendes geschrieben Zum zitierten Posting springen:
Es sind doch alle Primzahlen<65536 schon vorhanden, dann kann man auch sieben ( Ach, schooon wieder..) Mein segmentiertes Sieb mit 30030 ungeraden Zahlen ( 1..60061) , die teilerfremd zu 2,5,7,11,13 sind, siebt und zählt bis 1e9 in 3,104 Sekunden und bis 1,01e9 ( 1e7 Mio Zahlen mehr) in 32 ms mehr.
Das ist auch kein Problem, wie dem Sieb in MPArith erhält man Time [s]: 2.450, 50847534-te Primzahl = 999999937
Gammatester Threadstarter
ontopic starontopic starontopic starontopic starontopic starontopic starontopic starhalf ontopic star
Beiträge: 328
Erhaltene Danke: 101



BeitragVerfasst: Do 20.09.12 12:57 
Inzwischen habe ich die Anregung von Horst_H aufgegriffen und in der prime32-Funktion das MPArith-Sieb verwendet, hier der geänderte Pascalcode:
ausblenden 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:
function nprime32(k: longint): longint;
  {-Return the kth prime if 1 <= k <= 105097565, 0 otherwise}
var
  i,pk: longint;
  c,l1,l2: extended;
  sieve: TSieve;
const
  c1 = 2.006218;
  c2 = 2.004383;
  k1 = 57000000;
begin
  if (k<1or (k>105097565then nprime32 := 0
  else if k<=6542 then nprime32 := Primes16[k]
  else begin
    if k<k1 then c := c1
    else c := c2;
    l1 := ln(k);
    l2 := ln(l1);
    pk := trunc(k*(l1+l2-1.0 + (l2-c)/l1 - ((l2-6.0)*l2 + 11.0)/(2.0*l1*l1)));
    prime_sieve_init(sieve,pk);
    pk := prime_sieve_next(sieve);
    i  := primepi32(pk);
    while i<k do begin
      pk := prime_sieve_next(sieve);
      inc(i);
    end;
    prime_sieve_clear(sieve);
    nprime32 := pk;
  end;
end;
Das Ergebnis war doch ziemlich überraschend für mich! Es ist mir jetzt fast peinlich, aber ich habe mich bei den 85ms auf den frisch gestarteten T_Calc verlassen :oops: Offensichtlich kann man mit dem Warmlaufen des Codes, Cachefeinheiten etc. nicht vorsichtig genug sein! Mit einen Testprogramm (Pascal und Delphi6-Exe im Anhang), das (hoffentlich) erst Prozessor und Cache usw anwärmt und dann jeweils hundertmal rechnet, erhalte ich jetzt (New: mit Sieb, Old: mit findnextprime32):
ausblenden Quelltext
1:
2:
3:
4:
5:
6:
7:
8:
C:\TEST>t_primek.exe
Test of MP library version 1.22.28 (31/32 bit)   (c) W.Ehrhardt 2012
New [s]: 0.052  k= 95550393  prime(k)=1942823447
Old [s]: 0.058  k= 95550393  prime(k)=1942823447
New [s]: 0.057  k=105097565  prime(k)=2147483647
Old [s]: 0.061  k=105097565  prime(k)=2147483647
New [s]: 0.005  k= 10000000  prime(k)=179424673
Old [s]: 0.005  k= 10000000  prime(k)=179424673

Also ist der neue Code mit Sieb bei großen Zahlen etwa 10% schneller, aber insgesamt kann man wohl auch mit der alten Version leben.

Nochmal vielen Dank an alle für die Anregungen und Fragen.

Gruß Gammatester
Einloggen, um Attachments anzusehen!
Horst_H
ontopic starontopic starontopic starontopic starontopic starontopic starontopic starofftopic star
Beiträge: 1652
Erhaltene Danke: 243

WIN10,PuppyLinux
FreePascal,Lazarus
BeitragVerfasst: Fr 21.09.12 08:28 
Hallo,

OffTopic:
"Wahn!Wahn! Überall Wahn" ( Wagner, Meistersinger),
Oh shit, ich schaue in den Spiegel ;-)

OnTopic:
Ich war zu sehr im Zählen der Primzahlen und es sind mehr als 48 Primzahlen bis 210...
Ich habe doch tataächlich Probleme, den richtige Übertrag in das segmentierte Sieb zu berechnen.
Die Startzahl des Siebes bestimme ich zu:
StartZahl := (TestZahl DiV Sieblaenge)*Sieblaenge;
Den Uebertrag zu
Uebertrag := Primzahl -(Startzahl MOD PrimZahl);
Wenn der Rest Startzahl MOD PrimZahl = 11 ist, dann ist doch der Rest in das Sieb hinein = Primzahl- 11 oder nicht?

Gruß Horst
Gammatester Threadstarter
ontopic starontopic starontopic starontopic starontopic starontopic starontopic starhalf ontopic star
Beiträge: 328
Erhaltene Danke: 101



BeitragVerfasst: Fr 21.09.12 09:30 
user profile iconHorst_H hat folgendes geschrieben Zum zitierten Posting springen:

Ich war zu sehr im Zählen der Primzahlen und es sind mehr als 48 Primzahlen bis 210...
Lass gut sein so früh am morgen: primepi(210) = 46 :wink:

user profile iconHorst_H hat folgendes geschrieben Zum zitierten Posting springen:
Ich habe doch tataächlich Probleme, den richtige Übertrag in das segmentierte Sieb zu berechnen.
Die Startzahl des Siebes bestimme ich zu:
StartZahl := (TestZahl DiV Sieblaenge)*Sieblaenge;
Den Uebertrag zu
Uebertrag := Primzahl -(Startzahl MOD PrimZahl);
Wenn der Rest Startzahl MOD PrimZahl = 11 ist, dann ist doch der Rest in das Sieb hinein = Primzahl- 11 oder nicht?
Ist auch meine Auffassung, wenn Dein Sieb ein array[0..] ist und der Rest nicht 0 ist, allgemein also etwa so:
ausblenden Delphi-Quelltext
1:
2:
Uebertrag := (Startzahl MOD PrimZahl);
if Uebertrag<>0 then Uebertrag := Primzahl - Uebertrag;

Gruß Gammatester