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
|
(****************************************************************)
(*
** ARIBAS code to calculate pi to many decimal places
** author: Otto Forster <forster@rz.mathematik.uni-muenchen.de>
**
** Example call:
** ==> pi_chud(2000).
*)
(*--------------------------------------------------------------*)
(*
** Algorithmen zur Berechnung von pi und der Eulerzahl e
** auf viele Dezimalen (n <= 20000).
** pi_machin(n) berechnet pi nach der Machinschen Formel
** auf n Dezimalstellen genau;
** pi_agm(n) benutzt zur Berechnung das arithmetisch-geometrische
** Mittel.
** pi_chud(n) benutzt eine Methode von Ramanujam-Chudnowski
**
** euler(n) berechnet e auf n Dezimalstellen.
** Der Funktionswert ist jeweils eine ganze Zahl. Diese entspricht
** gerundet 10**n-mal pi bzw. e.
*)
(*--------------------------------------------------------------*)
(*
** Berechnet zz*log(2)
*)
function log2(zz)
var
x, u, k;
begin
x := 0;
k := 0;
u := (zz * 2**16 * 2) div 3
while u /= 0 do
x := x + u div (2*k + 1);
u := u div 9;
inc(k);
end;
return x div 2**16;
end.
(*------------------------------------------------------*)
(*
** Berechnet zz*arctan(1/n), wird von pi_machin benutzt
*)
function atan1(zz,n)
var
x, u, v, k, nn;
begin
x := 0;
k := 0;
nn := n*n;
u := zz div n;
while u /= 0 do
v := u div (2*k + 1);
if even(k) then
x := x + v;
else
x := x - v;
end;
u := u div nn;
inc(k);
end;
return x;
end.
(*------------------------------------------------------*)
(*
** Berechnet pi * 10**n nach der Machinschen Formel
**
** Beispiel-Aufruf: pi_machin(1000).
*)
function pi_machin(n)
var z1, x;
begin
z1 := 10**n * 2**16;
x := 16 * atan1(z1,5) - 4 * atan1(z1,239);
return x div 2**16;
end.
(*------------------------------------------------------*)
(*
** Berechnet exp(1) * 10**n
*)
function euler(n)
var zz, x, k;
begin
zz := 10**n * 2**16;
x := zz * 2;
k := 2;
while zz /= 0 do
zz := zz div k;
x := x + zz;
inc(k);
end;
return x div 2**16;
end.
(*------------------------------------------------------*)
(*
** Berechnet pi * 10**n,
** benutzt arithmetisch-geometrisches Mittel
** quadratische Konvergenz
**
** Beispiel-Aufruf: pi_agm(1000).
*)
function pi_agm(n)
var zz;
begin
zz := 10**n * 2**16;
return piaux(zz) div 2**16;
end.
(*------------------------------------------------------*)
(*
** Hilfsfunktion fuer pi_agm
*)
function piaux(zz)
var s, a, atemp, b, c, i;
begin
s := 0;
a := zz;
b := isqrt(zz * (zz div 2));
c := (a - b) div 2;
i := 1;
while c /= 0 do
writeln("eps(",i,") = ",c/zz);
s := s + (2**i * c * c) div zz;
atemp := a;
a := (a + b) div 2;
b := isqrt(atemp * b);
c := (a - b) div 2;
inc(i);
end;
return (4*a*a) div (zz - 2*s);
end.
(*------------------------------------------------------*)
(*
** Hilfsfunktion fuer pi_chud
*)
function Saux(zz)
const
k1 = 545140134;
k2 = 13591409;
k4 = 100100025;
k5 = 327843840;
var
A, n: integer;
S: integer;
begin
A := zz*k1;
S := A * k2;
n := 1;
while A > 0 do
A := A * ((6*n-5)*(6*n-3)*(6*n-1));
A := A div (n*n*n);
A := A div (k4*k5);
if even(n) then
S := S + A * (k2 + n*k1);
else
S := S - A * (k2 + n*k1);
end;
inc(n);
end;
return S div k1;
end;
(*--------------------------------------------------------*)
(*
** pi auf n Dezimalstellen nach Chudnowsky/Ramanujan
*)
function pi_chud(n: integer): integer;
const
k3 = 640320;
k6 = 53360;
var
zz: integer;
x: integer;
begin
zz := 2**16 * 10**n;
x := isqrt(zz*zz*k3)*k6;
x := (zz * x) div Saux(zz);
return (x div 2**16);
end;
(*--------------------------------------------------------*)
|
