summaryrefslogtreecommitdiff
diff options
context:
space:
mode:
authorJamesCook <>2010-09-10 19:59:58 (GMT)
committerLuite Stegeman <luite@luite.com>2010-09-10 19:59:58 (GMT)
commitb245cd924549f5780a829de89c999da9ddce8b41 (patch)
treecf4fa8a0a3b7e3b982ab54969b03e8c314dad26c
version 0.10.1
-rw-r--r--Setup.lhs5
-rw-r--r--roots.cabal40
-rw-r--r--src/Math/Root/Bracket.hs46
-rw-r--r--src/Math/Root/Finder.hs88
-rw-r--r--src/Math/Root/Finder/Bisection.hs38
-rw-r--r--src/Math/Root/Finder/Brent.hs143
-rw-r--r--src/Math/Root/Finder/Dekker.hs53
-rw-r--r--src/Math/Root/Finder/FalsePosition.hs46
-rw-r--r--src/Math/Root/Finder/InverseQuadratic.hs40
-rw-r--r--src/Math/Root/Finder/Newton.hs33
-rw-r--r--src/Math/Root/Finder/Ridders.hs58
-rw-r--r--src/Math/Root/Finder/Secant.hs46
12 files changed, 636 insertions, 0 deletions
diff --git a/Setup.lhs b/Setup.lhs
new file mode 100644
index 0000000..8193653
--- /dev/null
+++ b/Setup.lhs
@@ -0,0 +1,5 @@
+#!/usr/bin/env runhaskell
+
+> import Distribution.Simple
+> main = defaultMain
+
diff --git a/roots.cabal b/roots.cabal
new file mode 100644
index 0000000..b9cf1a8
--- /dev/null
+++ b/roots.cabal
@@ -0,0 +1,40 @@
+name: roots
+version: 0.1
+stability: experimental
+
+cabal-version: >= 1.6
+build-type: Simple
+
+author: James Cook <mokus@deepbondi.net>
+maintainer: James Cook <mokus@deepbondi.net>
+license: PublicDomain
+homepage: /dev/null
+
+category: Math, Numerical
+synopsis: Root-finding algorithms (1-dimensional)
+description: Framework for and a few implementations of
+ (1-dimensional) numerical root-finding algorithms.
+
+tested-with: GHC == 6.8.3,
+ GHC == 6.10.4,
+ GHC == 6.12.1, GHC == 6.12.3
+
+source-repository head
+ type: darcs
+ location: http://code.haskell.org/~mokus/roots
+
+Library
+ ghc-options: -Wall
+ hs-source-dirs: src
+ exposed-modules: Math.Root.Bracket
+ Math.Root.Finder
+ Math.Root.Finder.Bisection
+ Math.Root.Finder.Brent
+ Math.Root.Finder.Dekker
+ Math.Root.Finder.FalsePosition
+ Math.Root.Finder.InverseQuadratic
+ Math.Root.Finder.Newton
+ Math.Root.Finder.Ridders
+ Math.Root.Finder.Secant
+
+ build-depends: base >= 3 && <5, tagged
diff --git a/src/Math/Root/Bracket.hs b/src/Math/Root/Bracket.hs
new file mode 100644
index 0000000..be7bc5b
--- /dev/null
+++ b/src/Math/Root/Bracket.hs
@@ -0,0 +1,46 @@
+{-# OPTIONS_GHC -fno-warn-name-shadowing #-}
+module Math.Root.Bracket where
+
+-- |Predicate that returns 'True' whenever the given pair of points brackets
+-- a root of the given function.
+brackets :: (Eq a, Num b) => (a -> b) -> (a,a) -> Bool
+brackets f (x1,x2)
+ | x1 == x2 = f x1 == 0
+ | otherwise = signum (f x1) /= signum (f x2)
+
+-- |@bracket f x1 x2@: Given a function and an initial guessed range x1 to x2,
+-- this function expands the range geometrically until a root is bracketed by
+-- the returned values, returning a list of the successively expanded ranges.
+-- The list will be finite if and only if the sequence yields a bracketing pair.
+bracket :: (Fractional a, Num b, Ord b) =>
+ (a -> b) -> a -> a -> [(a, a)]
+bracket f x1 x2
+ | x1 == x2 = error "bracket: empty range"
+ | otherwise = go x1 (f x1) x2 (f x2)
+ where
+ factor = 1.618 -- golden ratio (close enough to it, anyway)
+ go x1 f1 x2 f2
+ | signum f1 /= signum f2 = [(x1, x2)]
+ | abs f1 < abs f2 = (x1, x2) : go x1' (f x1') x2 f2
+ | otherwise = (x1, x2) : go x1 f1 x2' (f x2')
+ where
+ x1' = x1 - factor * w
+ x2' = x2 + factor * w
+ w = x2 - x1
+
+-- |@subdivideAndBracket f x1 x2 n@: Given a function defined on the interval
+-- [x1,x2], subdivide the interval into n equally spaced segments and search
+-- for zero crossings of the function. The returned list will contain all
+-- bracketing pairs found.
+subdivideAndBracket :: (Num b, Fractional a, Integral c) =>
+ (a -> b) -> a -> a -> c -> [(a, a)]
+subdivideAndBracket f x1 x2 n =
+ [ (x1, x2)
+ | ((x1, y1), (x2, y2)) <- zip xys (tail xys)
+ , signum y1 /= signum y2
+ ]
+ where
+ dx = (x2 - x1) / fromIntegral n
+ xs = x1 : [x1 + dx * fromIntegral i | i <- [1..n]]
+ xys = map (\x -> (x, f x)) xs
+
diff --git a/src/Math/Root/Finder.hs b/src/Math/Root/Finder.hs
new file mode 100644
index 0000000..ed5cdf6
--- /dev/null
+++ b/src/Math/Root/Finder.hs
@@ -0,0 +1,88 @@
+{-# LANGUAGE MultiParamTypeClasses, ScopedTypeVariables, FlexibleContexts #-}
+module Math.Root.Finder where
+
+import Control.Monad.Instances ()
+import Data.Tagged
+
+-- |General interface for numerical root finders.
+class RootFinder r a b where
+ -- |@initRootFinder f x0 x1@: Initialize a root finder for the given
+ -- function with the initial bracketing interval (x0,x1).
+ initRootFinder :: (a -> b) -> a -> a -> r a b
+
+ -- |Step a root finder for the given function (which should generally
+ -- be the same one passed to @initRootFinder@), refining the finder's
+ -- estimate of the location of a root.
+ stepRootFinder :: (a -> b) -> r a b -> r a b
+
+ -- |Extract the finder's current estimate of the position of a root.
+ estimateRoot :: r a b -> a
+
+ -- |Extract the finder's current estimate of the upper bound of the
+ -- distance from @estimateRoot@ to an actual root in the function.
+ --
+ -- Generally, @estimateRoot r@ +- @estimateError r@ should bracket
+ -- a root of the function.
+ estimateError :: r a b -> a
+
+ -- |Test whether a root finding algorithm has converged to a given
+ -- relative accuracy.
+ converged :: (Num a, Ord a) => a -> r a b -> Bool
+ converged xacc r = abs (estimateError r) <= abs xacc
+
+ -- |Default number of steps after which root finding will be deemed
+ -- to have failed. Purely a convenience used to control the behavior
+ -- of built-in functions such as 'findRoot' and 'traceRoot'. The
+ -- default value is 250.
+ defaultNSteps :: Tagged (r a b) Int
+ defaultNSteps = Tagged 250
+
+-- |@traceRoot f x0 x1 mbEps@ initializes a root finder and repeatedly
+-- steps it, returning each step of the process in a list. When the algorithm
+-- terminates or the 'defaultNSteps' limit is exceeded, the list ends.
+-- Termination criteria depends on @mbEps@; if it is of the form @Just eps@
+-- then convergence to @eps@ is used (using the @converged@ method of the
+-- root finder). Otherwise, the trace is not terminated until subsequent
+-- states are equal (according to '=='). This is a stricter condition than
+-- convergence to 0; subsequent states may have converged to zero but as long
+-- as any internal state changes the trace will continue.
+traceRoot :: (Eq (r a b), RootFinder r a b, Num a, Ord a) =>
+ (a -> b) -> a -> a -> Maybe a -> [r a b]
+traceRoot f a b xacc = go nSteps start (stepRootFinder f start)
+ where
+ Tagged nSteps = (const :: Tagged a b -> a -> Tagged a b) defaultNSteps start
+ start = initRootFinder f a b
+
+ -- lookahead 1; if tracing with no convergence test, apply a
+ -- naive test to bail out if the root stops changing. This is
+ -- provided because that's not always the same as convergence to 0,
+ -- and the main purpose of this function is to watch what actually
+ -- happens inside the root finder.
+ go n x next
+ | maybe (x==next) (flip converged x) xacc = [x]
+ | n <= 0 = []
+ | otherwise = x : go (n-1) next (stepRootFinder f next)
+
+-- |@findRoot f x0 x1 eps@ initializes a root finder and repeatedly
+-- steps it. When the algorithm converges to @eps@ or the 'defaultNSteps'
+-- limit is exceeded, the current best guess is returned, with the @Right@
+-- constructor indicating successful convergence or the @Left@ constructor
+-- indicating failure to converge.
+findRoot :: (RootFinder r a b, Num a, Ord a) =>
+ (a -> b) -> a -> a -> a -> Either (r a b) (r a b)
+findRoot f a b xacc = go nSteps start
+ where
+ Tagged nSteps = (const :: Tagged a b -> a -> Tagged a b) defaultNSteps start
+ start = initRootFinder f a b
+
+ go n x
+ | converged xacc x = Right x
+ | n <= 0 = Left x
+ | otherwise = go (n-1) (stepRootFinder f x)
+
+-- |A useful constant: 'eps' is (for most 'RealFloat' types) the smallest
+-- positive number such that @1 + eps /= 1@.
+eps :: RealFloat a => a
+eps = eps'
+ where
+ eps' = encodeFloat 1 (1 - floatDigits eps')
diff --git a/src/Math/Root/Finder/Bisection.hs b/src/Math/Root/Finder/Bisection.hs
new file mode 100644
index 0000000..a511ad1
--- /dev/null
+++ b/src/Math/Root/Finder/Bisection.hs
@@ -0,0 +1,38 @@
+{-# LANGUAGE MultiParamTypeClasses, FlexibleInstances #-}
+module Math.Root.Finder.Bisection
+ ( Bisect, bisection
+ ) where
+
+import Math.Root.Finder
+
+-- |Bisect an interval in search of a root. At all times, @f (estimateRoot _)@
+-- is less than or equal to 0 and @f (estimateRoot _ + estimateError _)@ is
+-- greater than or equal to 0.
+data Bisect a b = Bisect { _bisX :: !a, _bisF :: !b , _bisDX :: !a }
+ deriving (Eq, Ord, Show)
+
+instance (Fractional a, Ord b, Num b) => RootFinder Bisect a b where
+ initRootFinder f x1 x2
+ | f1 < 0 = Bisect x1 f1 (x2-x1)
+ | otherwise = Bisect x2 f2 (x1-x2)
+ where f1 = f x1; f2 = f x2
+
+ stepRootFinder f orig@(Bisect x fx dx) = case fMid `compare` 0 of
+ LT -> Bisect xMid fMid dx2
+ EQ -> orig
+ GT -> Bisect x fx dx2
+ where
+ dx2 = dx * 0.5
+ xMid = x + dx2
+ fMid = f xMid
+
+ estimateRoot (Bisect x _ _) = x
+
+ estimateError (Bisect _ 0 _) = 0
+ estimateError (Bisect _ _ dx) = dx
+
+-- |Using bisection, return a root of a function known to lie between x1 and x2.
+-- The root will be refined till its accuracy is +-xacc. If convergence fails,
+-- returns the final state of the search.
+bisection :: (Ord a, Fractional a, Ord b, Num b) => (a -> b) -> a -> a -> a -> Either (Bisect a b) a
+bisection f x1 x2 xacc = fmap estimateRoot (findRoot f x1 x2 xacc)
diff --git a/src/Math/Root/Finder/Brent.hs b/src/Math/Root/Finder/Brent.hs
new file mode 100644
index 0000000..0acd5a8
--- /dev/null
+++ b/src/Math/Root/Finder/Brent.hs
@@ -0,0 +1,143 @@
+{-# LANGUAGE MultiParamTypeClasses, FlexibleInstances, FlexibleContexts #-}
+{-# OPTIONS_GHC -fno-warn-name-shadowing #-}
+module Math.Root.Finder.Brent
+ ( Brent
+ , brent
+ ) where
+
+import Math.Root.Bracket
+import Math.Root.Finder
+import Data.Maybe
+import Text.Printf
+
+-- Invariants:
+-- 1) B and C bracket the root
+-- 2) |f(B)| <= |f(C)|
+-- 3) min(f(B),f(C)) <= f(A) <= max(f(B),f(C))
+-- 4) e >= 0
+-- |Working state for Brent's root-finding method.
+data Brent a b = Brent
+ { brA :: !a
+ , brFA :: !b
+ , brB :: !a
+ , brFB :: !b
+ , brC :: !a
+ , brFC :: !b
+ , brE :: a
+ } deriving (Eq, Show)
+
+-- TODO: clean up this mess!
+instance (RealFloat a, Real b, Fractional b) => RootFinder Brent a b where
+ initRootFinder f x1 x2 = fixMagnitudes (Brent x1 f1 x2 f2 x1 f1 dx)
+ where f1 = f x1; f2 = f x2; dx = x2 - x1
+
+ stepRootFinder f r@(Brent a fa b fb c fc e)
+ | e >= 2 * min tol1 abs_s -- require that the method be making progress, overall
+ && 1.5 * m * signum s >= tol1 + abs_s -- require that the proposed step is getting closer to 'b' - specifically, s should be between 0 and 0.75*(c - b)
+ = advance s abs_s
+ | otherwise = advance m (abs (b - a))
+ where
+ -- Minimum step size to continue with inverse-quadratic interpolation
+ tol1 = eps * (abs b + 0.5)
+ abs_s = abs s
+
+ -- midpoint for bisection step
+ m = 0.5 * (c - b)
+
+ -- subdivision point for inverse quadratic interpolation step
+ s | fa /= fc && fa /= fb
+ = let a' = realToFrac (fa / (fc - fb))
+ b' = realToFrac (fb / (fc - fa))
+ c' = realToFrac (fc / (fb - fa))
+ in (a' * b' * c) - ((a' * c' + 1) * b) + (a * b' * c')
+ | otherwise
+ -- Fall back to linear interpolation when quadratic
+ -- interpolation will yield nonsensical results.
+ = (c - b) * realToFrac (fb / (fb - fc))
+
+ -- |Moves the current estimate by 'd' (or by tol1, whichever
+ -- is greater) and sets 'brE' to 'e', maintaining all invariants.
+ -- Ensuring that at least some tiny jump is made allows quick
+ -- discovery and termination in the case where the current best
+ -- estimate is already nearly on top of the root. Without such
+ -- a check, the method would repeatedly tighten the 'c' bound
+ -- by bisection every other step, which is really rather stupid
+ -- if 'b' is already sitting on a root.
+ advance d newE = update b' (f b') newE r
+ where
+ b' = if abs d > tol1 then b + d else b + tol1 * signum m
+
+
+ estimateRoot = brB
+ estimateError = brE
+ converged _ Brent{brFB = 0} = True
+ converged tol Brent{brB = b, brE = e} =
+ abs e <= 4 * eps * abs b + tol
+
+-- |@brent f x1 x2 xacc@: attempt to find a root of a function known to
+-- lie between x1 and x2, using Brent's method. The root will be refined
+-- till its accuracy is +-xacc. If convergence fails, returns the final
+-- state of the search.
+brent :: RealFloat a => (a -> a) -> a -> a -> a -> Either (Brent a a) a
+brent f x1 x2 xacc = fmap estimateRoot (findRoot f x1 x2 xacc)
+
+-- |Updates the state by incorporating a new estimate and setting 'brE',
+-- maintaining all invariants.
+update :: (Num a, Num b, Ord b) => a -> b -> a -> Brent a b -> Brent a b
+update b fb e r@Brent{brB = a, brFB = fa}
+ = fixMagnitudes (fixSigns r{brA = a, brFA = fa, brB = b, brFB = fb, brE = e})
+
+-- Establish invariant (1) that b and c bracket the root,
+-- based on precondition that (a,c) already does.
+--
+-- (a,c) brackets implies that either (b,c) or (a,b) brackets. In the
+-- former case, nothing needs to be done as (by construction) either fb is already
+-- between fa and fc or b is already between a and c (depending which kind of
+-- step was taken). In the latter case, discard C and use A in its place, because
+-- c and fc are both (by the existing invariants - (a,c) bracket, |f(c)| >= |f(a)|)
+-- outside the new region of interest.
+fixSigns :: (Num a, Num b, Ord b) => Brent a b -> Brent a b
+fixSigns br@Brent{ brA = a
+ , brFA = fa, brFB = fb, brFC = fc }
+ | (fb > 0 && fc > 0) || (fb < 0 && fc < 0)
+ = br { brC = a, brFC = fa }
+ | otherwise
+ = br
+
+-- Establish invariant (2) that |f(c)| >= |f(b)| and invariant (3) that
+-- 'fa' falls between fb and fc.
+fixMagnitudes :: (Num b, Ord b) => Brent a b -> Brent a b
+fixMagnitudes br@Brent{ brC = c, brB = b
+ , brFC = fc, brFB = fb }
+ | abs fc < abs fb
+ = br { brA = b, brFA = fb
+ , brB = c, brFB = fc
+ , brC = b, brFC = fb
+ }
+ | otherwise
+ = br
+
+-- |debugging function to show a nice trace of the progress of the algorithm
+_traceBrent :: (PrintfArg a, RealFloat a,
+ PrintfArg b, Ord b, Num b,
+ RootFinder Brent a b) =>
+ (a -> b) -> Maybe (a, a) -> IO ()
+_traceBrent f mbRange = do
+ xs <- sequence
+ [ put br >> return br
+ | br <- traceRoot f x0 x1 (Just eps)
+ ]
+
+ putStrLn "(converged)"
+ go (last xs)
+ where
+ (x0,x1) = fromMaybe (last (bracket f 0 1)) mbRange
+ put Brent{brA=a, brB=b, brC=c, brFA=fa, brFB=fb, brFC=fc} =
+ putStrLn . map fixPlus $
+ printf (concat (replicate 6 "%-+25g")) a b c fa fb fc
+ fixPlus '+' = ' '
+ fixPlus c = c
+ go x
+ | x == x' = return ()
+ | otherwise = put x >> go x'
+ where x' = stepRootFinder f x \ No newline at end of file
diff --git a/src/Math/Root/Finder/Dekker.hs b/src/Math/Root/Finder/Dekker.hs
new file mode 100644
index 0000000..00061cb
--- /dev/null
+++ b/src/Math/Root/Finder/Dekker.hs
@@ -0,0 +1,53 @@
+{-# LANGUAGE MultiParamTypeClasses, FlexibleInstances #-}
+module Math.Root.Finder.Dekker (Dekker, dekker) where
+
+import Math.Root.Finder
+
+-- fields: a fa b fb oldB oldFb
+-- invariants:
+-- 1) signum fA /= signum fB
+-- 2) abs fB <= abs fA
+-- 3) (oldB-a)*(oldB-b) >= 0
+data Dekker a b = Dekker !a !b !a !b !a !b deriving (Eq, Show)
+
+-- |@dekker f x1 x2 xacc@: attempt to find a root of a function known to
+-- lie between x1 and x2, using Dekker's method. The root will be refined
+-- till its accuracy is +-xacc. If convergence fails, returns the final
+-- state of the search.
+dekker :: RealFloat a => (a -> a) -> a -> a -> a -> Either (Dekker a a) a
+dekker f x1 x2 xacc = fmap estimateRoot (findRoot f x1 x2 xacc)
+
+instance (Fractional a, Ord a, Real b, Fractional b, Ord b) => RootFinder Dekker a b where
+ initRootFinder f x0 x1
+ | signum f0 == signum f1 = error "initRootFinder/Dekker: starting points do not (obviously) bracket a root"
+ | abs f0 <= abs f1 = Dekker x1 f1 x0 f0 x1 f1
+ | otherwise = Dekker x0 f0 x1 f1 x0 f0
+ where f0 = f x0; f1 = f x1
+
+ stepRootFinder f orig@(Dekker a _ b fb oldB oldFb)
+ | fb == 0 = orig
+ | oldFb /= fb
+ && s `between` (a,b) = step s (f s) orig
+ | otherwise = step m (f m) orig
+ where
+ s = b - (b * oldB) * realToFrac (fb / (fb - oldFb))
+ m = 0.5 * (a + b)
+
+ estimateRoot (Dekker _ _ b _ _ _) = b
+ estimateError (Dekker a _ b _ _ _) = a - b
+
+between :: Ord a => a -> (a,a) -> Bool
+a `between` (x,y) = (a > min x y) && (a < max x y)
+
+-- |Incorporates a new point, maintaining invariant 1, assuming invariant 3,
+-- and using 'accept' to restore invariant 2.
+step :: (Num b, Ord b) => a -> b -> Dekker a b -> Dekker a b
+step x fx orig@(Dekker a fa b fb _ _)
+ | signum fx /= signum fa = accept a fa x fx orig
+ | otherwise = accept x fx b fb orig
+
+-- |Re-establishes invariant 2 (abs fb <= abs fa) without affecting invariants 1 and 3.
+accept :: (Num b, Ord b) => a -> b -> a -> b -> Dekker a b -> Dekker a b
+accept a fa b fb (Dekker _ _ oldB oldFb _ _)
+ | abs fb <= abs fa = Dekker a fa b fb oldB oldFb
+ | otherwise = Dekker b fb a fa oldB oldFb \ No newline at end of file
diff --git a/src/Math/Root/Finder/FalsePosition.hs b/src/Math/Root/Finder/FalsePosition.hs
new file mode 100644
index 0000000..cf4a4be
--- /dev/null
+++ b/src/Math/Root/Finder/FalsePosition.hs
@@ -0,0 +1,46 @@
+{-# LANGUAGE MultiParamTypeClasses, FlexibleInstances #-}
+module Math.Root.Finder.FalsePosition
+ ( FalsePosition, falsePosition
+ ) where
+
+import Math.Root.Finder
+
+-- | @falsePosition f x1 x2 xacc@: Using the false-position method, return a
+-- root of a function known to lie between x1 and x2. The root is refined
+-- until its accuracy is += xacc.
+falsePosition :: (Ord a, Fractional a) => (a -> a) -> a -> a -> a -> Either (FalsePosition a a) a
+falsePosition f x1 x2 xacc = fmap estimateRoot (findRoot f x1 x2 xacc)
+
+-- |Iteratively refine a bracketing interval [x1, x2] of a root of f
+-- until total convergence (which may or may not ever be achieved) using
+-- the false-position method.
+data FalsePosition a b = FalsePosition
+ { fpRoot :: !a
+ , fpDX :: !a
+ , _fpXL :: !a
+ , _fpFL :: !a
+ , _fpXH :: !a
+ , _fpFH :: !a
+ } deriving (Eq, Show)
+
+instance (Fractional a, Ord a) => RootFinder FalsePosition a a where
+ initRootFinder f x1 x2
+ -- step once to compute first estimate
+ | f1 <= 0 && f2 >= 0
+ || f2 <= 0 && f1 >= 0 = stepRootFinder f $ FalsePosition 0 0 x2 f2 x1 f1
+ | otherwise = error "FalsePosition: given interval does not bracket a root"
+ where
+ f1 = f x1
+ f2 = f x2
+
+ stepRootFinder f (FalsePosition _ _ xl fl xh fh) = case compare fNew 0 of
+ LT -> FalsePosition xNew (xl - xNew) xNew fNew xh fh
+ EQ -> FalsePosition xNew 0 xNew fNew xNew fNew
+ GT -> FalsePosition xNew (xh - xNew) xl fl xNew fNew
+ where
+ dx = xh - xl
+ xNew = xl + dx * fl / (fl - fh)
+ fNew = f xNew
+
+ estimateRoot = fpRoot
+ estimateError = fpDX
diff --git a/src/Math/Root/Finder/InverseQuadratic.hs b/src/Math/Root/Finder/InverseQuadratic.hs
new file mode 100644
index 0000000..abfe7d0
--- /dev/null
+++ b/src/Math/Root/Finder/InverseQuadratic.hs
@@ -0,0 +1,40 @@
+{-# LANGUAGE
+ MultiParamTypeClasses,
+ FlexibleInstances
+ #-}
+module Math.Root.Finder.InverseQuadratic (InverseQuadratic, inverseQuadratic) where
+
+import Math.Root.Finder
+
+data InverseQuadratic a b = InverseQuadratic !a !b !a !b !a !b
+ deriving (Eq, Show)
+
+-- |@inverseQuadratic f x1 x2 xacc@: attempt to find a root of a function
+-- known to lie between x1 and x2, using the inverse quadratic interpolation
+-- method. The root will be refined till its accuracy is +-xacc. If
+-- convergence fails, returns the final state of the search.
+inverseQuadratic :: RealFloat a => (a -> a) -> a -> a -> a -> Either (InverseQuadratic a a) a
+inverseQuadratic f x1 x2 xacc = fmap estimateRoot (findRoot f x1 x2 xacc)
+
+instance (Fractional a, Ord a, Real b, Fractional b) => RootFinder InverseQuadratic a b where
+ initRootFinder f x1 x2 = InverseQuadratic x0 (f x0) x1 (f x1) x2 (f x2)
+ where x0 = 0.5 * (x1 + x2)
+ stepRootFinder f orig@(InverseQuadratic x0 f0 x1 f1 x2 f2)
+ | f1 /= f2 = InverseQuadratic newX newF x0 f0 x1 f1
+ | otherwise = orig
+ where
+ newX
+ | f0 /= f1 && f0 /= f2
+ = let a = realToFrac (f0 / (f2 - f1))
+ b = realToFrac (f1 / (f2 - f0))
+ c = realToFrac (f2 / (f1 - f0))
+ in (a * b * x2) - (a * c * x1) + (b * c * x0)
+ | otherwise
+ -- Fall back to secant method (linear interpolation)
+ -- when quadratic interpolation will yield nonsensical results.
+ = x1 - realToFrac f1 * (x1 - x2) / realToFrac (f1 - f2)
+ newF = f newX
+
+ estimateRoot (InverseQuadratic x0 _ _ _ _ _) = x0
+ estimateError (InverseQuadratic x0 _ x1 _ x2 _) =
+ maximum [x0, x1, x2] - minimum [x0, x1, x2]
diff --git a/src/Math/Root/Finder/Newton.hs b/src/Math/Root/Finder/Newton.hs
new file mode 100644
index 0000000..5ac836e
--- /dev/null
+++ b/src/Math/Root/Finder/Newton.hs
@@ -0,0 +1,33 @@
+{-# LANGUAGE MultiParamTypeClasses, FlexibleInstances #-}
+module Math.Root.Finder.Newton
+ ( Newton, newton
+ ) where
+
+import Math.Root.Finder
+
+data Newton a b = Newton
+ { newtRTN :: !a
+ , newtDX :: a
+ } deriving (Eq, Show)
+
+instance Fractional a => RootFinder Newton a (a,a) where
+ initRootFinder f x1 x2 = stepRootFinder f (Newton rtn undefined)
+ where
+ rtn = 0.5 * (x1 + x2)
+
+ stepRootFinder f Newton{newtRTN = rtn} = Newton (rtn - dx) dx
+ where
+ (y,dy) = f rtn
+ dx = y / dy
+
+ estimateRoot Newton{newtRTN = rtn} = rtn
+ estimateError Newton{newtDX = dx} = dx
+
+-- | @newton f x1 x2 xacc@: using Newton's method, return a root of a
+-- function known to lie between x1 and x2. The root is refined until its
+-- accuracy is += xacc.
+--
+-- The function passed should return a pair containing the value of the
+-- function and its derivative, respectively.
+newton :: (Ord a, Fractional a) => (a -> (a, a)) -> a -> a -> a -> Either (Newton a (a,a)) a
+newton f x1 x2 xacc = fmap estimateRoot (findRoot f x1 x2 xacc)
diff --git a/src/Math/Root/Finder/Ridders.hs b/src/Math/Root/Finder/Ridders.hs
new file mode 100644
index 0000000..03cafe9
--- /dev/null
+++ b/src/Math/Root/Finder/Ridders.hs
@@ -0,0 +1,58 @@
+{-# LANGUAGE MultiParamTypeClasses, FlexibleInstances #-}
+{-# OPTIONS_GHC -fno-warn-name-shadowing #-}
+module Math.Root.Finder.Ridders
+ ( RiddersMethod, ridders
+ ) where
+
+import Math.Root.Finder
+
+-- |@ridders f x1 x2 xacc@: attempt to find a root of a function known to
+-- lie between x1 and x2, using Ridders' method. The root will be refined
+-- till its accuracy is +-xacc. If convergence fails, returns the final
+-- state of the search.
+ridders :: (Ord a, Floating a) => (a -> a) -> a -> a -> a -> Either (RiddersMethod a a) a
+ridders f x1 x2 xacc = fmap estimateRoot (findRoot f x1 x2 xacc)
+
+data RiddersMethod a b
+ = ConvergedRidders !a
+ | RiddersMethod
+ { ridXL :: !a
+ , _ridFL :: !b
+ , ridXH :: !a
+ , _ridFH :: !b
+ } deriving (Eq, Show)
+
+instance (Floating a, Ord a) => RootFinder RiddersMethod a a where
+ initRootFinder f x1 x2
+ | f1 < 0 && f2 < 0
+ || f2 > 0 && f1 > 0 = error "riddersMethod: interval does not bracket a root"
+ | otherwise = RiddersMethod x1 f1 x2 f2
+ where
+ f1 = f x1
+ f2 = f x2
+ stepRootFinder _ orig@ConvergedRidders{} = orig
+ stepRootFinder f (RiddersMethod xl fl xh fh)
+ | signNEQ fm fNew = finish xNew fNew xm fm
+ | signNEQ fl fNew = finish xNew fNew xl fl
+ | signNEQ fh fNew = finish xNew fNew xh fh
+ | otherwise = error "RiddersMethod: encountered singularity"
+ where
+ xm = 0.5 * (xl + xh)
+ fm = f xm
+ s = sqrt (fm*fm - fl*fh)
+ xNew = xm + (xm-xl)*((if fl >= fh then id else negate) fm / s)
+ fNew = f xNew
+
+ signNEQ a b = a /= 0 && signum b /= signum a
+
+ finish xl fl xh fh
+ | xl == xh = ConvergedRidders xl
+ | fl == 0 = ConvergedRidders xl
+ | fh == 0 = ConvergedRidders xh
+ | otherwise = RiddersMethod xl fl xh fh
+
+ estimateRoot (ConvergedRidders x) = x
+ estimateRoot RiddersMethod{ridXL = x} = x
+
+ estimateError ConvergedRidders{} = 0
+ estimateError RiddersMethod{ridXL = xl, ridXH = xh} = xl - xh
diff --git a/src/Math/Root/Finder/Secant.hs b/src/Math/Root/Finder/Secant.hs
new file mode 100644
index 0000000..fc6072f
--- /dev/null
+++ b/src/Math/Root/Finder/Secant.hs
@@ -0,0 +1,46 @@
+{-# LANGUAGE MultiParamTypeClasses, FlexibleInstances #-}
+module Math.Root.Finder.Secant
+ ( SecantMethod, secant
+ ) where
+
+import Math.Root.Finder
+
+-- | @secant f x1 x2 xacc@: Using the secant method, return the root of a
+-- function thought to lie between x1 and x2. The root is refined until its
+-- accuracy is +-xacc.
+secant :: (Ord a, Fractional a) => (a -> a) -> a -> a -> a -> Either (SecantMethod a a) a
+secant f x1 x2 xacc = fmap estimateRoot (findRoot f x1 x2 xacc)
+
+-- |Iteratively refine 2 estimates x1, x2 of a root of f until total
+-- convergence (which may or may not ever be achieved) using the
+-- secant method.
+data SecantMethod a b
+ = ConvergedSecantMethod !a
+ | SecantMethod
+ { secDX :: !a
+ , secXL :: !a
+ , _secFL :: !b
+ , _secRTS :: !a
+ , _secFRTS :: !b
+ } deriving (Eq, Show)
+
+instance (Fractional a, Ord a) => RootFinder SecantMethod a a where
+ initRootFinder f x1 x2
+ | abs f1 < abs f2 = stepRootFinder f $ SecantMethod 0 x2 f2 x1 f1
+ | otherwise = stepRootFinder f $ SecantMethod 0 x1 f1 x2 f2
+ where f1 = f x1; f2 = f x2
+
+ stepRootFinder _ orig@ConvergedSecantMethod{} = orig
+ stepRootFinder f (SecantMethod _ xl fl rts fRts)
+ | fNew == 0 = ConvergedSecantMethod xNew
+ | otherwise = SecantMethod dx rts fRts xNew fNew
+ where
+ dx = (xl - rts) * fRts / (fRts - fl)
+ xNew = rts + dx
+ fNew = f xNew
+
+ estimateRoot (ConvergedSecantMethod x) = x
+ estimateRoot SecantMethod{secXL = x} = x
+
+ estimateError ConvergedSecantMethod{} = 0
+ estimateError SecantMethod{secDX = dx} = dx