1(****************************************************************) 2(* 3** ARIBAS code to calculate pi to many decimal places 4** author: Otto Forster <forster@mathematik.uni-muenchen.de> 5** 6** Example call: 7** ==> pi_chud(2000). 8*) 9(*--------------------------------------------------------------*) 10(* 11** Algorithmen zur Berechnung von pi und der Eulerzahl e 12** auf viele Dezimalen (n <= 20000). 13** pi_machin(n) berechnet pi nach der Machinschen Formel 14** auf n Dezimalstellen genau; 15** pi_agm(n) benutzt zur Berechnung das arithmetisch-geometrische 16** Mittel. 17** pi_chud(n) benutzt eine Methode von Ramanujam-Chudnowski 18** 19** euler(n) berechnet e auf n Dezimalstellen. 20** Der Funktionswert ist jeweils eine ganze Zahl. Diese entspricht 21** gerundet 10**n-mal pi bzw. e. 22*) 23(*--------------------------------------------------------------*) 24(* 25** Berechnet zz*log(2) 26*) 27function log2(zz) 28var 29 x, u, k; 30begin 31 x := 0; 32 k := 0; 33 u := (zz * 2**16 * 2) div 3 34 while u /= 0 do 35 x := x + u div (2*k + 1); 36 u := u div 9; 37 inc(k); 38 end; 39 return x div 2**16; 40end. 41(*------------------------------------------------------*) 42(* 43** Berechnet zz*arctan(1/n), wird von pi_machin benutzt 44*) 45function atan1(zz,n) 46var 47 x, u, v, k, nn; 48begin 49 x := 0; 50 k := 0; 51 nn := n*n; 52 u := zz div n; 53 while u /= 0 do 54 v := u div (2*k + 1); 55 if even(k) then 56 x := x + v; 57 else 58 x := x - v; 59 end; 60 u := u div nn; 61 inc(k); 62 end; 63 return x; 64end. 65(*------------------------------------------------------*) 66(* 67** Berechnet pi * 10**n nach der Machinschen Formel 68** 69** Beispiel-Aufruf: pi_machin(1000). 70*) 71function pi_machin(n) 72var z1, x; 73begin 74 z1 := 10**n * 2**16; 75 x := 16 * atan1(z1,5) - 4 * atan1(z1,239); 76 return x div 2**16; 77end. 78(*------------------------------------------------------*) 79(* 80** Berechnet exp(1) * 10**n 81*) 82function euler(n) 83var zz, x, k; 84begin 85 zz := 10**n * 2**16; 86 x := zz * 2; 87 k := 2; 88 while zz /= 0 do 89 zz := zz div k; 90 x := x + zz; 91 inc(k); 92 end; 93 return x div 2**16; 94end. 95(*------------------------------------------------------*) 96(* 97** Berechnet pi * 10**n, 98** benutzt arithmetisch-geometrisches Mittel 99** quadratische Konvergenz 100** 101** Beispiel-Aufruf: pi_agm(1000). 102*) 103function pi_agm(n) 104var zz; 105begin 106 zz := 10**n * 2**16; 107 return piaux(zz) div 2**16; 108end. 109(*------------------------------------------------------*) 110(* 111** Hilfsfunktion fuer pi_agm 112*) 113function piaux(zz) 114var s, a, atemp, b, c, i; 115begin 116 s := 0; 117 a := zz; 118 b := isqrt(zz * (zz div 2)); 119 c := (a - b) div 2; 120 i := 1; 121 while c /= 0 do 122 writeln("eps(",i,") = ",c/zz); 123 s := s + (2**i * c * c) div zz; 124 atemp := a; 125 a := (a + b) div 2; 126 b := isqrt(atemp * b); 127 c := (a - b) div 2; 128 inc(i); 129 end; 130 return (4*a*a) div (zz - 2*s); 131end. 132(*------------------------------------------------------*) 133(* 134** Hilfsfunktion fuer pi_chud 135*) 136function Saux(zz) 137const 138 k1 = 545140134; 139 k2 = 13591409; 140 k4 = 100100025; 141 k5 = 327843840; 142var 143 A, n: integer; 144 S: integer; 145begin 146 A := zz*k1; 147 S := A * k2; 148 n := 1; 149 while A > 0 do 150 A := A * ((6*n-5)*(6*n-3)*(6*n-1)); 151 A := A div (n*n*n); 152 A := A div (k4*k5); 153 if even(n) then 154 S := S + A * (k2 + n*k1); 155 else 156 S := S - A * (k2 + n*k1); 157 end; 158 inc(n); 159 end; 160 return S div k1; 161end; 162(*--------------------------------------------------------*) 163(* 164** pi auf n Dezimalstellen nach Chudnowsky/Ramanujan 165*) 166function pi_chud(n: integer): integer; 167const 168 k3 = 640320; 169 k6 = 53360; 170var 171 zz: integer; 172 x: integer; 173begin 174 zz := 2**16 * 10**n; 175 x := isqrt(zz*zz*k3)*k6; 176 x := (zz * x) div Saux(zz); 177 return (x div 2**16); 178end; 179(*--------------------------------------------------------*) 180 181