{-
	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@]	Defines the unit with which precision is measured, and operations on it.
-}
module Factory.Math.Precision(
-- * Types
-- ** Type-synonyms
	ConvergenceOrder,
	ConvergenceRate,
	DecimalDigits,
-- * Constants
	linearConvergence,
	quadraticConvergence,
	cubicConvergence,
	quarticConvergence,
-- * Functions
	getIterationsRequired,
	getTermsRequired,
	roundTo,
	promote,
	simplify
) where

import qualified	Data.Ratio

-- | The /order of convergence/; <https://en.wikipedia.org/wiki/Rate_of_convergence>.
type ConvergenceOrder	= Int

-- | The /rate of convergence/; <https://en.wikipedia.org/wiki/Rate_of_convergence>.
type ConvergenceRate	= Double

-- | A number of decimal digits; presumably positive.
type DecimalDigits	= Int

-- | /Linear/ convergence-rate; which may be qualified by the /rate of convergence/.
linearConvergence :: ConvergenceOrder
linearConvergence :: ConvergenceOrder
linearConvergence	= ConvergenceOrder
1

-- | /Quadratic/ convergence-rate.
quadraticConvergence :: ConvergenceOrder
quadraticConvergence :: ConvergenceOrder
quadraticConvergence	= ConvergenceOrder
2

-- | /Cubic/ convergence-rate.
cubicConvergence :: ConvergenceOrder
cubicConvergence :: ConvergenceOrder
cubicConvergence	= ConvergenceOrder
3

-- | /Quartic/ convergence-rate.
quarticConvergence :: ConvergenceOrder
quarticConvergence :: ConvergenceOrder
quarticConvergence	= ConvergenceOrder
4

-- | The predicted number of iterations, required to achieve a specific accuracy, at a given /order of convergence/.
getIterationsRequired :: Integral i
	=> ConvergenceOrder
	-> DecimalDigits	-- ^ The precision of the initial estimate.
	-> DecimalDigits	-- ^ The required precision.
	-> i
getIterationsRequired :: ConvergenceOrder -> ConvergenceOrder -> ConvergenceOrder -> i
getIterationsRequired ConvergenceOrder
convergenceOrder ConvergenceOrder
initialDecimalDigits ConvergenceOrder
requiredDecimalDigits
	| ConvergenceOrder
initialDecimalDigits ConvergenceOrder -> ConvergenceOrder -> Bool
forall a. Ord a => a -> a -> Bool
<= ConvergenceOrder
0	= [Char] -> i
forall a. HasCallStack => [Char] -> a
error ([Char] -> i) -> [Char] -> i
forall a b. (a -> b) -> a -> b
$ [Char]
"Factory.Math.Precision.getIterationsRequired:\tinsufficient 'initialDecimalDigits'; " [Char] -> [Char] -> [Char]
forall a. [a] -> [a] -> [a]
++ ConvergenceOrder -> [Char]
forall a. Show a => a -> [Char]
show ConvergenceOrder
initialDecimalDigits
	| Double
precisionRatio Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
<= Double
1		= i
0
	| Bool
otherwise			= Double -> i
forall a b. (RealFrac a, Integral b) => a -> b
ceiling (Double -> i) -> Double -> i
forall a b. (a -> b) -> a -> b
$ ConvergenceOrder -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral ConvergenceOrder
convergenceOrder Double -> Double -> Double
forall a. Floating a => a -> a -> a
`logBase` Double
precisionRatio
	where
		precisionRatio :: Double
		precisionRatio :: Double
precisionRatio	= ConvergenceOrder -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral ConvergenceOrder
requiredDecimalDigits Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ ConvergenceOrder -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral ConvergenceOrder
initialDecimalDigits

{- |
	* The predicted number of terms which must be extracted from a series,
	if it is to converge to the required accuracy,
	at the specified linear /convergence-rate/.

	* The /convergence-rate/ of a series, is the error in the series after summation of @(n+1)th@ terms,
	divided by the error after only @n@ terms, as the latter tends to infinity.
	As such, for a /convergent/ series (in which the error get smaller with successive terms), it's value lies in the range @0 .. 1@.

	* <https://en.wikipedia.org/wiki/Rate_of_convergence>.
-}
getTermsRequired :: Integral i
	=> ConvergenceRate
	-> DecimalDigits	-- ^ The additional number of correct decimal digits.
	-> i
getTermsRequired :: Double -> ConvergenceOrder -> i
getTermsRequired Double
_ ConvergenceOrder
0		= i
0
getTermsRequired Double
convergenceRate ConvergenceOrder
requiredDecimalDigits
	| Double
