```{-# LANGUAGE CPP #-}
{-
Copyright (C) 2011 Dr. Alistair Ward

This program is free software: you can redistribute it and/or modify
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@]	Determines the /Arithmetic-geometric mean/; <http://en.wikipedia.org/wiki/Arithmetic-geometric_mean>.
-}

module Factory.Math.ArithmeticGeometricMean(
-- * Types
-- ** Type-synonyms
ArithmeticMean,
GeometricMean,
AGM,
-- * Functions
convergeToAGM,
-- ** Accessors
getArithmeticMean,
getGeometricMean,
-- ** Predicates
isValid
) where

import			Control.Arrow((&&&))
import qualified	Data.Ratio
import qualified	Factory.Math.Precision	as Math.Precision
import qualified	Factory.Math.SquareRoot	as Math.SquareRoot

#if MIN_VERSION_parallel(3,0,0)
import qualified	Control.Parallel.Strategies
#endif

-- | The type of the /arithmetic mean/; <http://en.wikipedia.org/wiki/Arithmetic_mean>.
type ArithmeticMean	= Data.Ratio.Rational

-- | The type of the /geometric mean/; <http://en.wikipedia.org/wiki/Geometric_mean>.
type GeometricMean	= Data.Ratio.Rational

-- | Encapsulates both /arithmetic/ and /geometric/ means.
type AGM	= (ArithmeticMean, GeometricMean)

-- | Accessor.
{-# INLINE getArithmeticMean #-}
getArithmeticMean :: AGM -> ArithmeticMean
getArithmeticMean	= fst

-- | Accessor.
{-# INLINE getGeometricMean #-}
getGeometricMean :: AGM -> GeometricMean
getGeometricMean	= snd

-- | Returns an infinite list which converges on the /Arithmetic-geometric mean/.
convergeToAGM :: Math.SquareRoot.Algorithmic squareRootAlgorithm => squareRootAlgorithm -> Math.Precision.DecimalDigits -> AGM -> [AGM]
convergeToAGM squareRootAlgorithm decimalDigits agm
| decimalDigits <= 0	= error \$ "Factory.Math.ArithmeticGeometricMean.convergeToAGM:\tinvalid number of decimal digits; " ++ show decimalDigits
| not \$ isValid agm	= error \$ "Factory.Math.ArithmeticGeometricMean.convergeToAGM:\tboth means must be positive for a real geometric mean; " ++ show agm
| spread agm == 0	= repeat agm
| otherwise		= let
simplify :: Data.Ratio.Rational -> Data.Ratio.Rational
simplify	= Math.Precision.simplify (pred decimalDigits {-ignore single integral digit-})	--This makes a gigantic difference to performance.

findArithmeticMean :: AGM -> ArithmeticMean
findArithmeticMean	= (/ 2) . uncurry (+)

findGeometricMean :: AGM -> GeometricMean
findGeometricMean	= Math.SquareRoot.squareRoot squareRootAlgorithm decimalDigits . uncurry (*)
in iterate (
#if MIN_VERSION_parallel(3,0,0)
Control.Parallel.Strategies.withStrategy (
Control.Parallel.Strategies.parTuple2 Control.Parallel.Strategies.rdeepseq Control.Parallel.Strategies.rdeepseq
) .
#endif
(simplify . findArithmeticMean &&& simplify . findGeometricMean)
) agm

-- | Returns the bounds within which the 'AGM' has been constrained.