{- 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'.