-----------------------------------------------------------------------------
--
-- Module      :  Math.Roots.Bisection
-- Copyright   :  (c) 2013-16 Brian W Bush
-- License     :  MIT
--
-- Maintainer  :  Brian W Bush <consult@brianwbush.info>
-- Stability   :  Stable
-- Portability :  Portable
--
-- | Root finding using the bisection method.
--
-----------------------------------------------------------------------------


{-# LANGUAGE Safe #-}


module Math.Roots.Bisection (
-- * Functions.
  findRoot
) where


import Math.Roots (RootFinder)


-- | Find a root using the bisection method.
-- |
-- | Reference: William H. Press, Saul A. Teukolsky, William T. Vetterling, and Brian P. Flannery, Numerical Recipes in C: The Art of Scientific Computing>, Second Edition (New York: Cambridge Univ. Press, 1992).
findRoot :: (Fractional a, Num a, Ord a) =>
     Int          -- ^ Maximum number of bisections.
  -> a            -- ^ Absolute error tolerance.
  -> RootFinder a -- ^ The root finder.
findRoot jmax xacc func (x1, x2)
  | x1 >= x2               = Left "The root is not bracketed by an interval."
  | func x1 * func x2 >= 0 = Left "The root is not bracketed by the interval."
  | otherwise              =
    findRoot' jmax $
      if func x1 < 0
        then (x2 - x1, x1)
        else (x1 - x2, x2)
      where
        findRoot' 0 _ = Left "The maximum number of bisections was reached."
        findRoot' j (dx, rtb) =
          do
            let
              dx' = dx / 2
              xmid = rtb + dx'
              fmid = func xmid
              rtb' =
                if fmid <= 0
                  then xmid
                  else rtb
            if abs dx < xacc || fmid == 0
              then return rtb'
              else findRoot' (j - 1) (dx', rtb')