{-
	Copyright (C) 2011-2015 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 <https://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>.
	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

-- | 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 :: 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		-- The estimate was precise.
		| Bool
otherwise	= Result
x' Result -> Result -> Result
forall a. Num a => a -> a -> a
- Result
dx'	-- Correct the estimate.
		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	-- Identical to Newton-Raphson iteration.
			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'

{-
	* /Halley's/ method; <https://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 Algorithm
HalleysMethod operand
y Result
x
		| Result
dy Result -> Result -> Bool
forall a. Eq a => a -> a -> Bool
== Result
0	= Result
x		-- The estimate was precise.
		| Bool
otherwise	= Result
x Result -> Result -> Result
forall a. Num a => a -> a -> a
- Result
dx	-- Correct the estimate.
		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	-- Use the estimate to determine the error in 'y'.
			dydx :: Result
dydx	= Result
2 Result -> Result -> Result
forall a. Num a => a -> a -> a
* Result
x						-- The gradient, at the estimated value '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 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 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	-- Faster still.

	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	-- 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/Solving_quadratic_equations_with_continued_fractions>,
	<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>.
	<https://www.myreckonings.com/Dead_Reckoning/Online/Materials/General%20Method%20for%20Extracting%20Roots.pdf>

	* 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 :: 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

{- |
	* 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 :: [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) {-arbitrary-} Int
10 {-arbitrary-}
	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	-- (-4)^n
 ) [Integer
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
>		=====		===========
>		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 :: 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	-- 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 -> 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