{-
	Copyright (C) 2011 Dr. Alistair Ward

	This program is free software: you can redistribute it and/or modify
	it under the terms of the GNU General Public License as published by
	the Free Software Foundation, either version 3 of the License, or
	(at your option) any later version.

	This program is distributed in the hope that it will be useful,
	but WITHOUT ANY WARRANTY; without even the implied warranty of
	MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
	GNU General Public License for more details.

	You should have received a copy of the GNU General Public License
	along with this program.  If not, see <http://www.gnu.org/licenses/>.
-}
{- |
 [@AUTHOR@]	Dr. Alistair Ward

 [@DESCRIPTION@]	Defines the /Borwein/ series for /Pi/; <https://en.wikipedia.org/wiki/Borwein%27s_algorithm#Jonathan_Borwein_and_Peter_Borwein.27s_Version_.281993.29>
-}

module Factory.Math.Implementations.Pi.Borwein.Borwein1993(
-- * Constants
	series
) where

-- import		Control.Arrow((***))
import			Data.Ratio((%))
-- import		Factory.Data.PrimeFactors((>*<), (>/<), (>^))
-- import qualified	Factory.Data.PrimeFactors			as Data.PrimeFactors
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

-- | Defines the parameters of the /Borwein/ series.
series :: (Math.SquareRoot.Algorithmic squareRootAlgorithm, Math.Factorial.Algorithmic factorialAlgorithm) => Math.Implementations.Pi.Borwein.Series.Series squareRootAlgorithm factorialAlgorithm
series :: Series squareRootAlgorithm factorialAlgorithm
series = MkSeries :: forall squareRootAlgorithm factorialAlgorithm.
(squareRootAlgorithm
 -> factorialAlgorithm -> DecimalDigits -> (Rational, [Rational]))
