module Factory.Math.Implementations.Pi.Borwein.Borwein1993(
series
) where
import qualified Data.Ratio
import Data.Ratio((%))
import qualified Factory.Math.Factorial as Math.Factorial
import qualified Factory.Math.Implementations.Factorial as Math.Implementations.Factorial
import qualified Factory.Math.Implementations.Pi.Borwein.Series as Math.Implementations.Pi.Borwein.Series
import qualified Factory.Math.Power as Math.Power
import qualified Factory.Math.Precision as Math.Precision
import qualified Factory.Math.SquareRoot as Math.SquareRoot
series :: (Math.SquareRoot.Algorithm squareRootAlgorithm, Math.Factorial.Algorithm factorialAlgorithm) => Math.Implementations.Pi.Borwein.Series.Series squareRootAlgorithm factorialAlgorithm
series = Math.Implementations.Pi.Borwein.Series.MkSeries {
Math.Implementations.Pi.Borwein.Series.terms = \squareRootAlgorithm factorialAlgorithm decimalDigits -> let
simplify, squareRoot :: Data.Ratio.Rational -> Data.Ratio.Rational
simplify = Math.Precision.simplify (decimalDigits 1 )
squareRoot = simplify . Math.SquareRoot.squareRoot squareRootAlgorithm decimalDigits
sqrt5, a, b, c3 :: Data.Ratio.Rational
sqrt5 = squareRoot 5
a = 63365028312971999585426220 + sqrt5 * (28337702140800842046825600 + 384 * squareRoot (10891728551171178200467436212395209160385656017 + 4870929086578810225077338534541688721351255040 * sqrt5))
b = 7849910453496627210289749000 + 3510586678260932028965606400 * sqrt5 + 2515968 * squareRoot (3110 * (6260208323789001636993322654444020882161 + 2799650273060444296577206890718825190235 * sqrt5))
c3 = simplify . Math.Power.cube $ negate 214772995063512240 sqrt5 * (96049403338648032 + 1296 * squareRoot (10985234579463550323713318473 + 4912746253692362754607395912 * sqrt5))
in (
squareRoot $ negate c3,
zipWith (
\n power -> (
Math.Implementations.Factorial.risingFactorial (3 * n + 1) (3 * n) % Math.Power.cube (Math.Factorial.factorial factorialAlgorithm n)
) * (
(a + b * fromIntegral n) / power
)
) [0 :: Integer ..] $ iterate (* c3) 1
),
Math.Implementations.Pi.Borwein.Series.convergenceRate = 10 ** negate 50 --Empirical.
}