{-# LANGUAGE CPP #-}
{-
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 .
-}
{- |
[@AUTHOR@] Dr. Alistair Ward
[@DESCRIPTION@] Determines the /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 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/; .
type ArithmeticMean = Data.Ratio.Rational
-- | The type of the /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.
spread :: AGM -> Data.Ratio.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]