{-
	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/; <https://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/; <https://en.wikipedia.org/wiki/Arithmetic_mean>.
type ArithmeticMean	= Rational

-- | The type of the /geometric mean/; <https://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 :: AGM -> ArithmeticMean
getArithmeticMean	= AGM -> ArithmeticMean
forall a b. (a, b) -> a
fst

-- | Accessor.
{-# INLINE getGeometricMean #-}
getGeometricMean :: AGM -> GeometricMean
getGeometricMean :: AGM -> ArithmeticMean
getGeometricMean	= AGM -> ArithmeticMean
forall a b. (a, b) -> b
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 -> [AGM]
convergeToAGM squareRootAlgorithm
squareRootAlgorithm DecimalDigits
decimalDigits AGM
agm
	| DecimalDigits
decimalDigits DecimalDigits -> DecimalDigits -> Bool
forall a. Ord a => a -> a -> Bool
<= DecimalDigits
0	= [Char] -> [AGM]
forall a. HasCallStack => [Char] -> a
error ([Char] -> [AGM]) -> [Char] -> [AGM]
forall a b. (a -> b) -> a -> b
$ [Char]
"Factory.Math.ArithmeticGeometricMean.convergeToAGM:\tinvalid number of decimal digits; " [Char] -> [Char] -> [Char]
forall a. [a] -> [a] -> [a]
++ DecimalDigits -> [Char]
forall a. Show a => a -> [Char]
show DecimalDigits
decimalDigits
	| Bool -> Bool
not (Bool -> Bool) -> Bool -> Bool
forall a b. (a -> b) -> a -> b
$ AGM -> Bool
isValid AGM
agm	= [Char] -> [AGM]
forall a. HasCallStack => [Char] -> a
error ([Char] -> [AGM]) -> [Char] -> [AGM]
forall a b. (a -> b) -> a -> b
$ [Char]
"Factory.Math.ArithmeticGeometricMean.convergeToAGM:\tboth means must be positive for a real geometric mean; " [Char] -> [Char] -> [Char]
forall a. [a] -> [a] -> [a]
++ AGM -> [Char]
forall a. Show a => a -> [Char]
show AGM
agm
	| AGM -> ArithmeticMean
spread AGM
agm ArithmeticMean -> ArithmeticMean -> Bool
forall a. Eq a => a -> a -> Bool
== ArithmeticMean
0	= AGM -> [AGM]
forall a. a -> [a]
repeat AGM
agm
	| Bool
otherwise		= let
		simplify :: Rational -> Rational
		simplify :: ArithmeticMean -> ArithmeticMean
simplify	= DecimalDigits -> ArithmeticMean -> ArithmeticMean
forall operand.
RealFrac operand =>
DecimalDigits -> operand -> ArithmeticMean
Math.Precision.simplify (DecimalDigits -> DecimalDigits
forall a. Enum a => a -> a
pred DecimalDigits
decimalDigits {-ignore single integral digit-})	-- This makes a gigantic difference to performance.

		findArithmeticMean :: AGM -> ArithmeticMean
		findArithmeticMean :: AGM -> ArithmeticMean
findArithmeticMean	= (ArithmeticMean -> ArithmeticMean -> ArithmeticMean
forall a. Fractional a => a -> a -> a
/ ArithmeticMean
2) (ArithmeticMean -> ArithmeticMean)
-> (AGM -> ArithmeticMean) -> AGM -> ArithmeticMean
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (ArithmeticMean -> ArithmeticMean -> ArithmeticMean)
-> AGM -> ArithmeticMean
forall a b c. (a -> b -> c) -> (a, b) -> c
uncurry ArithmeticMean -> ArithmeticMean -> ArithmeticMean
forall a. Num a => a -> a -> a
(+)

		findGeometricMean :: AGM -> GeometricMean
		findGeometricMean :: AGM -> ArithmeticMean
findGeometricMean	= squareRootAlgorithm
-> DecimalDigits -> ArithmeticMean -> ArithmeticMean
forall algorithm operand.
(Algorithmic algorithm, Real operand, Show operand) =>
algorithm -> DecimalDigits -> operand -> ArithmeticMean
Math.SquareRoot.squareRoot squareRootAlgorithm
squareRootAlgorithm DecimalDigits
decimalDigits (ArithmeticMean -> ArithmeticMean)
-> (AGM -> ArithmeticMean) -> AGM -> ArithmeticMean
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (ArithmeticMean -> ArithmeticMean -> ArithmeticMean)
-> AGM -> ArithmeticMean
forall a b c. (a -> b -> c) -> (a, b) -> c
uncurry ArithmeticMean -> ArithmeticMean -> ArithmeticMean
forall a. Num a => a -> a -> a
(*)
	in (AGM -> AGM) -> AGM -> [AGM]
forall a. (a -> a) -> a -> [a]
iterate (
		Strategy AGM -> AGM -> AGM
forall a. Strategy a -> a -> a
Control.Parallel.Strategies.withStrategy (
			Strategy ArithmeticMean -> Strategy ArithmeticMean -> Strategy AGM
forall a b. Strategy a -> Strategy b -> Strategy (a, b)
Control.Parallel.Strategies.parTuple2 Strategy ArithmeticMean
forall a. NFData a => Strategy a
Control.Parallel.Strategies.rdeepseq Strategy ArithmeticMean
forall a. NFData a => Strategy a
Control.Parallel.Strategies.rdeepseq
		) (AGM -> AGM) -> (AGM -> AGM) -> AGM -> AGM
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (ArithmeticMean -> ArithmeticMean
simplify (ArithmeticMean -> ArithmeticMean)
-> (AGM -> ArithmeticMean) -> AGM -> ArithmeticMean
forall b c a. (b -> c) -> (a -> b) -> a -> c
. AGM -> ArithmeticMean
findArithmeticMean (AGM -> ArithmeticMean) -> (AGM -> ArithmeticMean) -> AGM -> AGM
forall (a :: * -> * -> *) b c c'.
Arrow a =>
a b c -> a b c' -> a b (c, c')
&&& ArithmeticMean -> ArithmeticMean
simplify (ArithmeticMean -> ArithmeticMean)
-> (AGM -> ArithmeticMean) -> AGM -> ArithmeticMean
forall b c a. (b -> c) -> (a -> b) -> a -> c
. AGM -> ArithmeticMean
findGeometricMean)
	) AGM
agm

-- | Returns the bounds within which the 'AGM' has been constrained.
spread :: AGM -> Rational
spread :: AGM -> ArithmeticMean
spread	= (ArithmeticMean -> ArithmeticMean -> ArithmeticMean)
-> AGM -> ArithmeticMean
forall a b c. (a -> b -> c) -> (a, b) -> c
uncurry (-)

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