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