convergenceRate Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
<= Double
0 Bool -> Bool -> Bool
|| Double
convergenceRate Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
>= Double
1	= [Char] -> i
forall a. HasCallStack => [Char] -> a
error ([Char] -> i) -> [Char] -> i
forall a b. (a -> b) -> a -> b
$ [Char]
"Factory.Math.Precision.getTermsRequired:\t(0 < convergence-rate < 1); " [Char] -> [Char] -> [Char]
forall a. [a] -> [a] -> [a]
++ Double -> [Char]
forall a. Show a => a -> [Char]
show Double
convergenceRate
	| ConvergenceOrder
requiredDecimalDigits ConvergenceOrder -> ConvergenceOrder -> Bool
forall a. Ord a => a -> a -> Bool
< ConvergenceOrder
0			= [Char] -> i
forall a. HasCallStack => [Char] -> a
error ([Char] -> i) -> [Char] -> i
forall a b. (a -> b) -> a -> b
$ [Char]
"Factory.Math.Precision.getTermsRequired:\t'requiredDecimalDigits' must be positive; " [Char] -> [Char] -> [Char]
forall a. [a] -> [a] -> [a]
++ ConvergenceOrder -> [Char]
forall a. Show a => a -> [Char]
show ConvergenceOrder
requiredDecimalDigits
	| Bool
otherwise					= Double -> i
forall a b. (RealFrac a, Integral b) => a -> b
ceiling (Double -> i) -> Double -> i
forall a b. (a -> b) -> a -> b
$ ConvergenceOrder -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral ConvergenceOrder
requiredDecimalDigits Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double -> Double
forall a. Num a => a -> a
negate (Double -> Double -> Double
forall a. Floating a => a -> a -> a
logBase Double
10 Double
convergenceRate)

-- | Rounds the specified number, to a positive number of 'DecimalDigits'.
roundTo :: (RealFrac a, Fractional f) => DecimalDigits -> a -> f
roundTo :: ConvergenceOrder -> a -> f
roundTo ConvergenceOrder
decimals = (f -> f -> f
forall a. Fractional a => a -> a -> a
/ Integer -> f
forall a. Num a => Integer -> a
fromInteger Integer
promotionFactor) (f -> f) -> (a -> f) -> a -> f
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Integer -> f
forall a. Num a => Integer -> a
fromInteger (Integer -> f) -> (a -> Integer) -> a -> f
forall b c a. (b -> c) -> (a -> b) -> a -> c
. a -> Integer
forall a b. (RealFrac a, Integral b) => a -> b
round (a -> Integer) -> (a -> a) -> a -> Integer
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (a -> a -> a
forall a. Num a => a -> a -> a
* Integer -> a
forall a. Num a => Integer -> a
fromInteger Integer
promotionFactor)	where
	promotionFactor :: Integer
	promotionFactor :: Integer
promotionFactor	= Integer
10 Integer -> ConvergenceOrder -> Integer
forall a b. (Num a, Integral b) => a -> b -> a
^ ConvergenceOrder
decimals

-- | Promotes the specified number, by a positive number of 'DecimalDigits'.
promote :: Num n => n -> DecimalDigits -> n
promote :: n -> ConvergenceOrder -> n
promote n
x	= (n -> n -> n
forall a. Num a => a -> a -> a
* n
x) (n -> n) -> (ConvergenceOrder -> n) -> ConvergenceOrder -> n
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (n
10 n -> ConvergenceOrder -> n
forall a b. (Num a, Integral b) => a -> b -> a
^)

{- |
	* Reduces a 'Rational' to the minimal form required for the specified number of /fractional/ decimal places;
	irrespective of the number of integral decimal places.

	* A 'Rational' approximation to an irrational number, may be very long, and provide an unknown excess precision.
	Whilst this doesn't sound harmful, it costs in performance and memory-requirement, and being unpredictable isn't actually useful.
-}
simplify :: RealFrac operand
	=> DecimalDigits	-- ^ The number of places after the decimal point, which are required.
	-> operand
	-> Rational
simplify :: ConvergenceOrder -> operand -> Rational
simplify ConvergenceOrder
decimalDigits operand
operand	= operand -> operand -> Rational
forall a. RealFrac a => a -> a -> Rational
Data.Ratio.approxRational operand
operand (operand -> Rational)
-> (operand -> operand) -> operand -> Rational
forall b c a. (b -> c) -> (a -> b) -> a -> c
. operand -> operand
forall a. Fractional a => a -> a
recip (operand -> Rational) -> operand -> Rational
forall a b. (a -> b) -> a -> b
$ operand
4 operand -> operand -> operand
forall a. Num a => a -> a -> a
* operand
10 operand -> ConvergenceOrder -> operand
forall a b. (Num a, Integral b) => a -> b -> a
^ ConvergenceOrder -> ConvergenceOrder
forall a. Enum a => a -> a
succ ConvergenceOrder
decimalDigits	-- Tolerate any error less than half the least significant digit required.