1/* Calculating Pi to multiple precision using advanced methods */
2
3/* Defined: PiMethod0(), PiMethod1(), PiMethod2(), PiBrentSalamin(), PiBorwein() */
4
5// Reference method: just use Newton's method all the time, no complicated logic to select precision steps. Slightly slower than method 1 but a lot simpler. This is implemented in Internal'Pi()
6
7PiMethod0() := [
8	Local(result, delta, k, Epsilon, prec, prec1, curprec);
9	prec := Builtin'Precision'Get();	// full required precision
10	prec1 := Ceil(N(prec/3));	// precision of the last-but-one iteration
11
12	/* initial approximation */
13	result := 3.14159265358979323846;
14	curprec := 20;
15	Builtin'Precision'Set(curprec);
16	For(k:=0, curprec < prec1, k:=k+1) [
17		curprec := Min(prec1, curprec * 3);
18		Builtin'Precision'Set(curprec);
19		Echo({"Iteration ", k, " setting precision to ", curprec});
20		result := Time(MathAdd(result, MathSin(result)));
21	];
22	// last iteration -- do by hand
23	Builtin'Precision'Set(prec);	// restore full precision
24	Echo("Iteration ", k, " setting precision to ", Builtin'Precision'Get());
25	result := Time(MathAdd(result, MathSin(result)));
26	Echo({"variable precision Newton's method: ", k, "iterations"});
27	result;
28];
29
30/* Brute-force method 1: start near 3.14159... and iterate using Nth order
31Newton method: x := x + ( sin(x) + 1/6*sin(x)^3 + 3/40*sin(x)^5 +
325/112*sin(x)^7 + ...) i.e. the Taylor series for arcsin but cut at a finite
33point. Convergence is of order of the next term, i.e. x^9, but need to evaluate
34Sin() at each step. However, we don't need to evaluate them to full precision
35each time, because each iteration will correct any accumulated errors. In fact,
36first iterations can be computed with significantly lower precision than the
37final result. This makes method1 the fastest for Yacas internal math. */
38
39PiMethod1() := [
40	Local(result, delta, deltasq, k, Epsilon, prec, curprec);
41	prec := Builtin'Precision'Get();
42	N([
43    /* initial approximation */
44    curprec := 20;
45    Builtin'Precision'Set(curprec);
46    result := 3.14159265358979323846;
47    /* right now we do all but the last iteration using the 8th order scheme, and the last iteration is done using the 2nd order scheme. However it would be faster to use a very high-order scheme first, then a smaller-order scheme, etc., because it's faster to do multiplications at low precision.
48    */
49    For(k:=0, curprec*3 < prec, k := k+1) [
50      curprec := Min(Ceil((prec/3)), curprec * 9);
51      Builtin'Precision'Set(curprec);
52      Echo("Iteration ", k, " setting precision to ", Builtin'Precision'Get());
53      delta := MathSin(result);
54      deltasq := (delta*delta);
55      result := Time(result + delta*(1 + deltasq*(1/6 + deltasq*(3/40 + deltasq*5/112))));
56    ];
57    // now do the last iteration
58    Builtin'Precision'Set(prec);
59    k := k+1;
60    Echo("Iteration ", k, " setting precision to ", Builtin'Precision'Get());
61    result := Time(MathAdd(result, MathSin(result)));
62    Echo({"8th order Newton's method: ", k, "iterations"});
63  ]);
64	result;
65];
66
67/* Brute-force method 2: evaluate full series for arctan */
68/* x0 := 3.14159... and Pi = x0 - ( tan(x0) - tan(x0)^3/3 + tan(x0)^5/5 +...) i.e. the Taylor series for arctan - go until it converges to Pi. Convergence is linear but unlike method 1, we don't need to evaluate Sin() and Cos() at every step, and we can start at a very good initial approximation to cut computing time.
69*/
70
71PiMethod2() := [
72	Local(result, delta, tansq, k, Epsilon);
73	N([
74    Epsilon := (2*10 ^ (-Builtin'Precision'Get()));
75
76    /* initial approximation */
77    result := 3.141592653589793;
78    delta := (-Tan(result));
79    tansq := (delta^2);
80    k := 0;
81
82    While(Abs(delta) > Epsilon) [
83  //		Echo(result);
84      result := (result + delta/(2*k+1));
85  //		Echo(delta, k);
86      delta := (-delta * tansq);
87      k := k+1;
88    ];
89    Echo({"Brute force method 2 (ArcTan series): ", k, "iterations"});
90  ]);
91	result;
92];
93
94/* Method due to Brent and Salamin (1976) */
95PiBrentSalamin() := [
96	Local(a, b, c, s, k, p, result, Epsilon);
97	Epsilon := N(2*10 ^ (-Builtin'Precision'Get()));
98
99	/* initialization */
100	a := 1; b := N(1/Sqrt(2)); s := N(1/2); k := 0;
101	/* this is just to make sure we stop - the algorithm does not require initialization of p */
102	p := 0; result := 1;
103	/* repeat until precision is saturated */
104	While(Abs(p-result) >= Epsilon) [
105		k := k+1;
106		result := p;
107		/* arithmetic and geometric mean */
108		{a, b} := {N((a+b)/2), N(Sqrt(a*b))};
109		/* difference between them is scaled by 2^k */
110		s := N(s - 2^k*(a^2-b^2));
111		p := N(2*a^2/s);
112	];
113	Echo({"Brent and Salamin's algorithm: ", k, "iterations"});
114
115	result;
116];
117
118/* Method due to Borwein (c. 1988) -- "quartic" method */
119PiBorwein() := [
120	Local(a, y, y4s, k, result, Epsilon);
121	Epsilon := N(2*10 ^ (-Builtin'Precision'Get()));
122
123	/* initialization */
124	a:=N(6-4*Sqrt(2)); y := N(Sqrt(2)-1); k := 0;
125	result := 0;
126	/* repeat until precision is saturated */
127	While(Abs(a-result) >= Epsilon) [
128		result := a;
129		/* precompute (1-y^4)^(1/4) */
130		y4s:=N(Sqrt(Sqrt(1-y^4)));
131		/* black magic */
132		y := N((1-y4s)/(1+y4s));
133		/* more black magic */
134		a := a*(1+y)^4-2^(2*k+3)*y*(1+y+y^2);
135		k := k+1;
136	];
137	/* {a} will converge to 1/Pi */
138	result := N(1/result);
139
140	Echo({"Borwein's quartic algorithm: ", k, "iterations"});
141	result;
142];
143
144// iterate x := x + Cos(x) + 1/6 *Cos(x)^3 + ... to converge to x=Pi/2
145PiMethod3() :=
146[
147	Local(result, delta, deltasq, k, order, prec, curprec);
148	order := 13;	// order of approximation
149	prec := Builtin'Precision'Get();
150  N([
151    /* initial approximation */
152    curprec := 20;
153    Builtin'Precision'Set(curprec);
154    result := 3.14159265358979323846*0.5;
155    // find optimal initial precision
156    For(k:=prec, k>=curprec, k:=Div(k,order)+2) True;
157    If(k<5, curprec:=5, curprec:=k);
158  //	Echo("initial precision", curprec);
159    // now k is the iteration counter
160    For(k:=0, curprec < prec, k := k+1) [
161    // at this iteration we know the result to curprec digits
162      curprec := Min(prec, curprec * order-2);	// 2 guard digits
163      Builtin'Precision'Set(curprec+2);
164      Echo("Iteration ", k, " setting precision to ", Builtin'Precision'Get());
165  //		Echo("old result=", MathCos(result));
166      Time([
167      delta := MathCos(result);
168      ]);
169      Time([
170      deltasq := MathMultiply(delta,delta);
171      ]);
172      result := Time(result + delta*(1 + deltasq*(1/6 + deltasq*(3/40 + deltasq*(5/112 + deltasq*(35/1152 + (deltasq*63)/2816))))));
173    ];
174    Echo({"Method 3, using Pi/2 and order", order, ":", k, "iterations"});
175  ]);
176	result*2;
177];
178
179PiChudnovsky() :=
180[	// use the Ramanujan series found by Chudnovsky brothers
181	Local(A, B, C, n, result, term);
182	A:=13591409; B:=545140134; C:=640320; // black magic, Rama, Rama, Ramanujan
183	prec := Builtin'Precision'Get();
184  N([
185    n:=Div(prec*479,6793)+1;	// n> P*Ln(10)/(3*Ln(C/12))
186    Echo({"Method: Chudnovsky, using ", n, " terms"});
187    Builtin'Precision'Set(prec+IntLog(n,10)+5);
188    result := (A+n*B);
189    While(n>0)
190    [
191  //	Echo(n,result);
192      result := A+(n-1)*B-24*(6*n-1)*(2*n-1)*(6*n-5) /(C*n)^3 *result;
193      n--;
194    ];
195    result := C/12*Sqrt(C)/Abs(result);
196  ]);
197	Builtin'Precision'Set(prec);
198	RoundTo(result,prec);
199];
200
201BenchmarkPi(prec) :=
202[
203	Local(result);
204	GlobalPush(Builtin'Precision'Get());
205	Builtin'Precision'Set(prec);
206
207	result := {
208		Time(MathPi()),
209		Time(PiMethod0()),
210		Time(PiMethod1()),
211		Time(PiMethod2()),
212		Time(PiMethod3()),
213//		Time(PiMethod4()),
214		Time(PiChudnovsky()),
215		Time(PiBrentSalamin()),
216		Time(PiBorwein()),
217	};
218	result := N(Sin(result));
219	Builtin'Precision'Set(GlobalPop());
220	result;
221];
222