Autor |
Beitrag |
GTA-Place
Beiträge: 5248
Erhaltene Danke: 2
WIN XP, IE 7, FF 2.0
Delphi 7, Lazarus
|
Verfasst: Fr 11.11.05 22:33
Hallo.
Wir haben es geschaft eine schnelle Primzahlfunktion zu schreiben.
Nun schaffen wir das doch bestimmt auch für eine Funktion zur Berechnung von Pi.
Ich habe folgende Funktion bei Wikipedia gefunden:
Delphi-Quelltext 1: 2: 3: 4: 5: 6: 7: 8: 9: 10: 11: 12: 13: 14: 15: 16: 17:
| var R, X, Y: Integer; KTreff: Real; QTreff: Real; begin R := StrToInt(RadiusEdit.Text);
KTreff := 0; QTreff := sqr(2 * r + 1);
for X := -R to R do for Y := -R to R do if sqr(X) + sqr(Y) <= sqr(R) then KTreff := KTreff + 1;
PiLab.Caption := FloatToStr(4 * KTreff / QTreff); end; |
Quelltext 1: 2:
| Pi bei Radius 10.000: 3,141 276 394 507 36 Dauer bei Radius 10.000: 4,828 Sekunden |
(Pi = 3,141 592 653 589)
Diese Funktion ist zwar bei einem Radius von 10.000 schon näh dran,
aber bei für einen Radius von 100.000 ist sie zu langsam.
Es gibt bei Wiki noch eine Funktion:
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:
| var DX, DY: Double; Innen: Integer; Gesamt: Integer; Tropfen: Integer; begin Randomize;
Innen := 0; Tropfen := StrToInt(RadiusEdit.Text); Gesamt := Tropfen;
while (Tropfen > 0) do begin DX := 2 * Random - 1; DY := 2 * Random - 1;
if (sqr(DX) + sqr(DY)) <= 1 then inc(Innen);
dec(Tropfen); end;
PiLab.Caption := FloatToStr(4 * Innen / Gesamt); end; |
Quelltext 1: 2:
| Pi bei 100.000.000 Tropfen: 3,141 684 88 Dauer bei 100.000.000 Tropfen: 7,625 Sekunden |
Bei jedem Durchgang dieser Funktion kommt ein anderes
Ergebnis raus. Nimmt man das Mittel z.B. von 20 Durchgängen
erhalten wir für Pi -> 3,141 610 564 (Dauer: 149 Sekunden).
Gruß
GTA-Place
_________________ "Wer Ego-Shooter Killerspiele nennt, muss konsequenterweise jeden Horrorstreifen als Killerfilm bezeichnen." (Zeit.de)
|
|
BenBE
Beiträge: 8721
Erhaltene Danke: 191
Win95, Win98SE, Win2K, WinXP
D1S, D3S, D4S, D5E, D6E, D7E, D9PE, D10E, D12P, DXEP, L0.9\FPC2.0
|
Verfasst: Sa 12.11.05 13:39
Die von Dir aufgeführten Algo's sind ja wohl nun auch die "rudimentärsten", die's gibt ... Damit kommst Du nicht weit ...
Es gibt für Pi direkte Reihen, die man berechnen kann.
Ein relativ bekannter Algo ist der hier:
Quelltext 1:
| Pi := 4 * Sum(N=0..Inf, (-1)^N / (2N+1)); |
Liefert ergebnisse auf N Stellen Genauigkeit mit einer Ausführungszeit von O(10^N) ... Also nix brauchbares ...
Effizienter funktioniert da folgende Formel:
Quelltext 1:
| Pi := Sum(I=0..Inf, (1/16)*(4/(8I+1)-2/(8I+4)-1/(8I+5)-1/(8I+6))); |
die laut Dokumentation [1] in O(n) zu N Stellen Genauigkeit führt).
[1] David Bailey, Peter Borwein and Simon Plouffe: On the rapid computation of various polylogarithmic constants.
Wer das genaue Dokument brauch, aus dem ich das hab, bitte bei mir melden.
_________________ Anyone who is capable of being elected president should on no account be allowed to do the job.
Ich code EdgeMonkey - In dubio pro Setting.
|
|
Horst_H
Beiträge: 1652
Erhaltene Danke: 243
WIN10,PuppyLinux
FreePascal,Lazarus
|
Verfasst: Sa 12.11.05 14:48
Hallo,
Die Formel von David Bailey, Peter Borwein and Simon Plouffe bezieht sich auf die Berechnung der n.ten Stelle von Pi zur Basis 16:
und die waere dann nur dieser Term:
i.te Stelle in Hex = 4/(8I+1)-2/(8I+4)-1/(8I+5)-1/(8I+6)
... HexaNth hexadecimal digit of Pi e.. aus db.uwaterloo.ca/~alo...math-faq/node38.html
Dort gibt es ja einige Ansaetze.
Gruss Horst
|
|
BenBE
Beiträge: 8721
Erhaltene Danke: 191
Win95, Win98SE, Win2K, WinXP
D1S, D3S, D4S, D5E, D6E, D7E, D9PE, D10E, D12P, DXEP, L0.9\FPC2.0
|
Verfasst: Sa 12.11.05 17:53
So, hab mal eine Implementierung ausprobiert (mit meiner BigNum2-Unit). Da diese aber nativ keine Brüche kann, musst ich das anderweitig implementieren. Ist also noch nicht optimal. ATM lieg ich bei 100 Stellen (Dezimal) bei etwa 60 Sekunden. Lass ich es mir zur Basis 256 (wie die Unit intern rechnet) ausgeben, brauch er < 1s (für den ungekürzten Bruch).
Derzeitige Probleme:
- Division zu langsam
- Konvertierung von Hex^2 --> Dez lahmar schig ...
Hier mal der Source vom Haupt-Programm:
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: 76: 77: 78: 79: 80: 81: 82: 83: 84: 85: 86: 87: 88: 89: 90: 91: 92: 93: 94: 95: 96: 97: 98: 99: 100: 101: 102: 103: 104: 105: 106: 107: 108: 109: 110: 111: 112: 113:
| program CalcPi;
{$APPTYPE CONSOLE}
uses Windows, BigNum2 in 'BigNum2.pas';
var BNZ_Pi: TBigNumber; BNN_Pi: TBigNumber; N: Integer; BNZ_N: TBigNumber; BNZ_N_P2: TBigNumber; BNZ_N_P3: TBigNumber; BNZ_N_P4: TBigNumber; BNC_15: TBigNumber; BNC_16: TBigNumber; BNC_47: TBigNumber; BNC_120: TBigNumber; BNC_151: TBigNumber; BNC_194: TBigNumber; BNC_512: TBigNumber; BNC_712: TBigNumber; BNC_1024: TBigNumber;
BNZ_Loop: TBigNumber; BNN_Loop: TBigNumber; BNZ_Temp: TBigNumber; BNN_Temp: TBigNumber; begin BNZ_Pi := BMD_StrToBigNum ('0', False); BNN_Pi := BMD_StrToBigNum ('1', False);
BNC_15:= BMD_StrToBigNum('15', false); BNC_16:= BMD_StrToBigNum('16', false); BNC_47:= BMD_StrToBigNum('47', false); BNC_120:= BMD_StrToBigNum('120', false); BNC_151:= BMD_StrToBigNum('151', false); BNC_194:= BMD_StrToBigNum('194', false); BNC_512:= BMD_StrToBigNum('512', false); BNC_712:= BMD_StrToBigNum('712', false); BNC_1024:= BMD_StrToBigNum('1024', false);
SetLength(BNZ_N, 4); For N := 0 TO 1000 do Begin PDWORD(@BNZ_N[0])^ := N; BNZ_N_P2 := BM_Multiply(BNZ_N, BNZ_N); BNZ_N_P3 := BM_Multiply(BNZ_N_P2, BNZ_N); BNZ_N_P4 := BM_Multiply(BNZ_N_P2, BNZ_N_P2);
BNZ_Loop := BM_Add(BM_Add(BM_Multiply(BNC_120, BNZ_N_P2), BM_Multiply(BNC_151, BNZ_N)), BNC_47); BNN_Loop := BM_Multiply(BM_Power(BNC_16, BNZ_N), BM_Add( BM_Add( BM_Add( BM_Multiply(BNC_512, BNZ_N_P4), BM_Multiply(BNC_1024, BNZ_N_P3)), BM_Add( BM_Multiply(BNC_712, BNZ_N_P2), BM_Multiply(BNC_194, BNZ_N))), BNC_15));
BNZ_Temp := BNZ_Pi; BNN_Temp := BNN_Pi;
BNZ_Pi := BM_Add(BM_Multiply(BNZ_Temp, BNN_Loop), BM_Multiply(BNZ_Loop, BNN_Temp)); BNN_Pi := BM_Multiply(BNN_Temp, BNN_Loop);
If N mod 10 = 0 Then Begin WriteLn('Schleifendurchlauf ', N, ':'); BM_GCDOptimize(BNZ_Pi, BNN_Pi);
If N mod 100 = 0 Then Begin WriteLn(BMD_BigNumToStr(BNZ_Pi, False), '/', BMD_BigNumToStr(BNN_Pi, False)); WriteLn; end; end; end; WriteLn('Fertig!'); ReadLn; end. |
Für Pi auf 100 Stellen sagt er (gekürzt) folgendes an:
Quelltext 1: 2: 3: 4: 5: 6: 7: 8: 9: 10: 11:
| 81776705870681452759637710119896159099056300650731684809020387290127889434752159 64160027329017114926448573863946585952299480108759888420967632238816400897010401 83717597478994548857737302371546864068859147022168987554730125622577804327767124 60515312450565509670918204687509789578432142848033892744657374401530802295636765 567320584177283902506718738017979080723730358253242933069 / 26030333938181939658353705200945844812527027714565700078970330710621216818398617 34230188021090311626569877955012660007341874288310243718319029892967189797885962 22179670306638387642732741682214277675120465124326346786632721359103997871558593 74162809132074725097909353004889334815140838013032498413165583929253576194224874 228203052332856355677249981064004796964724715720540160000 |
Könnte den Bruch bitte jemand mal überprüfen (müsste kleiner als der korrekte Wert sein).
Die BigNum2-Unit gibt's bei mir auf Anfrage.
_________________ Anyone who is capable of being elected president should on no account be allowed to do the job.
Ich code EdgeMonkey - In dubio pro Setting.
|
|
uall@ogc
Beiträge: 1826
Erhaltene Danke: 11
Win 2000 & VMware
Delphi 3 Prof, Delphi 7 Prof
|
Verfasst: Sa 12.11.05 18:03
laut derive stimmts
Einloggen, um Attachments anzusehen!
_________________ wer andern eine grube gräbt hat ein grubengrabgerät
- oder einfach zu viel zeit
|
|
Horst_H
Beiträge: 1652
Erhaltene Danke: 243
WIN10,PuppyLinux
FreePascal,Lazarus
|
Verfasst: So 13.11.05 15:02
Hallo,
zur zuegigen Basisumwandlung habe ich mal etwas fuer Turbo Pascal komponiert:
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: 76: 77: 78: 79: 80: 81: 82: 83: 84: 85: 86: 87: 88: 89: 90: 91: 92: 93: 94: 95: 96: 97: 98: 99: 100: 101: 102: 103: 104: 105: 106: 107: 108: 109: 110: 111: 112: 113: 114: 115: 116: 117: 118: 119: 120: 121: 122: 123: 124: 125: 126: 127: 128: 129: 130: 131: 132: 133: 134: 135: 136: 137: 138: 139: 140: 141: 142: 143: 144: 145: 146: 147: 148: 149: 150: 151: 152: 153: 154: 155: 156: 157: 158: 159: 160: 161: 162: 163: 164: 165: 166: 167: 168: 169: 170: 171: 172: 173: 174: 175: 176: 177: 178: 179: 180: 181: 182: 183: 184: 185: 186: 187: 188: 189: 190: 191: 192: 193: 194: 195: 196: 197: 198: 199: 200: 201: 202: 203: 204: 205: 206: 207: 208: 209: 210: 211: 212: 213: 214: 215: 216: 217: 218: 219: 220: 221: 222: 223: 224: 225: 226: 227: 228: 229:
| uses crt; const STELLEN_MAX = 1024; type tFeldIndex = 0..STELLEN_MAX+1; tZiffer = byte; tBasisZiffern = array[tFeldIndex] of tZiffer; var i : tFeldINdex; Basis1, Basis2, Basis1Stelle, Basis2Stelle : tZiffer;
Basis1ZiffernFeld, Basis2ZiffernFeld :tBasisZiffern;
Procedure ADD_Feld(var Feld1,Feld2 : tBasisZiffern; var MaxIndex : tFeldINdex; Basis : tZiffer);
var Uebertrag, Sum : tZiffer; i : tFeldIndex; Begin Uebertrag := 0; for i := 0 to MaxINdex do begin Sum := Feld1[i]+Feld2[i]+UeberTrag; If Sum < Basis then Uebertrag := 0 else Begin Sum := Sum-Basis; Uebertrag := 1; end; Feld1[i] := Sum; end; If Uebertrag = 1 then If MaxIndex < High(MaxINdex) then begin MaxINdex := MaxIndex+1; Feld1[MaxINdex] := 1 end else writeln('Konvertierungsueberlauf'); end;
Procedure Doppel_Feld(var Feld :tBasisZiffern; var MaxIndex:tFeldINdex; Basis : tZiffer);
var Uebertrag, Sum : tZiffer; i : tFeldIndex; Begin Uebertrag := 0; for i := 0 to MaxINdex do begin Sum := Feld[i] shl 1 + UeberTrag; If Sum < Basis then Uebertrag := 0 else Begin Sum := Sum-Basis; Uebertrag := 1; end; Feld[i] := Sum; end; If Uebertrag = 1 then If MaxIndex < High(MaxINdex) then begin MaxINdex := MaxIndex+1; Feld[MaxINdex] := 1 end else writeln('Konvertierungsueberlauf'); end;
procedure BasKonv( Bas1,Bas2 :tZiffer; var Feld1,Feld2 :tBasisZiffern); var Summ1Feld, Summ2Feld, pFeld1, pFeld2, Tausche :^tBasisZiffern;
Feld1_Index, Feld2_Index, I, J :tFeldINdex;
Vergleich, Bas1Ziff :tZiffer; begin If (Bas1 <2) or (Bas1 <2) then begin writeln('Die Basen sind unzulaessig'); EXIT; end; if Bas1 = Bas2 then begin For Feld1_Index :=High(tFeldIndex) downto 0 do Feld2[Feld1_Index] := Feld1[Feld1_Index]; EXIT; end; For Feld1_Index :=High(tFeldIndex) downto 0 do Feld2[Feld1_Index] := 0; for Feld1_Index := High(tFeldIndex) downto 0 do begin if Feld1[Feld1_Index]>=Bas1 then begin writeln('Falsche Basis1 in Feld1'); writeln(Feld1_Index:10,Feld1[Feld1_Index]:10); EXIT; end; if Feld1[Feld1_Index]<> 0 then break; end; if Feld1[Feld1_Index] <> 0 then begin new(Summ1Feld); new(Summ2Feld); Feld2_Index := 0; fillchar(Summ1Feld^,SizeOf(tBasisZiffern),#0); Summ1Feld^[0] := 1; fillchar(Summ2Feld^,SizeOf(tBasisZiffern),#0); For j := 0 to Feld1_Index do Begin Bas1Ziff := Feld1[j]; Vergleich := 1; Tausche := Summ2feld; Summ2Feld := Summ1Feld; Summ1Feld := Tausche; fillchar(Summ1Feld^,(Feld2_Index+1)*SizeOf(tZiffer),#0); repeat if (Bas1Ziff AND Vergleich) <> 0 then ADD_Feld(Feld2,Summ2Feld^,Feld2_Index,Bas2); if (Bas1 AND Vergleich)<>0 then ADD_Feld(Summ1Feld^,Summ2Feld^,Feld2_Index,Bas2); Doppel_Feld(Summ2Feld^,Feld2_Index,Bas2); Vergleich := Vergleich +Vergleich; until VerGleich>=Bas1; if (Bas1 AND Vergleich)<>0 then Begin Tausche := Summ2feld; Summ2Feld := Summ1Feld; Summ1Feld := Tausche end; end; Dispose(Summ2feld); Dispose(Summ1feld); end; end;
Begin clrscr; Basis1 := 16; For i := High(tFeldIndex) downto 256 do Basis1ZiffernFeld[i] := 0; for i := 255 downto 0 do Basis1ZiffernFeld[i] := Basis1-1; Basis1STelle := round(ln(Basis1)/ln(10)+0.4999999); for i := High(tFeldIndex) downto 0 do if Basis1ZiffernFeld[i] <> 0 then Break; for i := i downto 0 do write(Basis1ZiffernFeld[i]:Basis1Stelle,'|'); writeln(' zur Basis ',Basis1); repeat until keypressed; while keypressed do readkey; For Basis2 := 2 to 16 do begin BasKonv(Basis1,Basis2,Basis1ZiffernFeld,Basis2ZiffernFeld); Basis2STelle := round(ln(Basis2)/ln(10)+0.4999999); for i := High(tFeldIndex) downto 0 do if Basis2ZiffernFeld[i] <> 0 then Break; for i := i downto 0 do write(Basis2ZiffernFeld[i]:Basis2Stelle,'|'); writeln(' zur Basis ',Basis2); writeln('Weiter mit ENTER '); readln; end;
END. |
Es gibt auch eine Version mit statischen Feldern, die ADD_Feld und Doppel_Feld mit Assembler loest und fuer Ziffer auch longint nutzt.
Dann kann zur Basis 1e9 konvertieren , was erheblich schneller ist, als zur Basis 10.
Z.B.: eine 426 stellige Zahl zur Basis 256 (alles $FF) wird umgewandelt:
Quelltext 1: 2: 3: 4: 5: 6: 7: 8: 9: 10: 11: 12: 13: 14: 15: 16: 17: 18: 19: 20:
| 1026 Stellen zur Basis 10 CPU-Takte: 24490394 Weiter mit Tastendruck
513 Stellen zur Basis 100 CPU-Takte: 12410510 Weiter mit Tastendruck
342 Stellen zur Basis 1000 CPU-Takte: 8455497 Weiter mit Tastendruck
257 Stellen zur Basis 10000 CPU-Takte: 6481994 Weiter mit Tastendruck
206 Stellen zur Basis 100000 CPU-Takte: 5277828 Weiter mit Tastendruck
171 Stellen zur Basis 1000000 CPU-Takte: 4482436 Weiter mit Tastendruck
147 Stellen zur Basis 10000000 CPU-Takte: 3920387 Weiter mit Tastendruck
129 Stellen zur Basis 100000000 CPU-Takte: 3666084 Weiter mit Tastendruck
114 Stellen zur Basis 1000000000 CPU-Takte: 3169290 Weiter mit Tastendruck
bei einer riesigen Zahl (50000 Dezimale) 5556 Stellen zur Basis 1000000000 CPU-Takte: 6506234219 |
Gruss Horst
Zuletzt bearbeitet von Horst_H am So 13.11.05 16:26, insgesamt 2-mal bearbeitet
|
|
BenBE
Beiträge: 8721
Erhaltene Danke: 191
Win95, Win98SE, Win2K, WinXP
D1S, D3S, D4S, D5E, D6E, D7E, D9PE, D10E, D12P, DXEP, L0.9\FPC2.0
|
Verfasst: So 13.11.05 15:12
Mein Problem mit dem aktuellen Algo liegt an der Stelle, wo er Zähler und Nenner kürzt. Bei der Ausgabe in Durchlauf 100 brauch er ~50 Sekunden allein zum Kürzen. die eigentlichen 2 Divisionen danach führt er wieder <1s aus.
Ich hab etwa folgenden Call-Path:
Quelltext 1: 2: 3:
| BM_GCDOptimize BM_GCD (Euklidischer Algo) BM_Modulo (ein Aufruf ~1-2 Sekunden) |
Im Anhang mal beide benötigten Source-Files ...
Edit: Aktualisierte BigNum-Unit angehängt. Genaue Infos gibt's weiter unten.
Einloggen, um Attachments anzusehen!
_________________ Anyone who is capable of being elected president should on no account be allowed to do the job.
Ich code EdgeMonkey - In dubio pro Setting.
Zuletzt bearbeitet von BenBE am Mi 16.11.05 20:54, insgesamt 1-mal bearbeitet
|
|
Horst_H
Beiträge: 1652
Erhaltene Danke: 243
WIN10,PuppyLinux
FreePascal,Lazarus
|
Verfasst: So 13.11.05 15:57
Hallo,
wenn ich die Zeiten stoppe
Delphi-Quelltext 1: 2: 3: 4: 5: 6: 7: 8: 9: 10: 11: 12: 13: 14: 15: 16: 17: 18: 19:
| Begin t1 := time; BM_GCDOptimize(BNZ_Pi, BNN_Pi); t1 := time-t1; WriteLn('Schleifendurchlauf ', N, ':'+FormatDateTime(' hh:mm:ss.zzz',t1)); If N Mod 100 = 0 Then Begin t1:= time; WriteLn(BMD_BigNumToStr(BNZ_Pi, False), #13#10'/'#13#10, BMD_BigNumToStr(BNN_Pi, False), #13#10'='); t1 := t1-time; writeln(FormatDateTime(' hh:mm:ss.zzz',t1)); t1:= time; WriteLn(BMD_FractionToStr(BNZ_Pi, BNN_Pi)); t1 := t1-time; writeln(FormatDateTime(' hh:mm:ss.zzz',t1)); WriteLn; readln; End; End; |
komme ich auf andere Zeiten.
Quelltext 1: 2: 3: 4: 5: 6:
| Schleifendurchlauf 70: 00:00:01.490 Schleifendurchlauf 80: 00:00:01.370 Schleifendurchlauf 90: 00:00:01.760 Schleifendurchlauf 100: 00:00:02.470 8177..../26030......... 00:00:01.650 3.14..... 00:00:01.590 |
Gruss Horst
P.S.:
Vielleicht liegt es an sysutils, das auf 8 Bytegrenzen ausrichtet??
|
|
BenBE
Beiträge: 8721
Erhaltene Danke: 191
Win95, Win98SE, Win2K, WinXP
D1S, D3S, D4S, D5E, D6E, D7E, D9PE, D10E, D12P, DXEP, L0.9\FPC2.0
|
Verfasst: So 13.11.05 16:25
Jup. Die Zeiten sind mit Kürzen der Zahlen. Ich hatte die Werte ohne Kürzen angegeben, dann käme das nämlich mit <1 Sekunde hin.
Ich hab sowieso noch eine Optimierung für BM_Modulo nachzureichen, die gerade bei großen Zahlen, die dicht beieinander liegen wesentlich mehr Performance liefern sollte.
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:
| Function BM_Modulo(Num1, Num2: TBigNumber): TBigNumber; Var X: Integer; O: Byte; Begin BM_Optimize(Num1); If (Length(Num1) < Length(Num2)) OR ((Length(Num1) - Length(Num2) < 2) AND (Length(Num1) > 256)) Then Begin BM_Assign(Result, Num1); While BM_CompareNC(Result, Num2) Do Result := BM_Sub(Result, Num2); Exit; end;
SetLength(Result, 1); FillChar(Result[0], Length(Result), 0);
For X := Length(Num1) * 8 - 1 Downto 0 Do Begin O := Integer(0 <> Num1[X Div 8] And ($01 Shl (X Mod 8)));
Result := BM_Add(Result, Result); Result[0] := Result[0] + O;
If BM_CompareNC(Result, Num2) Then Result := BM_Sub(Result, Num2); End; End; |
Ich komm (mit Optimierung von BM_Modulo) auf folgende Werte:
Quelltext 1: 2: 3: 4: 5: 6: 7: 8: 9: 10: 11: 12: 13: 14: 15: 16:
| Schleifendurchlauf 0: 00:00:00.000 Bruch: 00:00:00.000 Zahl: 00:00:00.000 Schleifendurchlauf 10: 00:00:00.031 Schleifendurchlauf 20: 00:00:00.235 Schleifendurchlauf 30: 00:00:00.640 Schleifendurchlauf 40: 00:00:01.203 Schleifendurchlauf 50: 00:00:01.219 Schleifendurchlauf 60: 00:00:01.828 Schleifendurchlauf 70: 00:00:02.641 Schleifendurchlauf 80: 00:00:03.484 Schleifendurchlauf 90: 00:00:03.547 Schleifendurchlauf 100: 00:00:05.016 Bruch: 00:00:03.312 Zahl: 00:00:02.750 Schleifendurchlauf 110: 00:00:05.078 |
Teste mal, auf welche Werte Du kommst... Nutze einen AMD XP 1600+ bei 1,4 GHz ... RAM ist uninteressant, Verbrauch liegt bei ~16MB Virtuell.
_________________ Anyone who is capable of being elected president should on no account be allowed to do the job.
Ich code EdgeMonkey - In dubio pro Setting.
|
|
Horst_H
Beiträge: 1652
Erhaltene Danke: 243
WIN10,PuppyLinux
FreePascal,Lazarus
|
Verfasst: So 13.11.05 16:38
Hallo,
mit einem 1800 Duron applebred Delphi 7 PE
Quelltext 1: 2: 3: 4: 5: 6:
| Schleifendurchlauf 70: 00:00:01.040 Schleifendurchlauf 80: 00:00:01.320 Schleifendurchlauf 90: 00:00:01.700 Schleifendurchlauf 100: 00:00:02.470 Bruch = 00:00:01.650 3.141.. 00:00:01.370 |
Ja das war wohl nichts.Fast identische Zeiten.
Gruss Horst
P.S.:
Die Rechnungen muessten mit cardinal(das geht nur mit assembler,da man den Uebertrag verliert) oder zumindest word statt Byte erheblich schneller werden.
|
|
BenBE
Beiträge: 8721
Erhaltene Danke: 191
Win95, Win98SE, Win2K, WinXP
D1S, D3S, D4S, D5E, D6E, D7E, D9PE, D10E, D12P, DXEP, L0.9\FPC2.0
|
Verfasst: So 13.11.05 17:53
Hab noch eine Optimierung gefunden ... Allerdings noch keine Möglichkeit sie wirklich zu implementieren ...
Geht darum:
Wenn ich eine Zahl X durch Y teile, shiftet der aktuelle Algo zigmal umsonst, bis er das erste mal etwas subtrahieren kann. Effizienter wäre nun, die "Dummy"-Shifts automatisiert mit einem Move auszuführen und dann nur noch die verbleibenden Shifts in der Schleife abzuhandeln. Damit würde man einen haufen Rechenzeit sparen ...
@Horst: Das die Zeiten nahezu gleich sind, liegt daran, dass die Optimierung, wie Du in der If-Bedingung sicherlich sehen wirst, erst für Zahlen > 256^256 gedacht ist, da es dort sinn macht. Inwiefern das mit der Shifting-Lösung noch benötigt wird, weiß ich nicht.
_________________ Anyone who is capable of being elected president should on no account be allowed to do the job.
Ich code EdgeMonkey - In dubio pro Setting.
|
|
Horst_H
Beiträge: 1652
Erhaltene Danke: 243
WIN10,PuppyLinux
FreePascal,Lazarus
|
Verfasst: So 13.11.05 18:31
Hallo,
ein kleiner Ansporn hier
berechnet bei mir 1 Mio Stellen in 8.52 Sekunden.
Ach du Schande, wie soll man das erreichen, zumindest in 80 Sekunden...
Gruss Horst
|
|
BenBE
Beiträge: 8721
Erhaltene Danke: 191
Win95, Win98SE, Win2K, WinXP
D1S, D3S, D4S, D5E, D6E, D7E, D9PE, D10E, D12P, DXEP, L0.9\FPC2.0
|
Verfasst: So 13.11.05 19:06
k, hab mal noch etwas optimiert:
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: 76: 77: 78: 79: 80: 81: 82: 83: 84: 85: 86: 87: 88: 89: 90: 91: 92: 93: 94: 95: 96: 97: 98: 99: 100: 101: 102: 103: 104: 105: 106: 107: 108: 109: 110: 111: 112: 113: 114: 115: 116: 117: 118: 119: 120: 121:
| Function BM_Divide(Num1, Num2: TBigNumber): TBigNumber; Var X: Integer; O: Byte; Tmp: TBigNumber; Begin BM_Optimize(Num1);
SetLength(Result, 1); FillChar(Result[0], Length(Result), 0);
If Length(Num2) > 64 Then Begin SetLength(Tmp, High(Num2)); Move(Num1[Length(Num1) - High(Num2)], Tmp[0], Length(Tmp)); X := Length(Tmp) * 8; end else Begin SetLength(Tmp, 1); FillChar(Tmp[0], Length(Tmp), 0); X := 0; end; For X := Length(Num1) * 8 - 1 - X Downto 0 Do Begin O := Integer(0 <> Num1[X Div 8] And ($01 Shl (X Mod 8)));
Tmp := BM_Add(Tmp, Tmp); Tmp[0] := Tmp[0] + O;
Result := BM_Add(Result, Result); If BM_CompareNC(Tmp, Num2) Then Begin Tmp := BM_Sub(Tmp, Num2); Inc(Result[0]); End; End; End;
Function BM_Modulo(Num1, Num2: TBigNumber): TBigNumber; Var X: Integer; O: Byte; Begin BM_Optimize(Num1); If (Length(Num1) < Length(Num2)) OR ((Length(Num1) - Length(Num2) < 2) AND (Length(Num1) > 256)) Then Begin BM_Assign(Result, Num1); While BM_CompareNC(Result, Num2) Do Result := BM_Sub(Result, Num2); Exit; end;
If Length(Num2) > 64 Then Begin SetLength(Result, High(Num2)); Move(Num1[Length(Num1) - High(Num2)], Result[0], Length(Result)); X := Length(Result) * 8; end else Begin SetLength(Result, 1); FillChar(Result[0], Length(Result), 0); X := 0; end;
For X := Length(Num1) * 8 - 1 - X Downto 0 Do Begin O := Integer(0 <> Num1[X Div 8] And ($01 Shl (X Mod 8)));
Result := BM_Add(Result, Result); Result[0] := Result[0] + O;
If BM_CompareNC(Result, Num2) Then Result := BM_Sub(Result, Num2); End; End;
Function BM_DivMod(Var Num1: TBigNumber; Num2: TBigNumber): TBigNumber; Var X: Integer; O: Byte; Tmp: TBigNumber; Begin BM_Optimize(Num1);
SetLength(Tmp, 1); FillChar(Tmp[0], Length(Tmp), 0);
If Length(Num2) > 64 Then Begin SetLength(Result, High(Num2)); Move(Num1[Length(Num1) - High(Num2)], Result[0], Length(Result)); X := Length(Result) * 8; end else Begin SetLength(Result, 1); FillChar(Result[0], Length(Result), 0); X := 0; end;
For X := Length(Num1) * 8 - 1 - X Downto 0 Do Begin O := Integer(0 <> Num1[X Div 8] And ($01 Shl (X Mod 8)));
Result := BM_Add(Result, Result); Result[0] := Result[0] + O;
Tmp := BM_Add(Tmp, Tmp); If BM_CompareNC(Result, Num2) Then Begin Result := BM_Sub(Result, Num2); Inc(Tmp[0]); End; End;
BM_Assign(Num1, Tmp); End; |
Damit ist, wenn ich mal die Anzeige der Zahlen raus lasse, die Berechnung von 1000 Stellen in 11:15 erledigt ... Vorher >1h ...
Wäre jetzt noch die Optimierung der Konvertierung zu erledigen.
@Horst: Kannst du den TP Src mal in meine Unit einbauen? Festkodiert Quellbasis 256, Zielbasis 10. Für Hex hab ich bereits nen schnellen Algo
_________________ Anyone who is capable of being elected president should on no account be allowed to do the job.
Ich code EdgeMonkey - In dubio pro Setting.
|
|
Horst_H
Beiträge: 1652
Erhaltene Danke: 243
WIN10,PuppyLinux
FreePascal,Lazarus
|
Verfasst: Mo 14.11.05 19:29
Hallo,
ein grosses Zeitproblem ist die Implementation der Division, die bitweise durchgefuehrt wird.
Das tut nicht noetig, wie man hier nachlesen kann.
Da ist auch Wichtigkeit einer moeglichst grossen Basis beschrieben.
Zitat: | Using radix 1000, instead of ten, reduces the average division time by a factor of 9.
Radix 10,000 makes multiple-length division 16 times faster than decimal division. |
Mein Basiskonverter ist auch viel zu langsam, denn er benoetigt 3.6 Sekunden fuer 50000 Stellen.
PiFast43.exe wandelt 1 Mio Stellen in 0.44 Sekunden nach Dezimal.
Da muss grundsaetzlich anders vorgegangen werden
Gruss Horst
|
|
Horst_H
Beiträge: 1652
Erhaltene Danke: 243
WIN10,PuppyLinux
FreePascal,Lazarus
|
Verfasst: Di 15.11.05 15:34
Hallo,
ich habe mal hier ein paar Pascal Programme gefunden, wovon mir dieses ganz gut gefiel, das es recht kompakt ist.
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: 76: 77: 78: 79: 80: 81: 82: 83: 84: 85: 86: 87: 88: 89: 90: 91: 92: 93: 94: 95: 96: 97: 98: 99: 100: 101: 102: 103: 104: 105: 106: 107: 108: 109: 110: 111: 112: 113: 114: 115: 116: 117: 118: 119: 120: 121: 122: 123: 124: 125: 126: 127: 128: 129: 130: 131: 132: 133: 134: 135: 136: 137: 138: 139: 140: 141: 142: 143: 144: 145: 146: 147: 148: 149: 150: 151: 152: 153: 154: 155: 156: 157: 158: 159: 160: 161: 162: 163: 164: 165: 166: 167: 168: 169: 170: 171: 172: 173: 174: 175: 176: 177: 178: 179: 180: 181: 182: 183: 184: 185: 186:
| Program ProjectPi; {$apptype Console} uses sysutils;
const STellenZahl = 10000; lnBase = 22.180709878; MaxStellen= 1038+5; type Tint = array[1..MaxStellen+1] of Cardinal;
var i,k,AktStellen :Cardinal; a,b :Tint; t : tDateTime;
procedure divi(var A:Tint;y:Cardinal);assembler; asm PushAD MOV ESI,A MOV ECX,AktStellen MOV EBX,y XOR EDX,EDX MOV EAX,Dword Ptr[ESI]
@Loop: DIV EBX MOV Dword Ptr[ESI],EAX; ADD ESI,4 DEC ECX MOV EAX,Dword Ptr[ESI] JGE @Loop;
PopAD end;
procedure mult(var A:Tint;y:Cardinal);assembler; asm PushAD Mov ECX,MaxStellen; MOV ESI,A; XOR EAX,EAX Add EAX,EAX MOV EBX,y XOR EDI,EDI LEA ESI,[4*ECX+ESI] PUSHF @Loop: XOR EDX,EDX MOV EAX,Dword Ptr[ESI] MUL EBX POPF ADC EAX,EDI; MOV EDI,EDX; MOV [ESI],EAX; PUSHF SUB ESI,4 DEC ECX JGE @Loop;
POPF PopAD end;
procedure addi(var A,B:Tint);assembler; asm PushAD; MOV EDI,A; MOV ESI,B; MOV ECX,AktStellen CLC
@Loop: MOV EAX,[ESI+4*ECX] ADC [EDI+4*ECX],EAX; DEC ECX JGE @Loop
PopAD end;
procedure negate(A:Tint);Assembler asm PushAD
MOV ESI,A MOV ECX,AktStellen
@Loop: NOT [ESI+4*ECX] DEC ECX JG @Loop
NEG [ESI+4*ECX] POPAD end;
procedure atn(var A:Tint;m:cardinal); var v,delta : double; i,k:Cardinal; begin a[1]:=1; delta :=lnBase/ln(m)/2; v:= MaxStellen*delta; k:= trunc(v); v:= v-delta; AktStellen := 5; for i:=k downto 1 do begin while i<v do begin inc(AktStellen); v := v-delta; if AktStellen = MaxStellen then v := 0; end; divi(A,2*i+1); mult(A,2*i-1); divi(A,m*m); negate(A); end; divi(A,m); end;
procedure warten; begin writeln;writeln; write(i-1,' Stellen. '); if (i<MaxStellen-4) then write('Weiter mit der Enter-taste.') else write('Zurueck zum Programm mit der Enter-taste.'); writeln; readln;
end;
Begin writeln('Berechnet werden ',STellenZahl,' Stellen von Pi.'); t:=time; atn(A,172); mult(A,88);
atn(B,239); mult(B,51); addi(A,B);
atn(B,682); mult(B,32); addi(A,B);
atn(B,5357); mult(B,44); addi(A,B);
atn(B,12943); mult(B,68); addi(A,B);
mult(A,4);
t := time-t;
writeln(A[1],'.');
i := 0; k := 50; while i <STellenZahl do begin a[1]:=0; mult(A,100000); write(Format('%.5d',[a[1]])); i := i+5; If i = k then begin writeln(Format(' %6d',[k])); k:= k+50; end; end; writeln('Es dauerte '+FormatDateTime('hh:mm:ss.zzz',t));
readln; End. |
Das braucht etwa 2min 7s (Gauss 2 min 29,2 s) fuer 100000 genaue Nachkommastellen.
Viel ist auch so nicht moeglich, da in atn Basis*(2*k-1)< Cardinal bleiben muss.
Mehr Stellen bedeutet kleinere Basis, bis irgendwann doch eine Bibliothek fuer multiple digits Mupad,GMP(2.9 Sek fuer 1 mio Stellen von pi) und wie sie alle heissen, gebraucht wird.
Eine interessante und unverstaendliche(francais) Lektuere:
membres.lycos.fr/gersoo/docs/pi/pi.pdf
Gruss Horst
P.S.:
Neue Version z.t. in Assembler.Weniger globale Variable (ausser akt(uelle)Stellen(zahl) ).Laufzeit nur noch 28 Sekunden, wobei allein uf die Vergrosserung der Basis von 2^16 auf 2^32 einen Faktor von 2 bewirkt, da man nur die Haelfte an Rechenoperationen hat.
Der Abschnitt mult mit PushF,Popf gefaellt mir garnicht...
Zuletzt bearbeitet von Horst_H am Do 17.11.05 16:16, insgesamt 2-mal bearbeitet
|
|
BenBE
Beiträge: 8721
Erhaltene Danke: 191
Win95, Win98SE, Win2K, WinXP
D1S, D3S, D4S, D5E, D6E, D7E, D9PE, D10E, D12P, DXEP, L0.9\FPC2.0
|
Verfasst: Mi 16.11.05 21:03
Ich hab hier mal eine neue Version meiner BigNum-Unit angehängt. FFT-Multiplikation wird NOCH NICHT unterstützt, wenn ich eine passende Erklärung zur Implementation finde\bekomme, kann ich das nachreichen.
Neuerungen in der Unit sind eine Multiple Length Division sowie die Partial-Multiplikation, wenn der zweite Faktor nur ein Byte lang ist.
Damit berechnet mein Programm Pi auf 1200 Stellen in etwa 2:36. (Beschleunigung um Faktor ~4,3) selbst wenn ich die Anzeige der Teilergebnisse einschalte.
Bin für Optimierungen jeder Art offen.
@Horst: Ich schau mal in deinen Konvertierungs-Algo rein. Wenn Du für 3000 dez-Stellen <5 Sekunden brauchst guck ich mal, ob sich das übernehmen lohnt.
_________________ Anyone who is capable of being elected president should on no account be allowed to do the job.
Ich code EdgeMonkey - In dubio pro Setting.
|
|
Horst_H
Beiträge: 1652
Erhaltene Danke: 243
WIN10,PuppyLinux
FreePascal,Lazarus
|
Verfasst: Do 17.11.05 08:30
Hallo,
momentan sind es also 28 Sekunden fuer 10e5 Stellen, bei immer noch quadratischer Laufzeit.
Damit halte ich das Verfahren fuer normale Div,mult ausgereizt.
Es ist damit um den Faktor 875(=100*28 / 3.2 ) von piFast entfernt.
Auch die Basisumwandlung hat quadratische Laufzeit 50000 Stellen in 3.6 s bedeutet
0.013 s fuer 3000 Stellen oder 1280 s fuer 10e6 Stellen, aber wie Du oben siehst, braucht man das nicht unbedingt.
Gruss Horst
P.S:
Kleine Aenderung, um pushF,popF bei mult sozuwerden.
Ergibt bei mir nun eine Laufzeit von ~23.6 Sekunden.
Fuer 200000 Stellen quadratische Laufzeit
Es dauerte 00:01:34.984 = 2*2 *23.6 Sekunden
ohne Divisionen sind es nur 8.688 Sekunden, da bin ich aber sehr erstaunt.
Also geht ~91% der Rechenzeit auf die Divisionen (Ahtlo64 ist auch nocht langsamer als Athlon32)
Man koennte divi(m*m) vielleicht durch die Multiplkation mit dem Kehrwert ersetzen.
Dazu braucht man die Multplikation mit grossen Zahlen (Tint*Tint);
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: 76: 77: 78: 79: 80: 81: 82: 83: 84: 85: 86: 87: 88: 89: 90: 91: 92: 93: 94: 95: 96: 97: 98: 99: 100: 101: 102: 103: 104: 105: 106: 107: 108: 109: 110: 111: 112: 113: 114: 115: 116: 117: 118: 119: 120: 121: 122: 123: 124: 125: 126: 127: 128: 129: 130: 131: 132: 133: 134: 135: 136: 137: 138: 139: 140: 141: 142: 143: 144: 145: 146: 147: 148: 149: 150: 151: 152: 153: 154: 155: 156: 157: 158: 159: 160: 161: 162: 163: 164: 165: 166: 167: 168: 169: 170: 171: 172: 173: 174: 175: 176: 177: 178: 179: 180:
| Program ProjPi; {$apptype Console} uses sysutils;
const STellenZahl = 100000; lnBase = 22.180709878; MaxStellen= 10381+5; type Tint = array[1..MaxStellen+1] of Cardinal;
var i,k,AktStellen :Cardinal; a,b :Tint; t : tDateTime;
procedure divi(var A:Tint;y:Cardinal);assembler; asm PushAD MOV ESI,A MOV ECX,AktStellen MOV EBX,y XOR EDX,EDX MOV EAX,Dword Ptr[ESI]
@Loop: DIV EBX MOV Dword Ptr[ESI],EAX; ADD ESI,4 DEC ECX MOV EAX,Dword Ptr[ESI] JGE @Loop;
PopAD end;
procedure mult(var A:Tint;y:Cardinal);assembler; asm PushAD Mov ECX,AktStellen; MOV ESI,A; XOR EAX,EAX MOV EBX,y XOR EDI,EDI LEA ESI,[4*ECX+ESI] XOR EBP,EBP @Loop: XOR EDX,EDX MOV EAX,Dword Ptr[ESI] MUL EBX ADD EAX,EBP XOR EBP,EBP ADC EBP,0; ADD EAX,EDI; MOV EDI,EDX; MOV [ESI],EAX; ADC EBP,0; SUB ESI,4 DEC ECX JGE @Loop;
PopAD end;
procedure addi(var A,B:Tint);assembler; asm PushAD; MOV EDI,A; MOV ESI,B; MOV ECX,AktStellen CLC
@Loop: MOV EAX,[ESI+4*ECX] ADC [EDI+4*ECX],EAX; DEC ECX JGE @Loop
PopAD end;
procedure negate(A:Tint);Assembler; asm PushAD
MOV ESI,A MOV ECX,AktStellen
@Loop: NOT Dword Ptr [ESI+4*ECX] DEC ECX JG @Loop
NEG Dword Ptr [ESI+4*ECX] POPAD end;
procedure atn(var A:Tint;m:cardinal); var v,delta : double; i,k:Cardinal; begin a[1]:=1; delta :=lnBase/ln(m)/2; v:= MaxStellen*delta; k:= trunc(v); v:= v-delta; AktStellen := 5; for i:=k downto 1 do begin while i<v do begin inc(AktStellen); v := v-delta; if AktStellen = MaxStellen then v := 0; end; divi(A,2*i+1); mult(A,2*i-1); divi(A,m*m); negate(A); end; divi(A,m); end;
Begin writeln('Berechnet werden ',STellenZahl,' Stellen von Pi.'); t:=time;
atn(A,172); mult(A,88);
atn(B,239); mult(B,51); addi(A,B);
atn(B,682); mult(B,32); addi(A,B);
atn(B,5357); mult(B,44); addi(A,B);
atn(B,12943); mult(B,68); addi(A,B);
mult(A,4);
t := time-t;
writeln(A[1],'.');
i := 0; k := 50; while i <STellenZahl do begin a[1]:=0; mult(A,100000); write(Format('%.5d ',[a[1]])); i := i+5; If i = k then begin writeln(Format(' %6d',[k])); k:= k+50; end; end; writeln('Es dauerte '+FormatDateTime('hh:mm:ss.zzz',t));
readln; End. |
|
|
Alice
Beiträge: 120
|
Verfasst: Do 23.08.07 12:24
hi,
ich habe (endlich) mit hilfe von delphi7 + inline asm dieses erg. hinbekommen,
als formel diente hierbei die von den gebrüdern chudnovsky die
ansich schon zu einer der schnellsten zählt.
download link (235k):
rapidshare.com/files...8021/MaxPi2.zip.html
für 1m (1048576 stellen) benötigt es hier ca. 3.162sek.
cpu intel core 2 duo @ 3700mhz
// edit: k/sek. = errechnete nachkommastellen pro sekunde in K (1024)
cu
alice
Zuletzt bearbeitet von Alice am Do 23.08.07 12:40, insgesamt 1-mal bearbeitet
|
|
Chryzler
Beiträge: 1097
Erhaltene Danke: 2
|
Verfasst: Do 23.08.07 12:38
1. Bitte lade doch dein Progrämmchen hier im Forum hoch, anstatt Rapidshare & Co.
2. Screenshots eines Fensters gehen mit Alt + Druck.
3. Respekt! In Assembler krieg ich grad noch ne Integer-Addition hin, aber den ganzen Chudnovsky-Algo zu implementieren. Nicht schlecht.
4. Afaik ist der Chudnovsky-Algo der bisher schnellste.
|
|
|