```{-
Copyright (C) 2011-2015 Dr. Alistair Ward

This program is free software: you can redistribute it and/or modify
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@]	Implements 'Math.SquareRoot.Algorithmic' by a variety of methods.

[@CAVEAT@]

Caller may benefit from application of 'Math.Precision.simplify' before operating on the result;
which though of the required accuracy, may not be the most concise rational number satisfying that criterion.
-}
module Factory.Math.Implementations.SquareRoot(
-- * Types
-- ** Type-synonyms
--	ProblemSpecification,
Terms,
-- ** Data-types
Algorithm(..)
-- * Functions
--	squareRootByContinuedFraction,
--	squareRootByIteration,
--	squareRootByTaylorSeries,
--	taylorSeriesCoefficients
) 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

-- | The number of terms in a series.
type Terms	= Int

-- | The algorithms by which the /square-root/ has been implemented.
data Algorithm
= BakhshaliApproximation	-- ^ <https://en.wikipedia.org/wiki/Methods_of_computing_square_roots#Bakhshali_approximation>
| ContinuedFraction		-- ^ <https://en.wikipedia.org/wiki/Methods_of_computing_square_roots#Continued_fraction_expansion>.
| HalleysMethod			-- ^ <https://en.wikipedia.org/wiki/Halley%27s_method>.
| NewtonRaphsonIteration	-- ^ <https://en.wikipedia.org/wiki/Newton%27s_method>.
| TaylorSeries Terms		-- ^ <https://en.wikipedia.org/wiki/Methods_of_computing_square_roots#Taylor_series>.

instance Data.Default.Default Algorithm	where
def	= NewtonRaphsonIteration

-- | Returns an improved estimate for the /square-root/ of the specified value, to the required precision, using the supplied initial estimate..
type ProblemSpecification operand
= Math.SquareRoot.Estimate
-> Math.Precision.DecimalDigits	-- ^ The required precision.
-> operand			-- ^ The value for which to find the /square-root/.
-> 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		-- The estimate was precise.
| otherwise	= x' - dx'	-- Correct the estimate.
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	-- Identical to Newton-Raphson iteration.
dydx'	= 2 * x'
dx'	= Math.Power.square dx / dydx'

{-
* /Halley's/ method; <http://mathworld.wolfram.com/HalleysMethod.html>

>		X(n+1) = Xn - f(Xn) / [f'(Xn) - f''(Xn) * f(Xn) / 2 * f'(Xn)]
>			=> Xn - (Xn^2 - Y) / [2Xn - 2 * (Xn^2 - Y) / 2 * 2Xn] where Y = X^2, f(X) = X^2 - Y, f'(X) = 2X, f''(X) = 2
>			=> Xn - 1 / [2Xn / (Xn^2 - Y) - 1 / 2Xn]
-}
step HalleysMethod y x
| dy == 0	= x		-- The estimate was precise.
| otherwise	= x - dx	-- Correct the estimate.
where
dy, dydx, dx :: Math.SquareRoot.Result
dy	= negate \$ Math.SquareRoot.getDiscrepancy y x	-- Use the estimate to determine the error in 'y'.
dydx	= 2 * x						-- The gradient, at the estimated value 'x'.
dx	= recip \$ dydx / dy - recip dydx

--	step NewtonRaphsonIteration y x	= (x + toRational y / x) / 2		-- This is identical to the /Babylonian Method/.
--	step NewtonRaphsonIteration y x	= x / 2 + toRational y / (2 * x)	-- Faster.
step NewtonRaphsonIteration y x	= x / 2 + (toRational y / 2) / x	-- Faster still.

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 (TaylorSeries terms)	= terms	-- The order of convergence, per iteration, equals the number of terms in the series on each iteration.

{- |
* Uses /continued-fractions/, to iterate towards the principal /square-root/ of the specified positive integer;
<https://en.wikipedia.org/wiki/Generalized_continued_fraction#Roots_of_positive_numbers>,
<https://en.wikipedia.org/wiki/Methods_of_computing_square_roots#Continued_fraction_expansion>.

* The convergence <https://en.wikipedia.org/wiki/Rate_of_convergence> of the /continued-fraction/ is merely /1st order/ (linear).
-}
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

{- |
* The constant coefficients of the /Taylor-series/ for a /square-root/; <https://en.wikipedia.org/wiki/Methods_of_computing_square_roots#Taylor_series>.

* @ ((-1)^n * factorial(2*n)) / ((1 - 2*n) * 4^n * factorial(n^2)) @.
-}
taylorSeriesCoefficients :: Fractional f => [f]
taylorSeriesCoefficients	= zipWith (
\powers n	-> let
doubleN		= 2 * n
product'	= Data.PrimeFactors.product' (recip 2) {-arbitrary-} 10 {-arbitrary-}
in uncurry (/) . (
fromIntegral . product' *** fromIntegral . (* ((1 - doubleN) * powers)) . product'
) \$ Math.Implementations.Factorial.primeFactors doubleN >/< Math.Implementations.Factorial.primeFactors n >^ 2
) (
iterate (* negate 4) 1	-- (-4)^n
) [0 :: Integer ..]		-- n

{- |
* Returns the /Taylor-series/ for the /square-root/ of the specified value, to any requested number of terms.

* <https://en.wikipedia.org/wiki/Methods_of_computing_square_roots#Taylor_series>.

* The convergence of the series is merely /linear/,
in that each term increases the precision, by a constant number of decimal places, equal to the those in the original estimate.

* By feeding-back the improved estimate, to form a new series, the order of convergence, on each successive iteration,
becomes proportional to the number of terms;

>		Terms		Convergence
>		=====		===========
>		3 terms		/cubic/
-}
squareRootByTaylorSeries :: Real operand
=> Terms			-- ^ The number of terms of the infinite series, to evaluate.
-> operand			-- ^ The value for which the /square-root/ is required.
-> Math.SquareRoot.Result	-- ^ An initial estimate.
-> 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	-- Pedantically, this is the error in y, which is twice the magnitude of the error in x.

-- | Iterates from the estimated value, towards the /square-root/, a sufficient number of times to achieve the required accuracy.
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

```