module Factory.Math.Implementations.SquareRoot(
Terms,
Algorithm(..)
) where
import Control.Arrow((***))
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
import qualified ToolShed.Defaultable
type Terms = Int
data Algorithm
= BakhshaliApproximation
| ContinuedFraction
| HalleysMethod
| NewtonRaphsonIteration
| TaylorSeries Terms
deriving (Eq, Read, Show)
instance ToolShed.Defaultable.Defaultable Algorithm where
defaultValue = NewtonRaphsonIteration
type ProblemSpecification operand
= Math.SquareRoot.Estimate
-> Math.Precision.DecimalDigits
-> operand
-> Math.SquareRoot.Result
instance Math.SquareRoot.Algorithmic Algorithm where
squareRootFrom _ _ _ 0 = 0
squareRootFrom _ _ _ 1 = 1
squareRootFrom algorithm estimate@(x, decimalDigits) requiredDecimalDigits y
| decimalDigits >= requiredDecimalDigits = x
| requiredDecimalDigits <= 0 = error $ "Factory.Math.Implementations.SquareRoot.squareRootFrom:\tinvalid number of required decimal digits; " ++ show requiredDecimalDigits
| y < 0 = error $ "Factory.Math.Implementations.SquareRoot.squareRootFrom:\tthere's no real square-root of " ++ show y
| otherwise = (
case algorithm of
ContinuedFraction -> squareRootByContinuedFraction
_ -> squareRootByIteration algorithm
) estimate requiredDecimalDigits y
instance Math.SquareRoot.Iterator Algorithm where
step BakhshaliApproximation y x
| dy == 0 = x
| otherwise = x' dx'
where
dy, dydx, dx, x', dydx', dx' :: Math.SquareRoot.Result
dy = Math.SquareRoot.getDiscrepancy y x
dydx = 2 * x
dx = dy / dydx
x' = x + dx
dydx' = 2 * x'
dx' = Math.Power.square dx / dydx'
step HalleysMethod y x
| dy == 0 = x
| otherwise = x dx
where
dy, dydx, dx :: Math.SquareRoot.Result
dy = negate $ Math.SquareRoot.getDiscrepancy y x
dydx = 2 * x
dx = recip $ dydx / dy recip dydx
step NewtonRaphsonIteration y x = x / 2 + (toRational y / 2) / x
step (TaylorSeries terms) y x = squareRootByTaylorSeries terms y x
step algorithm _ _ = error $ "Factory.Math.Implementations.SquareRoot.step:\tinappropriate algorithm; " ++ show algorithm
convergenceOrder BakhshaliApproximation = Math.Precision.quarticConvergence
convergenceOrder ContinuedFraction = Math.Precision.linearConvergence
convergenceOrder HalleysMethod = Math.Precision.cubicConvergence
convergenceOrder NewtonRaphsonIteration = Math.Precision.quadraticConvergence
convergenceOrder (TaylorSeries terms) = terms
squareRootByContinuedFraction :: Real operand => ProblemSpecification operand
squareRootByContinuedFraction (initialEstimate, initialDecimalDigits) requiredDecimalDigits y = initialEstimate + (convergents initialEstimate !! Math.Precision.getTermsRequired (10 ^^ negate initialDecimalDigits) requiredDecimalDigits) where
convergents :: Math.SquareRoot.Result -> [Math.SquareRoot.Result]
convergents x = iterate ((Math.SquareRoot.getDiscrepancy y x /) . ((2 * x) +)) 0
taylorSeriesCoefficients :: Fractional f => [f]
taylorSeriesCoefficients = zipWith (
\powers n -> let
doubleN = 2 * n
product' = Data.PrimeFactors.product' (recip 2) 10
in uncurry (/) . (
fromIntegral . product' *** fromIntegral . (* ((1 doubleN) * powers)) . product'
) $ Math.Implementations.Factorial.primeFactors doubleN >/< Math.Implementations.Factorial.primeFactors n >^ 2
) (
iterate (* negate 4) 1
) [0 :: Integer ..]
squareRootByTaylorSeries :: Real operand
=> Terms
-> operand
-> Math.SquareRoot.Result
-> Math.SquareRoot.Result
squareRootByTaylorSeries _ _ 0 = error "Factory.Math.Implementations.SquareRoot.squareRootByTaylorSeries:\talgorithm can't cope with estimated value of zero."
squareRootByTaylorSeries terms y x
| terms < 2 = error $ "Factory.Math.Implementations.SquareRoot.squareRootByTaylorSeries:\tinvalid number of terms; " ++ show terms
| otherwise = Math.Summation.sumR' . take terms . zipWith (*) taylorSeriesCoefficients $ iterate (* relativeError) x
where
relativeError :: Math.SquareRoot.Result
relativeError = pred $ toRational y / Math.Power.square x
squareRootByIteration :: Real operand => Algorithm -> ProblemSpecification operand
squareRootByIteration algorithm (initialEstimate, initialDecimalDigits) requiredDecimalDigits y = iterate (Math.SquareRoot.step algorithm y) initialEstimate !! Math.Precision.getIterationsRequired (Math.SquareRoot.convergenceOrder algorithm) initialDecimalDigits requiredDecimalDigits