Entwickler-Ecke

Algorithmen, Optimierung und Assembler - Die Kreiszahl Pi


GTA-Place - Fr 11.11.05 22:33
Titel: Die Kreiszahl Pi
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


BenBE - 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.


Horst_H - 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 http://db.uwaterloo.ca/~alopez-o/math-faq/node38.html

Dort gibt es ja einige Ansaetze.

Gruss Horst


BenBE - 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 lahmarschig ...

Hier mal der Source vom Haupt-Programm:

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:
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; //Zähler
    BNN_Pi: TBigNumber; //Nenner

    N: Integer;
    BNZ_N: TBigNumber;  //Zählvariable als BigNumber
    BNZ_N_P2: TBigNumber;  //N^2
    BNZ_N_P3: TBigNumber;  //N^3
    BNZ_N_P4: TBigNumber;  //N^4

    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;   //Zähler des Summanden der Schleife
    BNN_Loop: TBigNumber;   //Nenner des Summanden der Schleife

    BNZ_Temp: TBigNumber;
    BNN_Temp: TBigNumber;
    
begin
    //Initialize Pi to zero
    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
        //Init the Loop data:
        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);

        //            4      2      1      1
        //Calculate ---- - ---- - ---- - ----
        //          8N+1   8N+4   8N+5   8N+6

        //          4*(8N+4)*(8N+5)*(8N+6)-2*(8N+1)*(8N+5)*(8N+6)-1*(8N+1)*(8N+4)*(8N+6)-1*(8N+1)*(8N+4)*(8N+5)
        //Calculate -------------------------------------------------------------------------------------------
        //          (8N+1)*(8N+4)*(8N+5)*(8N+6)
    
        //          960*N^2+1208*N+376
        //Calculate -------------------------------------
        //          4096*N^4+8192*N^3+5696*N^2+1552*N+120
        
        //          120*N^2+151*N+47
        //Calculate ---------------------------------
        //          512*N^4+1024*N^3+712*N^2+194*N+15

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

        //Add the values to Pi:
        //Z   A   C   A*D+C*B
        //- = - + - = -------
        //N   B   D     B*D

        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.


uall@ogc - Sa 12.11.05 18:03
laut derive stimmts


Horst_H - So 13.11.05 15:02
Hallo,

zur zuegigen Basisumwandlung habe ich mal etwas fuer Turbo Pascal komponiert:


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:
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);
{ Addiert zwei sehr lange Zahlen in einer beliebigen Basis >1}
var
  Uebertrag,
  Sum       : tZiffer;
  i         : tFeldIndex;
Begin
   Uebertrag := 0;
   {Von der kleinsten bis zur hoechsten Stellen die Ziffern addieren}
   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;
   {Falls ein Uebertrag in der hoechsten Stelle}
   {Korrektur des MaxIndex und speichern des Uebertrag's}
   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);
{ Verdoppeln einer sehr lange Zahl in einer beliebigen Basis >1}
var
  Uebertrag,
  Sum       : tZiffer;
  i         : tFeldIndex;
