{-
	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@]	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			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

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

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

instance ToolShed.Defaultable.Defaultable Algorithm	where
	defaultValue	= 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 + realToFrac y / x) / 2		--This is identical to the /Babylonian Method/.
--	step NewtonRaphsonIteration y x	= x / 2 + realToFrac y / (2 * x)	--Faster.
	step NewtonRaphsonIteration y x	= x / 2 + (realToFrac 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 NewtonRaphsonIteration	= Math.Precision.quadraticConvergence
	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;
	<http://en.wikipedia.org/wiki/Solving_quadratic_equations_with_continued_fractions>,
	<http://en.wikipedia.org/wiki/Generalized_continued_fraction#Roots_of_positive_numbers>,
	<http://en.wikipedia.org/wiki/Methods_of_computing_square_roots#Continued_fraction_expansion>.
	<http://www.myreckonings.com/Dead_Reckoning/Online/Materials/General%20Method%20for%20Extracting%20Roots.pdf>

	* The convergence <http://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/; <http://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.

	* <http://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
>		=====		===========
>		2 terms		/quadratic/
>		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 $ realToFrac 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