-> ConvergenceRate -> Series squareRootAlgorithm factorialAlgorithm
Math.Implementations.Pi.Borwein.Series.MkSeries {
	terms :: squareRootAlgorithm
-> factorialAlgorithm -> DecimalDigits -> (Rational, [Rational])
Math.Implementations.Pi.Borwein.Series.terms			= \squareRootAlgorithm
squareRootAlgorithm factorialAlgorithm
factorialAlgorithm DecimalDigits
decimalDigits -> let
		simplify, squareRoot :: Rational -> Rational
		simplify :: Rational -> Rational
simplify	= DecimalDigits -> Rational -> Rational
forall operand.
RealFrac operand =>
DecimalDigits -> operand -> Rational
Math.Precision.simplify (DecimalDigits -> Rational -> Rational)
-> DecimalDigits -> Rational -> Rational
forall a b. (a -> b) -> a -> b
$ DecimalDigits -> DecimalDigits
forall a. Enum a => a -> a
pred DecimalDigits
decimalDigits {-ignore single integral digit-}	-- This makes a gigantic difference to performance.
		squareRoot :: Rational -> Rational
squareRoot	= Rational -> Rational
simplify (Rational -> Rational)
-> (Rational -> Rational) -> Rational -> Rational
forall b c a. (b -> c) -> (a -> b) -> a -> c
. squareRootAlgorithm -> DecimalDigits -> Rational -> Rational
forall algorithm operand.
(Algorithmic algorithm, Real operand, Show operand) =>
algorithm -> DecimalDigits -> operand -> Rational
Math.SquareRoot.squareRoot squareRootAlgorithm
squareRootAlgorithm DecimalDigits
decimalDigits

		sqrt5, a, b, c3 :: Rational
		sqrt5 :: Rational
sqrt5	= Rational -> Rational
squareRoot Rational
5

		a :: Rational
a	= Rational
63365028312971999585426220 Rational -> Rational -> Rational
forall a. Num a => a -> a -> a
+ Rational
sqrt5 Rational -> Rational -> Rational
forall a. Num a => a -> a -> a
* (Rational
28337702140800842046825600 Rational -> Rational -> Rational
forall a. Num a => a -> a -> a
+ Rational
384 Rational -> Rational -> Rational
forall a. Num a => a -> a -> a
* Rational -> Rational
squareRoot (Rational
10891728551171178200467436212395209160385656017 Rational -> Rational -> Rational
forall a. Num a => a -> a -> a
+ Rational
4870929086578810225077338534541688721351255040 Rational -> Rational -> Rational
forall a. Num a => a -> a -> a
* Rational
sqrt5))
		b :: Rational
b	= Rational
7849910453496627210289749000 Rational -> Rational -> Rational
forall a. Num a => a -> a -> a
+ Rational
3510586678260932028965606400 Rational -> Rational -> Rational
forall a. Num a => a -> a -> a
* Rational
sqrt5 Rational -> Rational -> Rational
forall a. Num a => a -> a -> a
+ Rational
2515968 Rational -> Rational -> Rational
forall a. Num a => a -> a -> a
* Rational -> Rational
squareRoot (Rational
3110 Rational -> Rational -> Rational
forall a. Num a => a -> a -> a
* (Rational
6260208323789001636993322654444020882161 Rational -> Rational -> Rational
forall a. Num a => a -> a -> a
+ Rational
2799650273060444296577206890718825190235 Rational -> Rational -> Rational
forall a. Num a => a -> a -> a
* Rational
sqrt5))
		c3 :: Rational
c3	= Rational -> Rational
simplify (Rational -> Rational)
-> (Rational -> Rational) -> Rational -> Rational
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Rational -> Rational
forall n. Num n => n -> n
Math.Power.cube (Rational -> Rational) -> Rational -> Rational
forall a b. (a -> b) -> a -> b
$ Rational -> Rational
forall n. Num n => n -> n
negate Rational
214772995063512240 Rational -> Rational -> Rational
forall a. Num a => a -> a -> a
- Rational
sqrt5 Rational -> Rational -> Rational
forall a. Num a => a -> a -> a
* (Rational
96049403338648032 Rational -> Rational -> Rational
forall a. Num a => a -> a -> a
+ Rational
1296 Rational -> Rational -> Rational
forall a. Num a => a -> a -> a
* Rational -> Rational
squareRoot (Rational
10985234579463550323713318473 Rational -> Rational -> Rational
forall a. Num a => a -> a -> a
+ Rational
4912746253692362754607395912 Rational -> Rational -> Rational
forall a. Num a => a -> a -> a
* Rational
sqrt5))
	in (
		Rational -> Rational
squareRoot (Rational -> Rational) -> Rational -> Rational
forall a b. (a -> b) -> a -> b
$ Rational -> Rational
forall n. Num n => n -> n
negate Rational
c3,	-- The factor into which the series must be divided, to yield Pi.
		(Integer -> Rational -> Rational)
-> [Integer] -> [Rational] -> [Rational]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith (
{-
			\n power -> let
				product'	= Data.PrimeFactors.product' (recip 2) 10
			in uncurry (/) . (
				(* (a + b * fromIntegral n)) . fromIntegral . product' *** (* power) . fromIntegral . product'
			) $ Math.Implementations.Factorial.primeFactors (6 * n) >/< (
				Math.Implementations.Factorial.primeFactors (3 * n) >*< Math.Implementations.Factorial.primeFactors n >^ 3
			)
-}
			\Integer
n Rational
power -> (
				Integer -> Integer -> Integer
forall i. (Integral i, Show i) => i -> i -> i
Math.Implementations.Factorial.risingFactorial (Integer -> Integer
forall a. Enum a => a -> a
succ (Integer -> Integer) -> Integer -> Integer
forall a b. (a -> b) -> a -> b
$ Integer
3 Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
* Integer
n) (Integer
3 Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
* Integer
n) Integer -> Integer -> Rational
forall a. Integral a => a -> a -> Ratio a
% Integer -> Integer
forall n. Num n => n -> n
Math.Power.cube (factorialAlgorithm -> Integer -> Integer
forall algorithm i.
(Algorithmic algorithm, Integral i, Show i) =>
algorithm -> i -> i
Math.Factorial.factorial factorialAlgorithm
factorialAlgorithm Integer
n)
			) Rational -> Rational -> Rational
forall a. Num a => a -> a -> a
* (
				(Rational
a Rational -> Rational -> Rational
forall a. Num a => a -> a -> a
+ Rational
b Rational -> Rational -> Rational
forall a. Num a => a -> a -> a
* Integer -> Rational
forall a b. (Integral a, Num b) => a -> b
fromIntegral Integer
n) Rational -> Rational -> Rational
forall a. Fractional a => a -> a -> a
/ Rational
power
			)
		) [Integer
0 :: Integer ..] ([Rational] -> [Rational]) -> [Rational] -> [Rational]
forall a b. (a -> b) -> a -> b
$ (Rational -> Rational) -> Rational -> [Rational]
forall a. (a -> a) -> a -> [a]
iterate (Rational -> Rational -> Rational
forall a. Num a => a -> a -> a
* Rational
c3) Rational
1
	),
	convergenceRate :: ConvergenceRate
Math.Implementations.Pi.Borwein.Series.convergenceRate		= ConvergenceRate
10 ConvergenceRate -> ConvergenceRate -> ConvergenceRate
forall a. Floating a => a -> a -> a
** ConvergenceRate -> ConvergenceRate
forall n. Num n => n -> n
negate ConvergenceRate
50	-- Empirical.
}