diff --git a/Setup.lhs b/Setup.lhs
new file mode 100644
--- /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
--- /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
--- /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
--- /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
--- /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
--- /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
diff --git a/src/Math/Root/Finder/Dekker.hs b/src/Math/Root/Finder/Dekker.hs
new file mode 100644
--- /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
diff --git a/src/Math/Root/Finder/FalsePosition.hs b/src/Math/Root/Finder/FalsePosition.hs
new file mode 100644
--- /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
--- /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
--- /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
--- /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
--- /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