Begin
   Uebertrag := 0;
   {Von der kleinsten bis zur hoechsten Stellen die Ziffern addieren}
   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;
   {Falls ein Uebertrag in der hoechsten Stelle}
   {Korrektur des MaxIndex und speichern des Uebertrag's}
   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;
   {Keine Rechnung noetig? Zahl2 <- Zahl1}
   if Bas1 = Bas2 then
   begin
     For Feld1_Index :=High(tFeldIndex) downto 0 do
        Feld2[Feld1_Index] := Feld1[Feld1_Index];
     EXIT;
   end;
   {Feld2 loeschen -> Zahl ersteinmal auf Null setzen}
   For Feld1_Index :=High(tFeldIndex) downto 0 do
     Feld2[Feld1_Index] := 0;
   {Oberste Stelle<> 0 finden, dann Abbruch}
   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;
   {Nur wenn die Zahl ungleich 0 ist}
   if Feld1[Feld1_Index] <> 0 then
   begin
      {Umwandlung moeglich, erzeugen der temporaeren Summenfelder}
      new(Summ1Feld);
      new(Summ2Feld);
      Feld2_Index := 0;
      {Eine simple 1 erzeugen, fillchar ist schneller}
      {For i := 1 to STELLEN_MAX do
          Summ1Feld^[i] := 0;
      Summ1Feld^[0] := 1}
      fillchar(Summ1Feld^,SizeOf(tBasisZiffern),#0);
      Summ1Feld^[0] := 1;
      {Unbegingt Summ2Feld KOMPLETT loeschen,}
      {Da Feld2_Index oefter der geaendert wird}
      fillchar(Summ2Feld^,SizeOf(tBasisZiffern),#0);
      {Fuer alle Ziffern der Zahl zur Basis 1}
      For j := 0 to Feld1_Index do
      Begin
         Bas1Ziff := Feld1[j];
         Vergleich := 1;
         {Tausche Summ2Feld mit Summ1Feld ;-),Zeigertausch}
         Tausche   := Summ2feld;
         Summ2Feld := Summ1Feld;
         Summ1Feld := Tausche;
         {Ersetzt:}
         {For i := 0 To Feld2_Index do
         Summ2Feld^[i] := Summ1Feld^[i];}
         {Summ1feld wieder komplett loeschen}
         fillchar(Summ1Feld^,(Feld2_Index+1)*SizeOf(tZiffer),#0);
         {Ersetzt:}
         {for i := 0 to Feld2_Index do
             Summ1Feld^[i] := 0;}
         repeat
            {Zerlegung der Ziffer in ZweierPotenzen um eine}
            {schnellere Multiplikation zu simulieren}
            {*(2^n-1)  sind dann 2*n ADD + (max n Add Fuer Feld1)}
            {Also sollte Basis2 moeglichst gross gewaehlt werden}
            if (Bas1Ziff AND Vergleich) <> 0 then
               ADD_Feld(Feld2,Summ2Feld^,Feld2_Index,Bas2);
            {Dasselbe um Bas2^j zu erzeugen}
            if (Bas1 AND Vergleich)<>0 then
               ADD_Feld(Summ1Feld^,Summ2Feld^,Feld2_Index,Bas2);
            {Verdoppele SummenFeld2}
            Doppel_Feld(Summ2Feld^,Feld2_Index,Bas2);
            Vergleich := Vergleich +Vergleich;
         until VerGleich>=Bas1;
         {War Bas1 eine ZweierPotenz dann waere Summ1Feld = 0}
         if (Bas1 AND Vergleich)<>0 then
         Begin
           {For i := 0 To Feld2_Index do
               Summ1Feld^[i] := Summ2Feld^[i];}
           Tausche := Summ2feld;
           Summ2Feld := Summ1Feld;
           Summ1Feld := Tausche
         end;
         {Jetzt ist Summ1Feld = Summ1feld * Bas1}
      end;
      {Die temporaeren Felder wieder freigeben}
      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


BenBE - 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.


Horst_H - 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 - 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.


Horst_H - 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 - 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.


Horst_H - So 13.11.05 18:31
Hallo,

ein kleiner Ansporn hier [http://numbers.computation.free.fr/Constants/PiProgram/pifast.html]
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 - So 13.11.05 19:06
k, hab mal noch etwas optimiert:


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


Horst_H - 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 [http://citeseer.ist.psu.edu/cache/papers/cs/27835/http:zSzzSzwww.cs.ubc.cazSzlocalzSzreadingzSzproceedingszSzspe91-95zSzspezSz.zSzvol24zSzissue6zSzspe903.pdf/multiple-length-division-revisited.pdf] 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 - Di 15.11.05 15:34
Hallo,

ich habe mal hier [http://www.hjcaspar.de/hpxp/piprog.htm] ein paar Pascal Programme gefunden, wovon mir dieses ganz gut gefiel, das es recht kompakt ist.


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:
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;{10000 Stellen von pi mit der Stoermer-Formel}
 {$apptype Console}
 uses
   sysutils;{pi=Summe(f(i)*atn(1/m(i))) f(i),m(i) = Ganzzahlen}

 const
   STellenZahl = 10000;
   lnBase = 22.180709878;//ln(2^32)
                         //ln(10)/lnBase = 0.103810253
   MaxStellen= 1038+5;{STellenZahl*ln(10)/lnBase+5 Sicherheitsstellen hier}
 type
   //1.stelle=?[1] ist Vorkomma, der Rest Nachkommastelle
   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;             //Wie oft
     MOV ESI,A;             //Wo
     XOR EAX,EAX
     Add EAX,EAX            //clear Carry
     MOV EBX,y              //Multiplikator
     XOR EDI,EDI            //Ueberlauf des Produktes
     LEA ESI,[4*ECX+ESI]    //Genaue Adresse ohne Offset etc
     PUSHF
@Loop:
     XOR EDX,EDX            //Mul uebertrag = 0
     MOV EAX,Dword Ptr[ESI] //Multiplikand
     MUL EBX                //EAX*EBX->Produkt in EAX Ueberlauf in EDX
     POPF                   //Uebertrag der letzten Addition hervorholen
     ADC EAX,EDI;           //Ueberlauf addieren
     MOV EDI,EDX;           //Ueberlauf in EDX retten
     MOV [ESI],EAX;         //Produkt speichern
     PUSHF                  //Ueberlauf der Addition retten
     SUB ESI,4              //Adresse korrigieren
     DEC ECX                // 1 runterzaehlen
     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
 //Bestimmt 1-X
 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; {SicherheitsSTellen}
  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;
   {Stoermer}
   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:
http://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...


BenBE - Mi 16.11.05 21:03
Ich hab hier [http://www.delphi-forum.de/viewtopic.php?p=310260#310260] 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.


Horst_H - 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);

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:
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;{100000 Stellen von pi mit der Stoermer-Formel}
 {$apptype Console}
 uses
   sysutils;{pi=Summe(f(i)*atn(1/m(i))) f(i),m(i) = Ganzzahlen}

 const
   STellenZahl = 100000;
   lnBase = 22.180709878;//ln(2^32)
  //ln(10)/lnBase = 0.103810253
   MaxStellen= 10381+5;{STellenZahl*ln(10)/lnBase+5 Sicherheitsstellen hier}
 type
   //1.stelle=?[1] ist Vorkomma, der Rest Nachkommastelle
   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;    //Wie oft
     MOV ESI,A;             //Wo
     XOR EAX,EAX
     MOV EBX,y              //Multiplikator
     XOR EDI,EDI            //Ueberlauf des Produktes
     LEA ESI,[4*ECX+ESI]    //Genaue Adresse ohne Offset etc
     XOR EBP,EBP            //EBP Uebertrag
@Loop:
     XOR EDX,EDX            //Mul uebertrag = 0
     MOV EAX,Dword Ptr[ESI] //Multiplikand
     MUL EBX                //EAX*EBX->Produkt in EAX Ueberlauf in EDX
     ADD EAX,EBP            //Uebertrag der letzten Addition hervorholen
                            //Kann einen Ueberlauf erzeugen $FFF..F*$01=$FF..F
                            //ist noch nicht passiert bei 100000 Stellen
     XOR EBP,EBP            //Uebertrag vorbereiten
     ADC EBP,0;             //Uebertrag der 1. Addition retten
     ADD EAX,EDI;           //Ueberlauf addieren
     MOV EDI,EDX;           //Ueberlauf in EDX retten
     MOV [ESI],EAX;         //Produkt speichern
     ADC EBP,0;             //Uebertrag der 2. Addition retten
     SUB ESI,4              //Adresse korrigieren
     DEC ECX                // 1 runterzaehlen
     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;
 //Bestimmt 1-X
 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, gibt die Stellen die ueberhaupt jetzt berechnet werden muessen
  //zu Beginn 5 zum Ende alle ,wozu 0*x,0/x oder 0+x rechnen
  AktStellen := 5; {SicherheitsStellen}
  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;

//Hier geginnt's

Begin
   writeln('Berechnet werden ',STellenZahl,' Stellen von Pi.');
   t:=time;

   {Stoermer}
   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 - 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.

user defined image

download link (235k):

http://rapidshare.com/files/50748021/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


Chryzler - 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! :zustimm: 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.


Entwickler-Ecke.de based on phpBB
Copyright 2002 - 2011 by Tino Teuber, Copyright 2011 - 2025 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!