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

	* Exports a common interface for /square-root/ implementations.

	* Provides utilities for these implementations.
-}

module Factory.Math.SquareRoot(
-- * Type-classes
	Algorithmic(..),
	Iterator(..),
-- * Types
-- ** Type-synonyms
	Result,
	Estimate,
-- * Functions
	getAccuracy,
	getDiscrepancy,
	getEstimate,
--	rSqrt,
-- ** Predicates
	isPrecise
) where

import qualified	Factory.Math.Power	as Math.Power
import qualified	Factory.Math.Precision	as Math.Precision

-- | The result-type; actually, only the concrete return-type of 'Math.Precision.simplify', stops it being a polymorphic instance of 'Fractional'.
type Result	= Rational

-- | Contains an estimate for the /square-root/ of a value, and its accuracy.
type Estimate	= (Result, Math.Precision.DecimalDigits)

-- | Defines the methods expected of a /square-root/ algorithm.
class Algorithmic algorithm	where
	squareRootFrom	:: (Real operand, Show operand)
		=> algorithm
		-> Estimate			-- ^ An initial estimate from which to start.
		-> Math.Precision.DecimalDigits	-- ^ The required precision.
		-> operand			-- ^ The value for which to find the /square-root/.
		-> Result			-- ^ Returns an improved estimate of the /square-root/, found using the specified algorithm, accurate to at least the required number of decimal digits.

	squareRoot	:: (Real operand, Show operand)
		=> algorithm
		-> Math.Precision.DecimalDigits	-- ^ The required precision.
		-> operand			-- ^ The value for which to find the /square-root/.
		-> Result			-- ^ Returns an estimate of the /square-root/, found using the specified algorithm, accurate to at least the required number of decimal digits.
	squareRoot algorithm decimalDigits operand	= squareRootFrom algorithm (getEstimate operand) decimalDigits operand	-- Default implementation

-- | The interface required to iterate, from an estimate of the required value, to the next approximation.
class Iterator algorithm where
	step :: Real operand
		=> algorithm
		-> operand	-- ^ The value for which the /square-root/ is required; @y@.
		-> Result	-- ^ The current estimate; @x(n)@.
		-> Result	-- ^ An improved estimate; @x(n+1)@.

	convergenceOrder :: algorithm -> Math.Precision.ConvergenceOrder	-- ^ The ultimate ratio of successive terms as the iteration converges.

-- | Generalise 'sqrt' to operate on any 'Real' operand.
rSqrt :: Real operand => operand -> Double
rSqrt	= sqrt . realToFrac

-- | Uses 'Double'-precision floating-point arithmetic, to obtain an initial estimate for the /square-root/, and its accuracy.
getEstimate :: (Real operand, Show operand) => operand -> Estimate
getEstimate y
	| y < 0		= error $ "Factory.Math.SquareRoot.getEstimate:\tthere's no real square-root of " ++ show y
	| otherwise	= (Math.Precision.simplify decimalDigits {-doubles performance by roughly length of the Rational representation-} . toRational $ rSqrt y, decimalDigits)
	where
		decimalDigits :: Math.Precision.DecimalDigits
		decimalDigits	= 16	-- <http://en.wikipedia.org/wiki/IEEE_floating_point>.

{- |
	* The signed difference between the /square/ of an estimate for the /square-root/ of a value, and that value.

	* Positive when the estimate is too low.

	* CAVEAT: the magnitude is twice the error in the /square-root/.
-}
getDiscrepancy :: Real operand => operand -> Result -> Result
getDiscrepancy y x	= toRational y - Math.Power.square x

-- | True if the specified estimate for the /square-root/, is precise.
isPrecise :: Real operand => operand -> Result -> Bool
isPrecise y x	= getDiscrepancy y x == 0

{- |
	* For a given value and an estimate of its /square-root/,
	returns the number of decimals digits to which the /square-root/ is accurate; including the integral digits.

	* CAVEAT: the result returned for an exact match has been bodged.
-}
getAccuracy :: Real operand => operand -> Result -> Math.Precision.DecimalDigits
getAccuracy y x
	| absoluteError == 0	= maxBound	-- Bodge.
--	| otherwise		= length . takeWhile (< 1) $ iterate (* 10) relativeError	-- CAVEAT: too slow.
	| otherwise		= length $ show (round $ toRational y / absoluteError :: Integer)
	where
		absoluteError :: Result
		absoluteError	= abs (getDiscrepancy y x) / 2	-- NB: the magnitude of the error in 'y', is twice the error in its square-root, 'x'.