{-
	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@]	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,
	spread,
-- ** Accessors
	getArithmeticMean,
	getGeometricMean,
-- ** Predicates
	isValid
) where

import			Control.Arrow((&&&))
import qualified	Control.Parallel.Strategies
import qualified	Factory.Math.Precision	as Math.Precision
import qualified	Factory.Math.SquareRoot	as Math.SquareRoot

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

-- | The type of the /geometric mean/; <http://en.wikipedia.org/wiki/Geometric_mean>.
type GeometricMean	= 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 :: Rational -> 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 (
		Control.Parallel.Strategies.withStrategy (
			Control.Parallel.Strategies.parTuple2 Control.Parallel.Strategies.rdeepseq Control.Parallel.Strategies.rdeepseq
		) . (simplify . findArithmeticMean &&& simplify . findGeometricMean)
	) agm

-- | Returns the bounds within which the 'AGM' has been constrained.
spread :: AGM -> Rational
spread	= uncurry (-)

-- | Checks that both /means/ are positive, as required for the /geometric mean/ to be consistently /real/.
isValid :: AGM -> Bool
isValid (a, g)	= all (>= 0) [a, g]