11package math
22
33import (
4- "github.com/ericlagergren/decimal"
54 "math"
5+
6+ "github.com/ericlagergren/decimal"
67)
78
89var (
1314 _10939058860032000 = decimal .New (10939058860032000 , 0 )
1415)
1516
16- /*
17- piChudnovskyBrothers() return PI using binary splitting on the series definition of
18- Chudnovsky Brothers
19-
20- 426880*sqrt(10005)
21- Pi = --------------------------------
22- 13591409*aSum + 545140134*bSum
23-
24- where
25- 24(6n-5)(2n-1)(6n-1)
26- a_n = ----------------------
27- (640320^3)*n^3
28-
29- and
30- aSum = sum_(n=0)^n a_n
31- bSum = sum_(n=0)^n a_n*n
32-
33-
34-
35- a(n) = 1,
36- b(n) = 1,
37- p(0) = 1, p(n) = (13591409 + 545140134n)(5 − 6n)(2n − 1)(6n − 1),
38- q(0) = 1, q(n) = (n^3)(640320^3)/24 for n > 0
39-
40- */
41-
4217func getPiA (ctx decimal.Context ) func (n uint64 ) * decimal.Big {
4318 return func (n uint64 ) * decimal.Big {
4419 //returns A + Bn
@@ -97,13 +72,13 @@ func getPiQ(ctx decimal.Context) func(n uint64) *decimal.Big {
9772 }
9873
9974 var nn , tmp decimal.Big
100- //when n is super small we can be super fast
75+ // When n is super small we can be super fast.
10176 if n < 12 {
10277 tmp .SetUint64 (n * n * n * 10939058860032000 )
10378 return & tmp
10479 }
10580
106- // and until bit larger we can still speed up a portion
81+ // And until bit larger we can still speed up a portion.
10782 if n < 2642246 {
10883 tmp .SetUint64 (n * n * n )
10984 } else {
@@ -117,21 +92,42 @@ func getPiQ(ctx decimal.Context) func(n uint64) *decimal.Big {
11792 }
11893}
11994
95+ // piChudnovskyBrothers calculates PI using binary splitting on the series
96+ // definition of the Chudnovsky Brothers' algorithm for computing pi.
97+ //
98+ // 426880*sqrt(10005)
99+ // Pi = --------------------------------
100+ // 13591409*aSum + 545140134*bSum
101+ //
102+ // where
103+ // 24(6n-5)(2n-1)(6n-1)
104+ // a_n = ----------------------
105+ // (640320^3)*n^3
106+ //
107+ // and
108+ // aSum = sum_(n=0)^n a_n
109+ // bSum = sum_(n=0)^n a_n*n
110+ //
111+ // a(n) = 1,
112+ // b(n) = 1,
113+ // p(0) = 1, p(n) = (13591409 + 545140134n)(5 − 6n)(2n − 1)(6n − 1),
114+ // q(0) = 1, q(n) = (n^3)(640320^3)/24 for n > 0
115+ //
120116func piChudnovskyBrothers (z * decimal.Big , ctx decimal.Context ) * decimal.Big {
121- //since this algorithm's rate of convergence is static calculating the number of
122- // iteration required will always be faster than the dynamic method of binarysplit
117+ // Since this algorithm's rate of convergence is static, calculating the
118+ // number of iteration required will always be faster than the dynamic
119+ // method of binarysplit.
120+
121+ calcPrec := ctx .Precision + 16
122+ niters := uint64 (math .Ceil (float64 (calcPrec )/ 14 )) + 1
123+ ctx2 := decimal.Context {Precision : calcPrec }
124+
123125 var value decimal.Big
124- extraPrecision := 16
125- calculatingPrecision := ctx .Precision + extraPrecision
126- iterationNeeded := uint64 (math .Ceil (float64 (calculatingPrecision / 14.0 ))) + 1
127- ctx2 := decimal.Context {Precision : calculatingPrecision }
126+ BinarySplit (& value , ctx2 , 0 , niters , getPiA (ctx2 ), getPiP (ctx2 ), getPiB (), getPiQ (ctx2 ))
127+ s := Sqrt (decimal .WithPrecision (calcPrec ), _10005 )
128128
129129 var tmp decimal.Big
130-
131- BinarySplit (& value , ctx2 , 0 , iterationNeeded , getPiA (ctx2 ), getPiP (ctx2 ), getPiB (), getPiQ (ctx2 ))
132- s := Sqrt (decimal .WithPrecision (calculatingPrecision ), _10005 )
133130 ctx2 .Mul (& tmp , _426880 , s )
134131 ctx2 .Quo (z , & tmp , & value )
135-
136132 return z .Round (ctx .Precision )
137133}
0 commit comments