{-
	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 <https://www.gnu.org/licenses/>.
-}
{- |
 [@AUTHOR@]	Dr. Alistair Ward

 [@DESCRIPTION@]

	* Implements the /Brent-Salamin/ (AKA /Gauss-Legendre/) algorithm;
		<https://en.wikipedia.org/wiki/Gauss%E2%80%93Legendre_algorithm>,
		<https://mathworld.wolfram.com/Brent-SalaminFormula.html>,
		<https://www.pi314.net/eng/salamin.php>.

	* The precision of the result approximately doubles for each iteration.

 [@CAVEAT@]	Assumptions on the convergence-rate result in rounding-errors, when only a small number of digits are requested.
-}

module Factory.Math.Implementations.Pi.AGM.BrentSalamin(
-- * Functions
	openR
) where

import			Control.Arrow((&&&))
import qualified	Factory.Math.ArithmeticGeometricMean	as Math.ArithmeticGeometricMean
import qualified	Factory.Math.Power			as Math.Power
import qualified	Factory.Math.Precision			as Math.Precision
import qualified	Factory.Math.SquareRoot			as Math.SquareRoot

{- |
	* Returns /Pi/, accurate to the specified number of decimal digits.

	* This algorithm is based on the /arithmetic-geometric/ mean of @1@ and @(1 / sqrt 2)@,
	but there are many confusingly similar formulations.
	The algorithm I've used here, where @a@ is the /arithmetic mean/ and @g@ is the /geometric mean/, is equivalent to other common formulations:

>		pi = (a[N-1] + g[N-1])^2 / (1 - sum [2^n * (a[n] - g[n])^2])			where n = [0 .. N-1]
>		=> 4*a[N]^2 / (1 - sum [2^n * (a[n]^2 - 2*a[n]*g[n] + g[n]^2)])
>		=> 4*a[N]^2 / (1 - sum [2^n * (a[n]^2 + 2*a[n]*g[n] + g[n]^2 - 4*a[n]*g[n])])
>		=> 4*a[N]^2 / (1 - sum [2^n * ((a[n] + g[n])^2 - 4*a[n]*g[n])])
>		=> 4*a[N]^2 / (1 - sum [2^(n-1) * 4 * (a[n-1]^2 - g[n-1]^2)])			where n = [1 .. N]
>		=> 4*a[N]^2 / (1 - sum [2^(n+1) * (a[n-1]^2 - g[n-1]^2)])

-}
openR :: Math.SquareRoot.Algorithmic squareRootAlgorithm => squareRootAlgorithm -> Math.Precision.DecimalDigits -> Rational
openR :: squareRootAlgorithm -> DecimalDigits -> Rational
openR squareRootAlgorithm
squareRootAlgorithm DecimalDigits
decimalDigits	= (Rational -> Rational -> Rational)
-> (Rational, Rational) -> Rational
forall a b c. (a -> b -> c) -> (a, b) -> c
uncurry Rational -> Rational -> Rational
forall a. Fractional a => a -> a -> a
(/) ((Rational, Rational) -> Rational)
-> ([(Rational, 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.square (Rational -> Rational)
-> ([(Rational, Rational)] -> Rational)
-> [(Rational, Rational)]
-> Rational
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (Rational -> Rational -> Rational)
-> (Rational, Rational) -> Rational
forall a b c. (a -> b -> c) -> (a, b) -> c
uncurry Rational -> Rational -> Rational
forall a. Num a => a -> a -> a
(+) ((Rational, Rational) -> Rational)
-> ([(Rational, Rational)] -> (Rational, Rational))
-> [(Rational, Rational)]
-> Rational
forall b c a. (b -> c) -> (a -> b) -> a -> c
. [(Rational, Rational)] -> (Rational, Rational)
forall a. [a] -> a
last ([(Rational, Rational)] -> Rational)
-> ([(Rational, Rational)] -> Rational)
-> [(Rational, Rational)]
-> (Rational, Rational)
forall (a :: * -> * -> *) b c c'.
Arrow a =>
a b c -> a b c' -> a b (c, c')
&&& Rational -> Rational
forall n. Num n => n -> n
negate (Rational -> Rational)
-> ([(Rational, Rational)] -> Rational)
-> [(Rational, Rational)]
-> Rational
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Rational -> Rational
forall a. Enum a => a -> a
pred (Rational -> Rational)
-> ([(Rational, Rational)] -> Rational)
-> [(Rational, Rational)]
-> Rational
forall b c a. (b -> c) -> (a -> b) -> a -> c
. [Rational] -> Rational
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum ([Rational] -> Rational)
-> ([(Rational, Rational)] -> [Rational])
-> [(Rational, Rational)]
-> Rational
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (Rational -> Rational -> Rational)
-> [Rational] -> [Rational] -> [Rational]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith Rational -> Rational -> Rational
forall a. Num a => a -> a -> a
(*) ((Rational -> Rational) -> Rational -> [Rational]
forall a. (a -> a) -> a -> [a]
iterate (Rational -> Rational -> Rational
forall a. Num a => a -> a -> a
* Rational
2) Rational
1) ([Rational] -> [Rational])
-> ([(Rational, Rational)] -> [Rational])
-> [(Rational, Rational)]
-> [Rational]
forall b c a. (b -> c) -> (a -> b) -> a -> c
. ((Rational, Rational) -> Rational)
-> [(Rational, Rational)] -> [Rational]
forall a b. (a -> b) -> [a] -> [b]
map (Rational -> Rational
forall n. Num n => n -> n
Math.Power.square (Rational -> Rational)
-> ((Rational, Rational) -> Rational)
-> (Rational, Rational)
-> Rational
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (Rational, Rational) -> Rational
Math.ArithmeticGeometricMean.spread)
 ) ([(Rational, Rational)] -> (Rational, Rational))
-> ([(Rational, Rational)] -> [(Rational, Rational)])
-> [(Rational, Rational)]
-> (Rational, Rational)
forall b c a. (b -> c) -> (a -> b) -> a -> c
. DecimalDigits -> [(Rational, Rational)] -> [(Rational, Rational)]
forall a. DecimalDigits -> [a] -> [a]
take (
	DecimalDigits -> DecimalDigits -> DecimalDigits -> DecimalDigits
forall i.
Integral i =>
DecimalDigits -> DecimalDigits -> DecimalDigits -> i
Math.Precision.getIterationsRequired DecimalDigits
Math.Precision.quadraticConvergence DecimalDigits
1 DecimalDigits
decimalDigits
 ) ([(Rational, Rational)] -> Rational)
-> [(Rational, Rational)] -> Rational
forall a b. (a -> b) -> a -> b
$ squareRootAlgorithm
-> DecimalDigits -> (Rational, Rational) -> [(Rational, Rational)]
forall squareRootAlgorithm.
Algorithmic squareRootAlgorithm =>
squareRootAlgorithm
-> DecimalDigits -> (Rational, Rational) -> [(Rational, Rational)]
Math.ArithmeticGeometricMean.convergeToAGM squareRootAlgorithm
squareRootAlgorithm DecimalDigits
decimalDigits (Rational
1, squareRootAlgorithm -> DecimalDigits -> Rational -> Rational
forall algorithm operand.
(Algorithmic algorithm, Real operand, Show operand) =>
algorithm -> DecimalDigits -> operand -> Rational
Math.SquareRoot.squareRoot squareRootAlgorithm
squareRootAlgorithm DecimalDigits
decimalDigits (Rational -> Rational
forall a. Fractional a => a -> a
recip Rational
2 :: Rational))