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