module Factory.Math.Implementations.SquareRoot(
Terms,
Algorithm(..)
) where
import Control.Arrow((***))
import qualified Data.Default
import Factory.Data.PrimeFactors((>/<), (>^))
import qualified Factory.Data.PrimeFactors as Data.PrimeFactors
import qualified Factory.Math.Implementations.Factorial as Math.Implementations.Factorial
import qualified Factory.Math.Power as Math.Power
import qualified Factory.Math.Precision as Math.Precision
import qualified Factory.Math.SquareRoot as Math.SquareRoot
import qualified Factory.Math.Summation as Math.Summation
type Terms = Int
data Algorithm
= BakhshaliApproximation
| ContinuedFraction
| HalleysMethod
| NewtonRaphsonIteration
| TaylorSeries Terms
deriving (Algorithm -> Algorithm -> Bool
(Algorithm -> Algorithm -> Bool)
-> (Algorithm -> Algorithm -> Bool) -> Eq Algorithm
forall a. (a -> a -> Bool) -> (a -> a -> Bool) -> Eq a
/= :: Algorithm -> Algorithm -> Bool
$c/= :: Algorithm -> Algorithm -> Bool
== :: Algorithm -> Algorithm -> Bool
$c== :: Algorithm -> Algorithm -> Bool
Eq, ReadPrec [Algorithm]
ReadPrec Algorithm
Int -> ReadS Algorithm
ReadS [Algorithm]
(Int -> ReadS Algorithm)
-> ReadS [Algorithm]
-> ReadPrec Algorithm
-> ReadPrec [Algorithm]
-> Read Algorithm
forall a.
(Int -> ReadS a)
-> ReadS [a] -> ReadPrec a -> ReadPrec [a] -> Read a
readListPrec :: ReadPrec [Algorithm]
$creadListPrec :: ReadPrec [Algorithm]
readPrec :: ReadPrec Algorithm
$creadPrec :: ReadPrec Algorithm
readList :: ReadS [Algorithm]
$creadList :: ReadS [Algorithm]
readsPrec :: Int -> ReadS Algorithm
$creadsPrec :: Int -> ReadS Algorithm
Read, Int -> Algorithm -> ShowS
[Algorithm] -> ShowS
Algorithm -> String
(Int -> Algorithm -> ShowS)
-> (Algorithm -> String)
-> ([Algorithm] -> ShowS)
-> Show Algorithm
forall a.
(Int -> a -> ShowS) -> (a -> String) -> ([a] -> ShowS) -> Show a
showList :: [Algorithm] -> ShowS
$cshowList :: [Algorithm] -> ShowS
show :: Algorithm -> String
$cshow :: Algorithm -> String
showsPrec :: Int -> Algorithm -> ShowS
$cshowsPrec :: Int -> Algorithm -> ShowS
Show)
instance Data.Default.Default Algorithm where
def :: Algorithm
def = Algorithm
NewtonRaphsonIteration
type ProblemSpecification operand
= Math.SquareRoot.Estimate
-> Math.Precision.DecimalDigits
-> operand
-> Math.SquareRoot.Result
instance Math.SquareRoot.Algorithmic Algorithm where
squareRootFrom :: Algorithm -> Estimate -> Int -> operand -> Result
squareRootFrom Algorithm
_ Estimate
_ Int
_ operand
0 = Result
0
squareRootFrom Algorithm
_ Estimate
_ Int
_ operand
1 = Result
1
squareRootFrom Algorithm
algorithm estimate :: Estimate
estimate@(Result
x, Int
decimalDigits) Int
requiredDecimalDigits operand
y
| Int
decimalDigits Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
>= Int
requiredDecimalDigits = Result
x
| Int
requiredDecimalDigits Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
<= Int
0 = String -> Result
forall a. HasCallStack => String -> a
error (String -> Result) -> String -> Result
forall a b. (a -> b) -> a -> b
$ String
"Factory.Math.Implementations.SquareRoot.squareRootFrom:\tinvalid number of required decimal digits; " String -> ShowS
forall a. [a] -> [a] -> [a]
++ Int -> String
forall a. Show a => a -> String
show Int
requiredDecimalDigits
| operand
y operand -> operand -> Bool
forall a. Ord a => a -> a -> Bool
< operand
0 = String -> Result
forall a. HasCallStack => String -> a
error (String -> Result) -> String -> Result
forall a b. (a -> b) -> a -> b
$ String
"Factory.Math.Implementations.SquareRoot.squareRootFrom:\tthere's no real square-root of " String -> ShowS
forall a. [a] -> [a] -> [a]
++ operand -> String
forall a. Show a => a -> String
show operand
y
| Bool
otherwise = (
case Algorithm
algorithm of
Algorithm
ContinuedFraction -> Estimate -> Int -> operand -> Result
forall operand. Real operand => ProblemSpecification operand
squareRootByContinuedFraction
Algorithm
_ -> Algorithm -> Estimate -> Int -> operand -> Result
forall operand.
Real operand =>
Algorithm -> ProblemSpecification operand
squareRootByIteration Algorithm
algorithm
) Estimate
estimate Int
requiredDecimalDigits operand
y
instance Math.SquareRoot.Iterator Algorithm where
step :: Algorithm -> operand -> Result -> Result
step Algorithm
BakhshaliApproximation operand
y Result
x
| Result
dy Result -> Result -> Bool
forall a. Eq a => a -> a -> Bool
== Result
0 = Result
x
| Bool
otherwise = Result
x' Result -> Result -> Result
forall a. Num a => a -> a -> a
- Result
dx'
where
dy, dydx, dx, x', dydx', dx' :: Math.SquareRoot.Result
dy :: Result
dy = operand -> Result -> Result
forall operand. Real operand => operand -> Result -> Result
Math.SquareRoot.getDiscrepancy operand
y Result
x
dydx :: Result
dydx = Result
2 Result -> Result -> Result
forall a. Num a => a -> a -> a
* Result
x
dx :: Result
dx = Result
dy Result -> Result -> Result
forall a. Fractional a => a -> a -> a
/ Result
dydx
x' :: Result
x' = Result
x Result -> Result -> Result
forall a. Num a => a -> a -> a
+ Result
dx
dydx' :: Result
dydx' = Result
2 Result -> Result -> Result
forall a. Num a => a -> a -> a
* Result
x'
dx' :: Result
dx' = Result -> Result
forall n. Num n => n -> n
Math.Power.square Result
dx Result -> Result -> Result
forall a. Fractional a => a -> a -> a
/ Result
dydx'
step Algorithm
HalleysMethod operand
y Result
x
| Result
dy Result -> Result -> Bool
forall a. Eq a => a -> a -> Bool
== Result
0 = Result
x
| Bool
otherwise = Result
x Result -> Result -> Result
forall a. Num a => a -> a -> a
- Result
dx
where
dy, dydx, dx :: Math.SquareRoot.Result
dy :: Result
dy = Result -> Result
forall n. Num n => n -> n
negate (Result -> Result) -> Result -> Result
forall a b. (a -> b) -> a -> b
$ operand -> Result -> Result
forall operand. Real operand => operand -> Result -> Result
Math.SquareRoot.getDiscrepancy operand
y Result
x
dydx :: Result
dydx = Result
2 Result -> Result -> Result
forall a. Num a => a -> a -> a
* Result
x
dx :: Result
dx = Result -> Result
forall a. Fractional a => a -> a
recip (Result -> Result) -> Result -> Result
forall a b. (a -> b) -> a -> b
$ Result
dydx Result -> Result -> Result
forall a. Fractional a => a -> a -> a
/ Result
dy Result -> Result -> Result
forall a. Num a => a -> a -> a
- Result -> Result
forall a. Fractional a => a -> a
recip Result
dydx
step Algorithm
NewtonRaphsonIteration operand
y Result
x = Result
x Result -> Result -> Result
forall a. Fractional a => a -> a -> a
/ Result
2 Result -> Result -> Result
forall a. Num a => a -> a -> a
+ (operand -> Result
forall a. Real a => a -> Result
toRational operand
y Result -> Result -> Result
forall a. Fractional a => a -> a -> a
/ Result
2) Result -> Result -> Result
forall a. Fractional a => a -> a -> a
/ Result
x
step (TaylorSeries Int
terms) operand
y Result
x = Int -> operand -> Result -> Result
forall operand. Real operand => Int -> operand -> Result -> Result
squareRootByTaylorSeries Int
terms operand
y Result
x
step Algorithm
algorithm operand
_ Result
_ = String -> Result
forall a. HasCallStack => String -> a
error (String -> Result) -> String -> Result
forall a b. (a -> b) -> a -> b
$ String
"Factory.Math.Implementations.SquareRoot.step:\tinappropriate algorithm; " String -> ShowS
forall a. [a] -> [a] -> [a]
++ Algorithm -> String
forall a. Show a => a -> String
show Algorithm
algorithm
convergenceOrder :: Algorithm -> Int
convergenceOrder Algorithm
BakhshaliApproximation = Int
Math.Precision.quarticConvergence
convergenceOrder Algorithm
ContinuedFraction = Int
Math.Precision.linearConvergence
convergenceOrder Algorithm
HalleysMethod = Int
Math.Precision.cubicConvergence
convergenceOrder Algorithm
NewtonRaphsonIteration = Int
Math.Precision.quadraticConvergence
convergenceOrder (TaylorSeries Int
terms) = Int
terms
squareRootByContinuedFraction :: Real operand => ProblemSpecification operand
squareRootByContinuedFraction :: ProblemSpecification operand
squareRootByContinuedFraction (Result
initialEstimate, Int
initialDecimalDigits) Int
requiredDecimalDigits operand
y = Result
initialEstimate Result -> Result -> Result
forall a. Num a => a -> a -> a
+ (Result -> [Result]
convergents Result
initialEstimate [Result] -> Int -> Result
forall a. [a] -> Int -> a
!! ConvergenceRate -> Int -> Int
forall i. Integral i => ConvergenceRate -> Int -> i
Math.Precision.getTermsRequired (ConvergenceRate
10 ConvergenceRate -> Int -> ConvergenceRate
forall a b. (Fractional a, Integral b) => a -> b -> a
^^ Int -> Int
forall n. Num n => n -> n
negate Int
initialDecimalDigits) Int
requiredDecimalDigits) where
convergents :: Math.SquareRoot.Result -> [Math.SquareRoot.Result]
convergents :: Result -> [Result]
convergents Result
x = (Result -> Result) -> Result -> [Result]
forall a. (a -> a) -> a -> [a]
iterate ((operand -> Result -> Result
forall operand. Real operand => operand -> Result -> Result
Math.SquareRoot.getDiscrepancy operand
y Result
x Result -> Result -> Result
forall a. Fractional a => a -> a -> a
/) (Result -> Result) -> (Result -> Result) -> Result -> Result
forall b c a. (b -> c) -> (a -> b) -> a -> c
. ((Result
2 Result -> Result -> Result
forall a. Num a => a -> a -> a
* Result
x) Result -> Result -> Result
forall a. Num a => a -> a -> a
+)) Result
0
taylorSeriesCoefficients :: Fractional f => [f]
taylorSeriesCoefficients :: [f]
taylorSeriesCoefficients = (Integer -> Integer -> f) -> [Integer] -> [Integer] -> [f]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith (
\Integer
powers Integer
n -> let
doubleN :: Integer
doubleN = Integer
2 Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
* Integer
n
product' :: Factors Integer Integer -> Integer
product' = BisectionRatio -> Int -> Factors Integer Integer -> Integer
forall base exponent.
(Num base, Integral exponent) =>
BisectionRatio -> Int -> Factors base exponent -> base
Data.PrimeFactors.product' (BisectionRatio -> BisectionRatio
forall a. Fractional a => a -> a
recip BisectionRatio
2) Int
10
in (f -> f -> f) -> (f, f) -> f
forall a b c. (a -> b -> c) -> (a, b) -> c
uncurry f -> f -> f
forall a. Fractional a => a -> a -> a
(/) ((f, f) -> f)
-> ((Factors Integer Integer, Factors Integer Integer) -> (f, f))
-> (Factors Integer Integer, Factors Integer Integer)
-> f
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (
Integer -> f
forall a b. (Integral a, Num b) => a -> b
fromIntegral (Integer -> f)
-> (Factors Integer Integer -> Integer)
-> Factors Integer Integer
-> f
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Factors Integer Integer -> Integer
product' (Factors Integer Integer -> f)
-> (Factors Integer Integer -> f)
-> (Factors Integer Integer, Factors Integer Integer)
-> (f, f)
forall (a :: * -> * -> *) b c b' c'.
Arrow a =>
a b c -> a b' c' -> a (b, b') (c, c')
*** Integer -> f
forall a b. (Integral a, Num b) => a -> b
fromIntegral (Integer -> f)
-> (Factors Integer Integer -> Integer)
-> Factors Integer Integer
-> f
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
* ((Integer
1 Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
- Integer
doubleN) Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
* Integer
powers)) (Integer -> Integer)
-> (Factors Integer Integer -> Integer)
-> Factors Integer Integer
-> Integer
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Factors Integer Integer -> Integer
product'
) ((Factors Integer Integer, Factors Integer Integer) -> f)
-> (Factors Integer Integer, Factors Integer Integer) -> f
forall a b. (a -> b) -> a -> b
$ Integer -> Factors Integer Integer
forall base. Integral base => base -> Factors base base
Math.Implementations.Factorial.primeFactors Integer
doubleN Factors Integer Integer
-> Factors Integer Integer
-> (Factors Integer Integer, Factors Integer Integer)
forall base exponent.
(Integral base, Integral exponent) =>
Factors base exponent
-> Factors base exponent
-> (Factors base exponent, Factors base exponent)
>/< Integer -> Factors Integer Integer
forall base. Integral base => base -> Factors base base
Math.Implementations.Factorial.primeFactors Integer
n Factors Integer Integer -> Integer -> Factors Integer Integer
forall exponent base.
Num exponent =>
Factors base exponent -> exponent -> Factors base exponent
>^ Integer
2
) (
(Integer -> Integer) -> Integer -> [Integer]
forall a. (a -> a) -> a -> [a]
iterate (Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
* Integer -> Integer
forall n. Num n => n -> n
negate Integer
4) Integer
1
) [Integer
0 :: Integer ..]
squareRootByTaylorSeries :: Real operand
=> Terms
-> operand
-> Math.SquareRoot.Result
-> Math.SquareRoot.Result
squareRootByTaylorSeries :: Int -> operand -> Result -> Result
squareRootByTaylorSeries Int
_ operand
_ Result
0 = String -> Result
forall a. HasCallStack => String -> a
error String
"Factory.Math.Implementations.SquareRoot.squareRootByTaylorSeries:\talgorithm can't cope with estimated value of zero."
squareRootByTaylorSeries Int
terms operand
y Result
x
| Int
terms Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
< Int
2 = String -> Result
forall a. HasCallStack => String -> a
error (String -> Result) -> String -> Result
forall a b. (a -> b) -> a -> b
$ String
"Factory.Math.Implementations.SquareRoot.squareRootByTaylorSeries:\tinvalid number of terms; " String -> ShowS
forall a. [a] -> [a] -> [a]
++ Int -> String
forall a. Show a => a -> String
show Int
terms
| Bool
otherwise = [Result] -> Result
forall i. Integral i => [Ratio i] -> Ratio i
Math.Summation.sumR' ([Result] -> Result)
-> ([Result] -> [Result]) -> [Result] -> Result
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Int -> [Result] -> [Result]
forall a. Int -> [a] -> [a]
take Int
terms ([Result] -> [Result])
-> ([Result] -> [Result]) -> [Result] -> [Result]
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (Result -> Result -> Result) -> [Result] -> [Result] -> [Result]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith Result -> Result -> Result
forall a. Num a => a -> a -> a
(*) [Result]
forall f. Fractional f => [f]
taylorSeriesCoefficients ([Result] -> Result) -> [Result] -> Result
forall a b. (a -> b) -> a -> b
$ (Result -> Result) -> Result -> [Result]
forall a. (a -> a) -> a -> [a]
iterate (Result -> Result -> Result
forall a. Num a => a -> a -> a
* Result
relativeError) Result
x
where
relativeError :: Math.SquareRoot.Result
relativeError :: Result
relativeError = Result -> Result
forall a. Enum a => a -> a
pred (Result -> Result) -> Result -> Result
forall a b. (a -> b) -> a -> b
$ operand -> Result
forall a. Real a => a -> Result
toRational operand
y Result -> Result -> Result
forall a. Fractional a => a -> a -> a
/ Result -> Result
forall n. Num n => n -> n
Math.Power.square Result
x
squareRootByIteration :: Real operand => Algorithm -> ProblemSpecification operand
squareRootByIteration :: Algorithm -> ProblemSpecification operand
squareRootByIteration Algorithm
algorithm (Result
initialEstimate, Int
initialDecimalDigits) Int
requiredDecimalDigits operand
y = (Result -> Result) -> Result -> [Result]
forall a. (a -> a) -> a -> [a]
iterate (Algorithm -> operand -> Result -> Result
forall algorithm operand.
(Iterator algorithm, Real operand) =>
algorithm -> operand -> Result -> Result
Math.SquareRoot.step Algorithm
algorithm operand
y) Result
initialEstimate [Result] -> Int -> Result
forall a. [a] -> Int -> a
!! Int -> Int -> Int -> Int
forall i. Integral i => Int -> Int -> Int -> i
Math.Precision.getIterationsRequired (Algorithm -> Int
forall algorithm. Iterator algorithm => algorithm -> Int
Math.SquareRoot.convergenceOrder Algorithm
algorithm) Int
initialDecimalDigits Int
requiredDecimalDigits