{-# LANGUAGE BangPatterns       #-}
{-# LANGUAGE DeriveDataTypeable #-}
{-# LANGUAGE DeriveFoldable     #-}
{-# LANGUAGE DeriveGeneric      #-}
{-# LANGUAGE DeriveTraversable  #-}
{-# LANGUAGE TypeFamilies       #-}
-- |
-- Module    : Numeric.RootFinding
-- Copyright : (c) 2011 Bryan O'Sullivan, 2018 Alexey Khudyakov
-- License   : BSD3
--
-- Maintainer  : bos@serpentine.com
-- Stability   : experimental
-- Portability : portable
--
-- Haskell functions for finding the roots of real functions of real
-- arguments. These algorithms are iterative so we provide both
-- function returning root (or failure to find root) and list of
-- iterations.
module Numeric.RootFinding
    ( -- * Data types
      Root(..)
    , fromRoot
    , Tolerance(..)
    , withinTolerance
    , IterationStep(..)
    , findRoot
    -- * Ridders algorithm
    , RiddersParam(..)
    , ridders
    , riddersIterations
    , RiddersStep(..)
    -- * Newton-Raphson algorithm
    , NewtonParam(..)
    , newtonRaphson
    , newtonRaphsonIterations
    , NewtonStep(..)
    -- * References
    -- $references
    ) where

import Control.Applicative              (Alternative(..))
import Control.Monad                    (MonadPlus(..), ap)
import Control.DeepSeq                  (NFData(..))
import Data.Data                        (Data, Typeable)
import Data.Default.Class
import GHC.Generics                     (Generic)
import Numeric.MathFunctions.Comparison (within,eqRelErr)
import Numeric.MathFunctions.Constants  (m_epsilon)



----------------------------------------------------------------
-- Data types
----------------------------------------------------------------

-- | The result of searching for a root of a mathematical function.
data Root a
  = NotBracketed
    -- ^ The function does not have opposite signs when
    -- evaluated at the lower and upper bounds of the search.
  | SearchFailed
    -- ^ The search failed to converge to within the given
    -- error tolerance after the given number of iterations.
  | Root !a
    -- ^ A root was successfully found.
  deriving (Root a -> Root a -> Bool
(Root a -> Root a -> Bool)
-> (Root a -> Root a -> Bool) -> Eq (Root a)
forall a. Eq a => Root a -> Root a -> Bool
forall a. (a -> a -> Bool) -> (a -> a -> Bool) -> Eq a
$c== :: forall a. Eq a => Root a -> Root a -> Bool
== :: Root a -> Root a -> Bool
$c/= :: forall a. Eq a => Root a -> Root a -> Bool
/= :: Root a -> Root a -> Bool
Eq, ReadPrec [Root a]
ReadPrec (Root a)
Int -> ReadS (Root a)
ReadS [Root a]
(Int -> ReadS (Root a))
-> ReadS [Root a]
-> ReadPrec (Root a)
-> ReadPrec [Root a]
-> Read (Root a)
forall a. Read a => ReadPrec [Root a]
forall a. Read a => ReadPrec (Root a)
forall a. Read a => Int -> ReadS (Root a)
forall a. Read a => ReadS [Root a]
forall a.
(Int -> ReadS a)
-> ReadS [a] -> ReadPrec a -> ReadPrec [a] -> Read a
$creadsPrec :: forall a. Read a => Int -> ReadS (Root a)
readsPrec :: Int -> ReadS (Root a)
$creadList :: forall a. Read a => ReadS [Root a]
readList :: ReadS [Root a]
$creadPrec :: forall a. Read a => ReadPrec (Root a)
readPrec :: ReadPrec (Root a)
$creadListPrec :: forall a. Read a => ReadPrec [Root a]
readListPrec :: ReadPrec [Root a]
Read, Int -> Root a -> ShowS
[Root a] -> ShowS
Root a -> String
(Int -> Root a -> ShowS)
-> (Root a -> String) -> ([Root a] -> ShowS) -> Show (Root a)
forall a. Show a => Int -> Root a -> ShowS
forall a. Show a => [Root a] -> ShowS
forall a. Show a => Root a -> String
forall a.
(Int -> a -> ShowS) -> (a -> String) -> ([a] -> ShowS) -> Show a
$cshowsPrec :: forall a. Show a => Int -> Root a -> ShowS
showsPrec :: Int -> Root a -> ShowS
$cshow :: forall a. Show a => Root a -> String
show :: Root a -> String
$cshowList :: forall a. Show a => [Root a] -> ShowS
showList :: [Root a] -> ShowS
Show, Typeable, Typeable (Root a)
Typeable (Root a) =>
(forall (c :: * -> *).
 (forall d b. Data d => c (d -> b) -> d -> c b)
 -> (forall g. g -> c g) -> Root a -> c (Root a))
-> (forall (c :: * -> *).
    (forall b r. Data b => c (b -> r) -> c r)
    -> (forall r. r -> c r) -> Constr -> c (Root a))
-> (Root a -> Constr)
-> (Root a -> DataType)
-> (forall (t :: * -> *) (c :: * -> *).
    Typeable t =>
    (forall d. Data d => c (t d)) -> Maybe (c (Root a)))
-> (forall (t :: * -> * -> *) (c :: * -> *).
    Typeable t =>
    (forall d e. (Data d, Data e) => c (t d e)) -> Maybe (c (Root a)))
-> ((forall b. Data b => b -> b) -> Root a -> Root a)
-> (forall r r'.
    (r -> r' -> r)
    -> r -> (forall d. Data d => d -> r') -> Root a -> r)
-> (forall r r'.
    (r' -> r -> r)
    -> r -> (forall d. Data d => d -> r') -> Root a -> r)
-> (forall u. (forall d. Data d => d -> u) -> Root a -> [u])
-> (forall u. Int -> (forall d. Data d => d -> u) -> Root a -> u)
-> (forall (m :: * -> *).
    Monad m =>
    (forall d. Data d => d -> m d) -> Root a -> m (Root a))
-> (forall (m :: * -> *).
    MonadPlus m =>
    (forall d. Data d => d -> m d) -> Root a -> m (Root a))
-> (forall (m :: * -> *).
    MonadPlus m =>
    (forall d. Data d => d -> m d) -> Root a -> m (Root a))
-> Data (Root a)
Root a -> Constr
Root a -> DataType
(forall b. Data b => b -> b) -> Root a -> Root a
forall a. Data a => Typeable (Root a)
forall a. Data a => Root a -> Constr
forall a. Data a => Root a -> DataType
forall a.
Data a =>
(forall b. Data b => b -> b) -> Root a -> Root a
forall a u.
Data a =>
Int -> (forall d. Data d => d -> u) -> Root a -> u
forall a u. Data a => (forall d. Data d => d -> u) -> Root a -> [u]
forall a r r'.
Data a =>
(r -> r' -> r) -> r -> (forall d. Data d => d -> r') -> Root a -> r
forall a r r'.
Data a =>
(r' -> r -> r) -> r -> (forall d. Data d => d -> r') -> Root a -> r
forall a (m :: * -> *).
(Data a, Monad m) =>
(forall d. Data d => d -> m d) -> Root a -> m (Root a)
forall a (m :: * -> *).
(Data a, MonadPlus m) =>
(forall d. Data d => d -> m d) -> Root a -> m (Root a)
forall a (c :: * -> *).
Data a =>
(forall b r. Data b => c (b -> r) -> c r)
-> (forall r. r -> c r) -> Constr -> c (Root a)
forall a (c :: * -> *).
Data a =>
(forall d b. Data d => c (d -> b) -> d -> c b)
-> (forall g. g -> c g) -> Root a -> c (Root a)
forall a (t :: * -> *) (c :: * -> *).
(Data a, Typeable t) =>
(forall d. Data d => c (t d)) -> Maybe (c (Root a))
forall a (t :: * -> * -> *) (c :: * -> *).
(Data a, Typeable t) =>
(forall d e. (Data d, Data e) => c (t d e)) -> Maybe (c (Root a))
forall a.
Typeable a =>
(forall (c :: * -> *).
 (forall d b. Data d => c (d -> b) -> d -> c b)
 -> (forall g. g -> c g) -> a -> c a)
-> (forall (c :: * -> *).
    (forall b r. Data b => c (b -> r) -> c r)
    -> (forall r. r -> c r) -> Constr -> c a)
-> (a -> Constr)
-> (a -> DataType)
-> (forall (t :: * -> *) (c :: * -> *).
    Typeable t =>
    (forall d. Data d => c (t d)) -> Maybe (c a))
-> (forall (t :: * -> * -> *) (c :: * -> *).
    Typeable t =>
    (forall d e. (Data d, Data e) => c (t d e)) -> Maybe (c a))
-> ((forall b. Data b => b -> b) -> a -> a)
-> (forall r r'.
    (r -> r' -> r) -> r -> (forall d. Data d => d -> r') -> a -> r)
-> (forall r r'.
    (r' -> r -> r) -> r -> (forall d. Data d => d -> r') -> a -> r)
-> (forall u. (forall d. Data d => d -> u) -> a -> [u])
-> (forall u. Int -> (forall d. Data d => d -> u) -> a -> u)
-> (forall (m :: * -> *).
    Monad m =>
    (forall d. Data d => d -> m d) -> a -> m a)
-> (forall (m :: * -> *).
    MonadPlus m =>
    (forall d. Data d => d -> m d) -> a -> m a)
-> (forall (m :: * -> *).
    MonadPlus m =>
    (forall d. Data d => d -> m d) -> a -> m a)
-> Data a
forall u. Int -> (forall d. Data d => d -> u) -> Root a -> u
forall u. (forall d. Data d => d -> u) -> Root a -> [u]
forall r r'.
(r -> r' -> r) -> r -> (forall d. Data d => d -> r') -> Root a -> r
forall r r'.
(r' -> r -> r) -> r -> (forall d. Data d => d -> r') -> Root a -> r
forall (m :: * -> *).
Monad m =>
(forall d. Data d => d -> m d) -> Root a -> m (Root a)
forall (m :: * -> *).
MonadPlus m =>
(forall d. Data d => d -> m d) -> Root a -> m (Root a)
forall (c :: * -> *).
(forall b r. Data b => c (b -> r) -> c r)
-> (forall r. r -> c r) -> Constr -> c (Root a)
forall (c :: * -> *).
(forall d b. Data d => c (d -> b) -> d -> c b)
-> (forall g. g -> c g) -> Root a -> c (Root a)
forall (t :: * -> *) (c :: * -> *).
Typeable t =>
(forall d. Data d => c (t d)) -> Maybe (c (Root a))
forall (t :: * -> * -> *) (c :: * -> *).
Typeable t =>
(forall d e. (Data d, Data e) => c (t d e)) -> Maybe (c (Root a))
$cgfoldl :: forall a (c :: * -> *).
Data a =>
(forall d b. Data d => c (d -> b) -> d -> c b)
-> (forall g. g -> c g) -> Root a -> c (Root a)
gfoldl :: forall (c :: * -> *).
(forall d b. Data d => c (d -> b) -> d -> c b)
-> (forall g. g -> c g) -> Root a -> c (Root a)
$cgunfold :: forall a (c :: * -> *).
Data a =>
(forall b r. Data b => c (b -> r) -> c r)
-> (forall r. r -> c r) -> Constr -> c (Root a)
gunfold :: forall (c :: * -> *).
(forall b r. Data b => c (b -> r) -> c r)
-> (forall r. r -> c r) -> Constr -> c (Root a)
$ctoConstr :: forall a. Data a => Root a -> Constr
toConstr :: Root a -> Constr
$cdataTypeOf :: forall a. Data a => Root a -> DataType
dataTypeOf :: Root a -> DataType
$cdataCast1 :: forall a (t :: * -> *) (c :: * -> *).
(Data a, Typeable t) =>
(forall d. Data d => c (t d)) -> Maybe (c (Root a))
dataCast1 :: forall (t :: * -> *) (c :: * -> *).
Typeable t =>
(forall d. Data d => c (t d)) -> Maybe (c (Root a))
$cdataCast2 :: forall a (t :: * -> * -> *) (c :: * -> *).
(Data a, Typeable t) =>
(forall d e. (Data d, Data e) => c (t d e)) -> Maybe (c (Root a))
dataCast2 :: forall (t :: * -> * -> *) (c :: * -> *).
Typeable t =>
(forall d e. (Data d, Data e) => c (t d e)) -> Maybe (c (Root a))
$cgmapT :: forall a.
Data a =>
(forall b. Data b => b -> b) -> Root a -> Root a
gmapT :: (forall b. Data b => b -> b) -> Root a -> Root a
$cgmapQl :: forall a r r'.
Data a =>
(r -> r' -> r) -> r -> (forall d. Data d => d -> r') -> Root a -> r
gmapQl :: forall r r'.
(r -> r' -> r) -> r -> (forall d. Data d => d -> r') -> Root a -> r
$cgmapQr :: forall a r r'.
Data a =>
(r' -> r -> r) -> r -> (forall d. Data d => d -> r') -> Root a -> r
gmapQr :: forall r r'.
(r' -> r -> r) -> r -> (forall d. Data d => d -> r') -> Root a -> r
$cgmapQ :: forall a u. Data a => (forall d. Data d => d -> u) -> Root a -> [u]
gmapQ :: forall u. (forall d. Data d => d -> u) -> Root a -> [u]
$cgmapQi :: forall a u.
Data a =>
Int -> (forall d. Data d => d -> u) -> Root a -> u
gmapQi :: forall u. Int -> (forall d. Data d => d -> u) -> Root a -> u
$cgmapM :: forall a (m :: * -> *).
(Data a, Monad m) =>
(forall d. Data d => d -> m d) -> Root a -> m (Root a)
gmapM :: forall (m :: * -> *).
Monad m =>
(forall d. Data d => d -> m d) -> Root a -> m (Root a)
$cgmapMp :: forall a (m :: * -> *).
(Data a, MonadPlus m) =>
(forall d. Data d => d -> m d) -> Root a -> m (Root a)
gmapMp :: forall (m :: * -> *).
MonadPlus m =>
(forall d. Data d => d -> m d) -> Root a -> m (Root a)
$cgmapMo :: forall a (m :: * -> *).
(Data a, MonadPlus m) =>
(forall d. Data d => d -> m d) -> Root a -> m (Root a)
gmapMo :: forall (m :: * -> *).
MonadPlus m =>
(forall d. Data d => d -> m d) -> Root a -> m (Root a)
Data, (forall m. Monoid m => Root m -> m)
-> (forall m a. Monoid m => (a -> m) -> Root a -> m)
-> (forall m a. Monoid m => (a -> m) -> Root a -> m)
-> (forall a b. (a -> b -> b) -> b -> Root a -> b)
-> (forall a b. (a -> b -> b) -> b -> Root a -> b)
-> (forall b a. (b -> a -> b) -> b -> Root a -> b)
-> (forall b a. (b -> a -> b) -> b -> Root a -> b)
-> (forall a. (a -> a -> a) -> Root a -> a)
-> (forall a. (a -> a -> a) -> Root a -> a)
-> (forall a. Root a -> [a])
-> (forall a. Root a -> Bool)
-> (forall a. Root a -> Int)
-> (forall a. Eq a => a -> Root a -> Bool)
-> (forall a. Ord a => Root a -> a)
-> (forall a. Ord a => Root a -> a)
-> (forall a. Num a => Root a -> a)
-> (forall a. Num a => Root a -> a)
-> Foldable Root
forall a. Eq a => a -> Root a -> Bool
forall a. Num a => Root a -> a
forall a. Ord a => Root a -> a
forall m. Monoid m => Root m -> m
forall a. Root a -> Bool
forall a. Root a -> Int
forall a. Root a -> [a]
forall a. (a -> a -> a) -> Root a -> a
forall m a. Monoid m => (a -> m) -> Root a -> m
forall b a. (b -> a -> b) -> b -> Root a -> b
forall a b. (a -> b -> b) -> b -> Root a -> b
forall (t :: * -> *).
(forall m. Monoid m => t m -> m)
-> (forall m a. Monoid m => (a -> m) -> t a -> m)
-> (forall m a. Monoid m => (a -> m) -> t a -> m)
-> (forall a b. (a -> b -> b) -> b -> t a -> b)
-> (forall a b. (a -> b -> b) -> b -> t a -> b)
-> (forall b a. (b -> a -> b) -> b -> t a -> b)
-> (forall b a. (b -> a -> b) -> b -> t a -> b)
-> (forall a. (a -> a -> a) -> t a -> a)
-> (forall a. (a -> a -> a) -> t a -> a)
-> (forall a. t a -> [a])
-> (forall a. t a -> Bool)
-> (forall a. t a -> Int)
-> (forall a. Eq a => a -> t a -> Bool)
-> (forall a. Ord a => t a -> a)
-> (forall a. Ord a => t a -> a)
-> (forall a. Num a => t a -> a)
-> (forall a. Num a => t a -> a)
-> Foldable t
$cfold :: forall m. Monoid m => Root m -> m
fold :: forall m. Monoid m => Root m -> m
$cfoldMap :: forall m a. Monoid m => (a -> m) -> Root a -> m
foldMap :: forall m a. Monoid m => (a -> m) -> Root a -> m
$cfoldMap' :: forall m a. Monoid m => (a -> m) -> Root a -> m
foldMap' :: forall m a. Monoid m => (a -> m) -> Root a -> m
$cfoldr :: forall a b. (a -> b -> b) -> b -> Root a -> b
foldr :: forall a b. (a -> b -> b) -> b -> Root a -> b
$cfoldr' :: forall a b. (a -> b -> b) -> b -> Root a -> b
foldr' :: forall a b. (a -> b -> b) -> b -> Root a -> b
$cfoldl :: forall b a. (b -> a -> b) -> b -> Root a -> b
foldl :: forall b a. (b -> a -> b) -> b -> Root a -> b
$cfoldl' :: forall b a. (b -> a -> b) -> b -> Root a -> b
foldl' :: forall b a. (b -> a -> b) -> b -> Root a -> b
$cfoldr1 :: forall a. (a -> a -> a) -> Root a -> a
foldr1 :: forall a. (a -> a -> a) -> Root a -> a
$cfoldl1 :: forall a. (a -> a -> a) -> Root a -> a
foldl1 :: forall a. (a -> a -> a) -> Root a -> a
$ctoList :: forall a. Root a -> [a]
toList :: forall a. Root a -> [a]
$cnull :: forall a. Root a -> Bool
null :: forall a. Root a -> Bool
$clength :: forall a. Root a -> Int
length :: forall a. Root a -> Int
$celem :: forall a. Eq a => a -> Root a -> Bool
elem :: forall a. Eq a => a -> Root a -> Bool
$cmaximum :: forall a. Ord a => Root a -> a
maximum :: forall a. Ord a => Root a -> a
$cminimum :: forall a. Ord a => Root a -> a
minimum :: forall a. Ord a => Root a -> a
$csum :: forall a. Num a => Root a -> a
sum :: forall a. Num a => Root a -> a
$cproduct :: forall a. Num a => Root a -> a
product :: forall a. Num a => Root a -> a
Foldable, Functor Root
Foldable Root
(Functor Root, Foldable Root) =>
(forall (f :: * -> *) a b.
 Applicative f =>
 (a -> f b) -> Root a -> f (Root b))
-> (forall (f :: * -> *) a.
    Applicative f =>
    Root (f a) -> f (Root a))
-> (forall (m :: * -> *) a b.
    Monad m =>
    (a -> m b) -> Root a -> m (Root b))
-> (forall (m :: * -> *) a. Monad m => Root (m a) -> m (Root a))
-> Traversable Root
forall (t :: * -> *).
(Functor t, Foldable t) =>
(forall (f :: * -> *) a b.
 Applicative f =>
 (a -> f b) -> t a -> f (t b))
-> (forall (f :: * -> *) a. Applicative f => t (f a) -> f (t a))
-> (forall (m :: * -> *) a b.
    Monad m =>
    (a -> m b) -> t a -> m (t b))
-> (forall (m :: * -> *) a. Monad m => t (m a) -> m (t a))
-> Traversable t
forall (m :: * -> *) a. Monad m => Root (m a) -> m (Root a)
forall (f :: * -> *) a. Applicative f => Root (f a) -> f (Root a)
forall (m :: * -> *) a b.
Monad m =>
(a -> m b) -> Root a -> m (Root b)
forall (f :: * -> *) a b.
Applicative f =>
(a -> f b) -> Root a -> f (Root b)
$ctraverse :: forall (f :: * -> *) a b.
Applicative f =>
(a -> f b) -> Root a -> f (Root b)
traverse :: forall (f :: * -> *) a b.
Applicative f =>
(a -> f b) -> Root a -> f (Root b)
$csequenceA :: forall (f :: * -> *) a. Applicative f => Root (f a) -> f (Root a)
sequenceA :: forall (f :: * -> *) a. Applicative f => Root (f a) -> f (Root a)
$cmapM :: forall (m :: * -> *) a b.
Monad m =>
(a -> m b) -> Root a -> m (Root b)
mapM :: forall (m :: * -> *) a b.
Monad m =>
(a -> m b) -> Root a -> m (Root b)
$csequence :: forall (m :: * -> *) a. Monad m => Root (m a) -> m (Root a)
sequence :: forall (m :: * -> *) a. Monad m => Root (m a) -> m (Root a)
Traversable, (forall a b. (a -> b) -> Root a -> Root b)
-> (forall a b. a -> Root b -> Root a) -> Functor Root
forall a b. a -> Root b -> Root a
forall a b. (a -> b) -> Root a -> Root b
forall (f :: * -> *).
(forall a b. (a -> b) -> f a -> f b)
-> (forall a b. a -> f b -> f a) -> Functor f
$cfmap :: forall a b. (a -> b) -> Root a -> Root b
fmap :: forall a b. (a -> b) -> Root a -> Root b
$c<$ :: forall a b. a -> Root b -> Root a
<$ :: forall a b. a -> Root b -> Root a
Functor, (forall x. Root a -> Rep (Root a) x)
-> (forall x. Rep (Root a) x -> Root a) -> Generic (Root a)
forall x. Rep (Root a) x -> Root a
forall x. Root a -> Rep (Root a) x
forall a.
(forall x. a -> Rep a x) -> (forall x. Rep a x -> a) -> Generic a
forall a x. Rep (Root a) x -> Root a
forall a x. Root a -> Rep (Root a) x
$cfrom :: forall a x. Root a -> Rep (Root a) x
from :: forall x. Root a -> Rep (Root a) x
$cto :: forall a x. Rep (Root a) x -> Root a
to :: forall x. Rep (Root a) x -> Root a
Generic)

instance (NFData a) => NFData (Root a) where
    rnf :: Root a -> ()
rnf Root a
NotBracketed = ()
    rnf Root a
SearchFailed = ()
    rnf (Root a
a)     = a -> ()
forall a. NFData a => a -> ()
rnf a
a

instance Applicative Root where
    pure :: forall a. a -> Root a
pure  = a -> Root a
forall a. a -> Root a
Root
    <*> :: forall a b. Root (a -> b) -> Root a -> Root b
(<*>) = Root (a -> b) -> Root a -> Root b
forall (m :: * -> *) a b. Monad m => m (a -> b) -> m a -> m b
ap

instance Monad Root where
    Root a
NotBracketed >>= :: forall a b. Root a -> (a -> Root b) -> Root b
>>= a -> Root b
_ = Root b
forall a. Root a
NotBracketed
    Root a
SearchFailed >>= a -> Root b
_ = Root b
forall a. Root a
SearchFailed
    Root a
a       >>= a -> Root b
f = a -> Root b
f a
a
    return :: forall a. a -> Root a
return = a -> Root a
forall a. a -> Root a
forall (f :: * -> *) a. Applicative f => a -> f a
pure

instance MonadPlus Root where
    mzero :: forall a. Root a
mzero = Root a
forall a. Root a
forall (f :: * -> *) a. Alternative f => f a
empty
    mplus :: forall a. Root a -> Root a -> Root a
mplus = Root a -> Root a -> Root a
forall a. Root a -> Root a -> Root a
forall (f :: * -> *) a. Alternative f => f a -> f a -> f a
(<|>)

instance Alternative Root where
    empty :: forall a. Root a
empty = Root a
forall a. Root a
NotBracketed
    r :: Root a
r@Root{}     <|> :: forall a. Root a -> Root a -> Root a
<|> Root a
_            = Root a
r
    Root a
_            <|> r :: Root a
r@Root{}     = Root a
r
    Root a
NotBracketed <|> Root a
r            = Root a
r
    Root a
r            <|> Root a
NotBracketed = Root a
r
    Root a
_            <|> Root a
r            = Root a
r

-- | Returns either the result of a search for a root, or the default
-- value if the search failed.
fromRoot :: a                 -- ^ Default value.
         -> Root a            -- ^ Result of search for a root.
         -> a
fromRoot :: forall a. a -> Root a -> a
fromRoot a
_ (Root a
a) = a
a
fromRoot a
a Root a
_        = a
a


-- | Error tolerance for finding root. It describes when root finding
--   algorithm should stop trying to improve approximation.
data Tolerance
  = RelTol !Double
    -- ^ Relative error tolerance. Given @RelTol ε@ two values are
    --   considered approximately equal if
    --   \[ \frac{|a - b|}{|\operatorname{max}(a,b)} < \varepsilon \]
  | AbsTol !Double
    -- ^ Absolute error tolerance. Given @AbsTol δ@ two values are
    --   considered approximately equal if \[ |a - b| < \delta \].
    --   Note that @AbsTol 0@ could be used to require to find
    --   approximation within machine precision.
  deriving (Tolerance -> Tolerance -> Bool
(Tolerance -> Tolerance -> Bool)
-> (Tolerance -> Tolerance -> Bool) -> Eq Tolerance
forall a. (a -> a -> Bool) -> (a -> a -> Bool) -> Eq a
$c== :: Tolerance -> Tolerance -> Bool
== :: Tolerance -> Tolerance -> Bool
$c/= :: Tolerance -> Tolerance -> Bool
/= :: Tolerance -> Tolerance -> Bool
Eq, ReadPrec [Tolerance]
ReadPrec Tolerance
Int -> ReadS Tolerance
ReadS [Tolerance]
(Int -> ReadS Tolerance)
-> ReadS [Tolerance]
-> ReadPrec Tolerance
-> ReadPrec [Tolerance]
-> Read Tolerance
forall a.
(Int -> ReadS a)
-> ReadS [a] -> ReadPrec a -> ReadPrec [a] -> Read a
$creadsPrec :: Int -> ReadS Tolerance
readsPrec :: Int -> ReadS Tolerance
$creadList :: ReadS [Tolerance]
readList :: ReadS [Tolerance]
$creadPrec :: ReadPrec Tolerance
readPrec :: ReadPrec Tolerance
$creadListPrec :: ReadPrec [Tolerance]
readListPrec :: ReadPrec [Tolerance]
Read, Int -> Tolerance -> ShowS
[Tolerance] -> ShowS
Tolerance -> String
(Int -> Tolerance -> ShowS)
-> (Tolerance -> String)
-> ([Tolerance] -> ShowS)
-> Show Tolerance
forall a.
(Int -> a -> ShowS) -> (a -> String) -> ([a] -> ShowS) -> Show a
$cshowsPrec :: Int -> Tolerance -> ShowS
showsPrec :: Int -> Tolerance -> ShowS
$cshow :: Tolerance -> String
show :: Tolerance -> String
$cshowList :: [Tolerance] -> ShowS
showList :: [Tolerance] -> ShowS
Show, Typeable, Typeable Tolerance
Typeable Tolerance =>
(forall (c :: * -> *).
 (forall d b. Data d => c (d -> b) -> d -> c b)
 -> (forall g. g -> c g) -> Tolerance -> c Tolerance)
-> (forall (c :: * -> *).
    (forall b r. Data b => c (b -> r) -> c r)
    -> (forall r. r -> c r) -> Constr -> c Tolerance)
-> (Tolerance -> Constr)
-> (Tolerance -> DataType)
-> (forall (t :: * -> *) (c :: * -> *).
    Typeable t =>
    (forall d. Data d => c (t d)) -> Maybe (c Tolerance))
-> (forall (t :: * -> * -> *) (c :: * -> *).
    Typeable t =>
    (forall d e. (Data d, Data e) => c (t d e)) -> Maybe (c Tolerance))
-> ((forall b. Data b => b -> b) -> Tolerance -> Tolerance)
-> (forall r r'.
    (r -> r' -> r)
    -> r -> (forall d. Data d => d -> r') -> Tolerance -> r)
-> (forall r r'.
    (r' -> r -> r)
    -> r -> (forall d. Data d => d -> r') -> Tolerance -> r)
-> (forall u. (forall d. Data d => d -> u) -> Tolerance -> [u])
-> (forall u.
    Int -> (forall d. Data d => d -> u) -> Tolerance -> u)
-> (forall (m :: * -> *).
    Monad m =>
    (forall d. Data d => d -> m d) -> Tolerance -> m Tolerance)
-> (forall (m :: * -> *).
    MonadPlus m =>
    (forall d. Data d => d -> m d) -> Tolerance -> m Tolerance)
-> (forall (m :: * -> *).
    MonadPlus m =>
    (forall d. Data d => d -> m d) -> Tolerance -> m Tolerance)
-> Data Tolerance
Tolerance -> Constr
Tolerance -> DataType
(forall b. Data b => b -> b) -> Tolerance -> Tolerance
forall a.
Typeable a =>
(forall (c :: * -> *).
 (forall d b. Data d => c (d -> b) -> d -> c b)
 -> (forall g. g -> c g) -> a -> c a)
-> (forall (c :: * -> *).
    (forall b r. Data b => c (b -> r) -> c r)
    -> (forall r. r -> c r) -> Constr -> c a)
-> (a -> Constr)
-> (a -> DataType)
-> (forall (t :: * -> *) (c :: * -> *).
    Typeable t =>
    (forall d. Data d => c (t d)) -> Maybe (c a))
-> (forall (t :: * -> * -> *) (c :: * -> *).
    Typeable t =>
    (forall d e. (Data d, Data e) => c (t d e)) -> Maybe (c a))
-> ((forall b. Data b => b -> b) -> a -> a)
-> (forall r r'.
    (r -> r' -> r) -> r -> (forall d. Data d => d -> r') -> a -> r)
-> (forall r r'.
    (r' -> r -> r) -> r -> (forall d. Data d => d -> r') -> a -> r)
-> (forall u. (forall d. Data d => d -> u) -> a -> [u])
-> (forall u. Int -> (forall d. Data d => d -> u) -> a -> u)
-> (forall (m :: * -> *).
    Monad m =>
    (forall d. Data d => d -> m d) -> a -> m a)
-> (forall (m :: * -> *).
    MonadPlus m =>
    (forall d. Data d => d -> m d) -> a -> m a)
-> (forall (m :: * -> *).
    MonadPlus m =>
    (forall d. Data d => d -> m d) -> a -> m a)
-> Data a
forall u. Int -> (forall d. Data d => d -> u) -> Tolerance -> u
forall u. (forall d. Data d => d -> u) -> Tolerance -> [u]
forall r r'.
(r -> r' -> r)
-> r -> (forall d. Data d => d -> r') -> Tolerance -> r
forall r r'.
(r' -> r -> r)
-> r -> (forall d. Data d => d -> r') -> Tolerance -> r
forall (m :: * -> *).
Monad m =>
(forall d. Data d => d -> m d) -> Tolerance -> m Tolerance
forall (m :: * -> *).
MonadPlus m =>
(forall d. Data d => d -> m d) -> Tolerance -> m Tolerance
forall (c :: * -> *).
(forall b r. Data b => c (b -> r) -> c r)
-> (forall r. r -> c r) -> Constr -> c Tolerance
forall (c :: * -> *).
(forall d b. Data d => c (d -> b) -> d -> c b)
-> (forall g. g -> c g) -> Tolerance -> c Tolerance
forall (t :: * -> *) (c :: * -> *).
Typeable t =>
(forall d. Data d => c (t d)) -> Maybe (c Tolerance)
forall (t :: * -> * -> *) (c :: * -> *).
Typeable t =>
(forall d e. (Data d, Data e) => c (t d e)) -> Maybe (c Tolerance)
$cgfoldl :: forall (c :: * -> *).
(forall d b. Data d => c (d -> b) -> d -> c b)
-> (forall g. g -> c g) -> Tolerance -> c Tolerance
gfoldl :: forall (c :: * -> *).
(forall d b. Data d => c (d -> b) -> d -> c b)
-> (forall g. g -> c g) -> Tolerance -> c Tolerance
$cgunfold :: forall (c :: * -> *).
(forall b r. Data b => c (b -> r) -> c r)
-> (forall r. r -> c r) -> Constr -> c Tolerance
gunfold :: forall (c :: * -> *).
(forall b r. Data b => c (b -> r) -> c r)
-> (forall r. r -> c r) -> Constr -> c Tolerance
$ctoConstr :: Tolerance -> Constr
toConstr :: Tolerance -> Constr
$cdataTypeOf :: Tolerance -> DataType
dataTypeOf :: Tolerance -> DataType
$cdataCast1 :: forall (t :: * -> *) (c :: * -> *).
Typeable t =>
(forall d. Data d => c (t d)) -> Maybe (c Tolerance)
dataCast1 :: forall (t :: * -> *) (c :: * -> *).
Typeable t =>
(forall d. Data d => c (t d)) -> Maybe (c Tolerance)
$cdataCast2 :: forall (t :: * -> * -> *) (c :: * -> *).
Typeable t =>
(forall d e. (Data d, Data e) => c (t d e)) -> Maybe (c Tolerance)
dataCast2 :: forall (t :: * -> * -> *) (c :: * -> *).
Typeable t =>
(forall d e. (Data d, Data e) => c (t d e)) -> Maybe (c Tolerance)
$cgmapT :: (forall b. Data b => b -> b) -> Tolerance -> Tolerance
gmapT :: (forall b. Data b => b -> b) -> Tolerance -> Tolerance
$cgmapQl :: forall r r'.
(r -> r' -> r)
-> r -> (forall d. Data d => d -> r') -> Tolerance -> r
gmapQl :: forall r r'.
(r -> r' -> r)
-> r -> (forall d. Data d => d -> r') -> Tolerance -> r
$cgmapQr :: forall r r'.
(r' -> r -> r)
-> r -> (forall d. Data d => d -> r') -> Tolerance -> r
gmapQr :: forall r r'.
(r' -> r -> r)
-> r -> (forall d. Data d => d -> r') -> Tolerance -> r
$cgmapQ :: forall u. (forall d. Data d => d -> u) -> Tolerance -> [u]
gmapQ :: forall u. (forall d. Data d => d -> u) -> Tolerance -> [u]
$cgmapQi :: forall u. Int -> (forall d. Data d => d -> u) -> Tolerance -> u
gmapQi :: forall u. Int -> (forall d. Data d => d -> u) -> Tolerance -> u
$cgmapM :: forall (m :: * -> *).
Monad m =>
(forall d. Data d => d -> m d) -> Tolerance -> m Tolerance
gmapM :: forall (m :: * -> *).
Monad m =>
(forall d. Data d => d -> m d) -> Tolerance -> m Tolerance
$cgmapMp :: forall (m :: * -> *).
MonadPlus m =>
(forall d. Data d => d -> m d) -> Tolerance -> m Tolerance
gmapMp :: forall (m :: * -> *).
MonadPlus m =>
(forall d. Data d => d -> m d) -> Tolerance -> m Tolerance
$cgmapMo :: forall (m :: * -> *).
MonadPlus m =>
(forall d. Data d => d -> m d) -> Tolerance -> m Tolerance
gmapMo :: forall (m :: * -> *).
MonadPlus m =>
(forall d. Data d => d -> m d) -> Tolerance -> m Tolerance
Data, (forall x. Tolerance -> Rep Tolerance x)
-> (forall x. Rep Tolerance x -> Tolerance) -> Generic Tolerance
forall x. Rep Tolerance x -> Tolerance
forall x. Tolerance -> Rep Tolerance x
forall a.
(forall x. a -> Rep a x) -> (forall x. Rep a x -> a) -> Generic a
$cfrom :: forall x. Tolerance -> Rep Tolerance x
from :: forall x. Tolerance -> Rep Tolerance x
$cto :: forall x. Rep Tolerance x -> Tolerance
to :: forall x. Rep Tolerance x -> Tolerance
Generic)

-- | Check that two values are approximately equal. In addition to
--   specification values are considered equal if they're within 1ulp
--   of precision. No further improvement could be done anyway.
withinTolerance :: Tolerance -> Double -> Double -> Bool
withinTolerance :: Tolerance -> Double -> Double -> Bool
withinTolerance Tolerance
_ Double
a Double
b
  | Int -> Double -> Double -> Bool
within Int
1 Double
a Double
b = Bool
True
withinTolerance (RelTol Double
eps) Double
a Double
b = Double -> Double -> Double -> Bool
eqRelErr Double
eps Double
a Double
b
withinTolerance (AbsTol Double
tol) Double
a Double
b = Double -> Double
forall a. Num a => a -> a
abs (Double
a Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
b) Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
< Double
tol

-- | Type class for checking whether iteration converged already.
class IterationStep a where
  -- | Return @Just root@ is current iteration converged within
  --   required error tolerance. Returns @Nothing@ otherwise.
  matchRoot :: Tolerance -> a -> Maybe (Root Double)

-- | Find root in lazy list of iterations.
findRoot :: IterationStep a
  => Int                        -- ^ Maximum
  -> Tolerance                  -- ^ Error tolerance
  -> [a]
  -> Root Double
findRoot :: forall a. IterationStep a => Int -> Tolerance -> [a] -> Root Double
findRoot Int
maxN Tolerance
tol = Int -> [a] -> Root Double
go Int
0
  where
    go :: Int -> [a] -> Root Double
go !Int
i [a]
_  | Int
i Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
>= Int
maxN = Root Double
forall a. Root a
SearchFailed
    go !Int
_ []             = Root Double
forall a. Root a
SearchFailed
    go  Int
i (a
x:[a]
xs)  = case Tolerance -> a -> Maybe (Root Double)
forall a. IterationStep a => Tolerance -> a -> Maybe (Root Double)
matchRoot Tolerance
tol a
x of
      Just Root Double
r  -> Root Double
r
      Maybe (Root Double)
Nothing -> Int -> [a] -> Root Double
go (Int
iInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
1) [a]
xs
{-# INLINABLE  findRoot #-}
{-# SPECIALIZE findRoot :: Int -> Tolerance -> [RiddersStep] -> Root Double #-}
{-# SPECIALIZE findRoot :: Int -> Tolerance -> [NewtonStep]  -> Root Double #-}


----------------------------------------------------------------
-- Attaching information to roots
----------------------------------------------------------------

-- | Parameters for 'ridders' root finding
data RiddersParam = RiddersParam
  { RiddersParam -> Int
riddersMaxIter :: !Int
    -- ^ Maximum number of iterations. Default = 100
  , RiddersParam -> Tolerance
riddersTol     :: !Tolerance
    -- ^ Error tolerance for root approximation. Default is relative
    --   error 4·ε, where ε is machine precision.
  }
  deriving (RiddersParam -> RiddersParam -> Bool
(RiddersParam -> RiddersParam -> Bool)
-> (RiddersParam -> RiddersParam -> Bool) -> Eq RiddersParam
forall a. (a -> a -> Bool) -> (a -> a -> Bool) -> Eq a
$c== :: RiddersParam -> RiddersParam -> Bool
== :: RiddersParam -> RiddersParam -> Bool
$c/= :: RiddersParam -> RiddersParam -> Bool
/= :: RiddersParam -> RiddersParam -> Bool
Eq, ReadPrec [RiddersParam]
ReadPrec RiddersParam
Int -> ReadS RiddersParam
ReadS [RiddersParam]
(Int -> ReadS RiddersParam)
-> ReadS [RiddersParam]
-> ReadPrec RiddersParam
-> ReadPrec [RiddersParam]
-> Read RiddersParam
forall a.
(Int -> ReadS a)
-> ReadS [a] -> ReadPrec a -> ReadPrec [a] -> Read a
$creadsPrec :: Int -> ReadS RiddersParam
readsPrec :: Int -> ReadS RiddersParam
$creadList :: ReadS [RiddersParam]
readList :: ReadS [RiddersParam]
$creadPrec :: ReadPrec RiddersParam
readPrec :: ReadPrec RiddersParam
$creadListPrec :: ReadPrec [RiddersParam]
readListPrec :: ReadPrec [RiddersParam]
Read, Int -> RiddersParam -> ShowS
[RiddersParam] -> ShowS
RiddersParam -> String
(Int -> RiddersParam -> ShowS)
-> (RiddersParam -> String)
-> ([RiddersParam] -> ShowS)
-> Show RiddersParam
forall a.
(Int -> a -> ShowS) -> (a -> String) -> ([a] -> ShowS) -> Show a
$cshowsPrec :: Int -> RiddersParam -> ShowS
showsPrec :: Int -> RiddersParam -> ShowS
$cshow :: RiddersParam -> String
show :: RiddersParam -> String
$cshowList :: [RiddersParam] -> ShowS
showList :: [RiddersParam] -> ShowS
Show, Typeable, Typeable RiddersParam
Typeable RiddersParam =>
(forall (c :: * -> *).
 (forall d b. Data d => c (d -> b) -> d -> c b)
 -> (forall g. g -> c g) -> RiddersParam -> c RiddersParam)
-> (forall (c :: * -> *).
    (forall b r. Data b => c (b -> r) -> c r)
    -> (forall r. r -> c r) -> Constr -> c RiddersParam)
-> (RiddersParam -> Constr)
-> (RiddersParam -> DataType)
-> (forall (t :: * -> *) (c :: * -> *).
    Typeable t =>
    (forall d. Data d => c (t d)) -> Maybe (c RiddersParam))
-> (forall (t :: * -> * -> *) (c :: * -> *).
    Typeable t =>
    (forall d e. (Data d, Data e) => c (t d e))
    -> Maybe (c RiddersParam))
-> ((forall b. Data b => b -> b) -> RiddersParam -> RiddersParam)
-> (forall r r'.
    (r -> r' -> r)
    -> r -> (forall d. Data d => d -> r') -> RiddersParam -> r)
-> (forall r r'.
    (r' -> r -> r)
    -> r -> (forall d. Data d => d -> r') -> RiddersParam -> r)
-> (forall u. (forall d. Data d => d -> u) -> RiddersParam -> [u])
-> (forall u.
    Int -> (forall d. Data d => d -> u) -> RiddersParam -> u)
-> (forall (m :: * -> *).
    Monad m =>
    (forall d. Data d => d -> m d) -> RiddersParam -> m RiddersParam)
-> (forall (m :: * -> *).
    MonadPlus m =>
    (forall d. Data d => d -> m d) -> RiddersParam -> m RiddersParam)
-> (forall (m :: * -> *).
    MonadPlus m =>
    (forall d. Data d => d -> m d) -> RiddersParam -> m RiddersParam)
-> Data RiddersParam
RiddersParam -> Constr
RiddersParam -> DataType
(forall b. Data b => b -> b) -> RiddersParam -> RiddersParam
forall a.
Typeable a =>
(forall (c :: * -> *).
 (forall d b. Data d => c (d -> b) -> d -> c b)
 -> (forall g. g -> c g) -> a -> c a)
-> (forall (c :: * -> *).
    (forall b r. Data b => c (b -> r) -> c r)
    -> (forall r. r -> c r) -> Constr -> c a)
-> (a -> Constr)
-> (a -> DataType)
-> (forall (t :: * -> *) (c :: * -> *).
    Typeable t =>
    (forall d. Data d => c (t d)) -> Maybe (c a))
-> (forall (t :: * -> * -> *) (c :: * -> *).
    Typeable t =>
    (forall d e. (Data d, Data e) => c (t d e)) -> Maybe (c a))
-> ((forall b. Data b => b -> b) -> a -> a)
-> (forall r r'.
    (r -> r' -> r) -> r -> (forall d. Data d => d -> r') -> a -> r)
-> (forall r r'.
    (r' -> r -> r) -> r -> (forall d. Data d => d -> r') -> a -> r)
-> (forall u. (forall d. Data d => d -> u) -> a -> [u])
-> (forall u. Int -> (forall d. Data d => d -> u) -> a -> u)
-> (forall (m :: * -> *).
    Monad m =>
    (forall d. Data d => d -> m d) -> a -> m a)
-> (forall (m :: * -> *).
    MonadPlus m =>
    (forall d. Data d => d -> m d) -> a -> m a)
-> (forall (m :: * -> *).
    MonadPlus m =>
    (forall d. Data d => d -> m d) -> a -> m a)
-> Data a
forall u. Int -> (forall d. Data d => d -> u) -> RiddersParam -> u
forall u. (forall d. Data d => d -> u) -> RiddersParam -> [u]
forall r r'.
(r -> r' -> r)
-> r -> (forall d. Data d => d -> r') -> RiddersParam -> r
forall r r'.
(r' -> r -> r)
-> r -> (forall d. Data d => d -> r') -> RiddersParam -> r
forall (m :: * -> *).
Monad m =>
(forall d. Data d => d -> m d) -> RiddersParam -> m RiddersParam
forall (m :: * -> *).
MonadPlus m =>
(forall d. Data d => d -> m d) -> RiddersParam -> m RiddersParam
forall (c :: * -> *).
(forall b r. Data b => c (b -> r) -> c r)
-> (forall r. r -> c r) -> Constr -> c RiddersParam
forall (c :: * -> *).
(forall d b. Data d => c (d -> b) -> d -> c b)
-> (forall g. g -> c g) -> RiddersParam -> c RiddersParam
forall (t :: * -> *) (c :: * -> *).
Typeable t =>
(forall d. Data d => c (t d)) -> Maybe (c RiddersParam)
forall (t :: * -> * -> *) (c :: * -> *).
Typeable t =>
(forall d e. (Data d, Data e) => c (t d e))
-> Maybe (c RiddersParam)
$cgfoldl :: forall (c :: * -> *).
(forall d b. Data d => c (d -> b) -> d -> c b)
-> (forall g. g -> c g) -> RiddersParam -> c RiddersParam
gfoldl :: forall (c :: * -> *).
(forall d b. Data d => c (d -> b) -> d -> c b)
-> (forall g. g -> c g) -> RiddersParam -> c RiddersParam
$cgunfold :: forall (c :: * -> *).
(forall b r. Data b => c (b -> r) -> c r)
-> (forall r. r -> c r) -> Constr -> c RiddersParam
gunfold :: forall (c :: * -> *).
(forall b r. Data b => c (b -> r) -> c r)
-> (forall r. r -> c r) -> Constr -> c RiddersParam
$ctoConstr :: RiddersParam -> Constr
toConstr :: RiddersParam -> Constr
$cdataTypeOf :: RiddersParam -> DataType
dataTypeOf :: RiddersParam -> DataType
$cdataCast1 :: forall (t :: * -> *) (c :: * -> *).
Typeable t =>
(forall d. Data d => c (t d)) -> Maybe (c RiddersParam)
dataCast1 :: forall (t :: * -> *) (c :: * -> *).
Typeable t =>
(forall d. Data d => c (t d)) -> Maybe (c RiddersParam)
$cdataCast2 :: forall (t :: * -> * -> *) (c :: * -> *).
Typeable t =>
(forall d e. (Data d, Data e) => c (t d e))
-> Maybe (c RiddersParam)
dataCast2 :: forall (t :: * -> * -> *) (c :: * -> *).
Typeable t =>
(forall d e. (Data d, Data e) => c (t d e))
-> Maybe (c RiddersParam)
$cgmapT :: (forall b. Data b => b -> b) -> RiddersParam -> RiddersParam
gmapT :: (forall b. Data b => b -> b) -> RiddersParam -> RiddersParam
$cgmapQl :: forall r r'.
(r -> r' -> r)
-> r -> (forall d. Data d => d -> r') -> RiddersParam -> r
gmapQl :: forall r r'.
(r -> r' -> r)
-> r -> (forall d. Data d => d -> r') -> RiddersParam -> r
$cgmapQr :: forall r r'.
(r' -> r -> r)
-> r -> (forall d. Data d => d -> r') -> RiddersParam -> r
gmapQr :: forall r r'.
(r' -> r -> r)
-> r -> (forall d. Data d => d -> r') -> RiddersParam -> r
$cgmapQ :: forall u. (forall d. Data d => d -> u) -> RiddersParam -> [u]
gmapQ :: forall u. (forall d. Data d => d -> u) -> RiddersParam -> [u]
$cgmapQi :: forall u. Int -> (forall d. Data d => d -> u) -> RiddersParam -> u
gmapQi :: forall u. Int -> (forall d. Data d => d -> u) -> RiddersParam -> u
$cgmapM :: forall (m :: * -> *).
Monad m =>
(forall d. Data d => d -> m d) -> RiddersParam -> m RiddersParam
gmapM :: forall (m :: * -> *).
Monad m =>
(forall d. Data d => d -> m d) -> RiddersParam -> m RiddersParam
$cgmapMp :: forall (m :: * -> *).
MonadPlus m =>
(forall d. Data d => d -> m d) -> RiddersParam -> m RiddersParam
gmapMp :: forall (m :: * -> *).
MonadPlus m =>
(forall d. Data d => d -> m d) -> RiddersParam -> m RiddersParam
$cgmapMo :: forall (m :: * -> *).
MonadPlus m =>
(forall d. Data d => d -> m d) -> RiddersParam -> m RiddersParam
gmapMo :: forall (m :: * -> *).
MonadPlus m =>
(forall d. Data d => d -> m d) -> RiddersParam -> m RiddersParam
Data, (forall x. RiddersParam -> Rep RiddersParam x)
-> (forall x. Rep RiddersParam x -> RiddersParam)
-> Generic RiddersParam
forall x. Rep RiddersParam x -> RiddersParam
forall x. RiddersParam -> Rep RiddersParam x
forall a.
(forall x. a -> Rep a x) -> (forall x. Rep a x -> a) -> Generic a
$cfrom :: forall x. RiddersParam -> Rep RiddersParam x
from :: forall x. RiddersParam -> Rep RiddersParam x
$cto :: forall x. Rep RiddersParam x -> RiddersParam
to :: forall x. Rep RiddersParam x -> RiddersParam
Generic)

instance Default RiddersParam where
  def :: RiddersParam
def = RiddersParam
        { riddersMaxIter :: Int
riddersMaxIter = Int
100
        , riddersTol :: Tolerance
riddersTol     = Double -> Tolerance
RelTol (Double
4 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
m_epsilon)
        }

-- | Single Ridders step. It's a bracket of root
data RiddersStep
  = RiddersStep   !Double !Double
  -- ^ Ridders step. Parameters are bracket for the root
  | RiddersBisect !Double !Double
  -- ^ Bisection step. It's fallback which is taken when Ridders
  --   update takes us out of bracket
  | RiddersRoot   !Double
  -- ^ Root found
  | RiddersNoBracket
  -- ^ Root is not bracketed
  deriving (RiddersStep -> RiddersStep -> Bool
(RiddersStep -> RiddersStep -> Bool)
-> (RiddersStep -> RiddersStep -> Bool) -> Eq RiddersStep
forall a. (a -> a -> Bool) -> (a -> a -> Bool) -> Eq a
$c== :: RiddersStep -> RiddersStep -> Bool
== :: RiddersStep -> RiddersStep -> Bool
$c/= :: RiddersStep -> RiddersStep -> Bool
/= :: RiddersStep -> RiddersStep -> Bool
Eq, ReadPrec [RiddersStep]
ReadPrec RiddersStep
Int -> ReadS RiddersStep
ReadS [RiddersStep]
(Int -> ReadS RiddersStep)
-> ReadS [RiddersStep]
-> ReadPrec RiddersStep
-> ReadPrec [RiddersStep]
-> Read RiddersStep
forall a.
(Int -> ReadS a)
-> ReadS [a] -> ReadPrec a -> ReadPrec [a] -> Read a
$creadsPrec :: Int -> ReadS RiddersStep
readsPrec :: Int -> ReadS RiddersStep
$creadList :: ReadS [RiddersStep]
readList :: ReadS [RiddersStep]
$creadPrec :: ReadPrec RiddersStep
readPrec :: ReadPrec RiddersStep
$creadListPrec :: ReadPrec [RiddersStep]
readListPrec :: ReadPrec [RiddersStep]
Read, Int -> RiddersStep -> ShowS
[RiddersStep] -> ShowS
RiddersStep -> String
(Int -> RiddersStep -> ShowS)
-> (RiddersStep -> String)
-> ([RiddersStep] -> ShowS)
-> Show RiddersStep
forall a.
(Int -> a -> ShowS) -> (a -> String) -> ([a] -> ShowS) -> Show a
$cshowsPrec :: Int -> RiddersStep -> ShowS
showsPrec :: Int -> RiddersStep -> ShowS
$cshow :: RiddersStep -> String
show :: RiddersStep -> String
$cshowList :: [RiddersStep] -> ShowS
showList :: [RiddersStep] -> ShowS
Show, Typeable, Typeable RiddersStep
Typeable RiddersStep =>
(forall (c :: * -> *).
 (forall d b. Data d => c (d -> b) -> d -> c b)
 -> (forall g. g -> c g) -> RiddersStep -> c RiddersStep)
-> (forall (c :: * -> *).
    (forall b r. Data b => c (b -> r) -> c r)
    -> (forall r. r -> c r) -> Constr -> c RiddersStep)
-> (RiddersStep -> Constr)
-> (RiddersStep -> DataType)
-> (forall (t :: * -> *) (c :: * -> *).
    Typeable t =>
    (forall d. Data d => c (t d)) -> Maybe (c RiddersStep))
-> (forall (t :: * -> * -> *) (c :: * -> *).
    Typeable t =>
    (forall d e. (Data d, Data e) => c (t d e))
    -> Maybe (c RiddersStep))
-> ((forall b. Data b => b -> b) -> RiddersStep -> RiddersStep)
-> (forall r r'.
    (r -> r' -> r)
    -> r -> (forall d. Data d => d -> r') -> RiddersStep -> r)
-> (forall r r'.
    (r' -> r -> r)
    -> r -> (forall d. Data d => d -> r') -> RiddersStep -> r)
-> (forall u. (forall d. Data d => d -> u) -> RiddersStep -> [u])
-> (forall u.
    Int -> (forall d. Data d => d -> u) -> RiddersStep -> u)
-> (forall (m :: * -> *).
    Monad m =>
    (forall d. Data d => d -> m d) -> RiddersStep -> m RiddersStep)
-> (forall (m :: * -> *).
    MonadPlus m =>
    (forall d. Data d => d -> m d) -> RiddersStep -> m RiddersStep)
-> (forall (m :: * -> *).
    MonadPlus m =>
    (forall d. Data d => d -> m d) -> RiddersStep -> m RiddersStep)
-> Data RiddersStep
RiddersStep -> Constr
RiddersStep -> DataType
(forall b. Data b => b -> b) -> RiddersStep -> RiddersStep
forall a.
Typeable a =>
(forall (c :: * -> *).
 (forall d b. Data d => c (d -> b) -> d -> c b)
 -> (forall g. g -> c g) -> a -> c a)
-> (forall (c :: * -> *).
    (forall b r. Data b => c (b -> r) -> c r)
    -> (forall r. r -> c r) -> Constr -> c a)
-> (a -> Constr)
-> (a -> DataType)
-> (forall (t :: * -> *) (c :: * -> *).
    Typeable t =>
    (forall d. Data d => c (t d)) -> Maybe (c a))
-> (forall (t :: * -> * -> *) (c :: * -> *).
    Typeable t =>
    (forall d e. (Data d, Data e) => c (t d e)) -> Maybe (c a))
-> ((forall b. Data b => b -> b) -> a -> a)
-> (forall r r'.
    (r -> r' -> r) -> r -> (forall d. Data d => d -> r') -> a -> r)
-> (forall r r'.
    (r' -> r -> r) -> r -> (forall d. Data d => d -> r') -> a -> r)
-> (forall u. (forall d. Data d => d -> u) -> a -> [u])
-> (forall u. Int -> (forall d. Data d => d -> u) -> a -> u)
-> (forall (m :: * -> *).
    Monad m =>
    (forall d. Data d => d -> m d) -> a -> m a)
-> (forall (m :: * -> *).
    MonadPlus m =>
    (forall d. Data d => d -> m d) -> a -> m a)
-> (forall (m :: * -> *).
    MonadPlus m =>
    (forall d. Data d => d -> m d) -> a -> m a)
-> Data a
forall u. Int -> (forall d. Data d => d -> u) -> RiddersStep -> u
forall u. (forall d. Data d => d -> u) -> RiddersStep -> [u]
forall r r'.
(r -> r' -> r)
-> r -> (forall d. Data d => d -> r') -> RiddersStep -> r
forall r r'.
(r' -> r -> r)
-> r -> (forall d. Data d => d -> r') -> RiddersStep -> r
forall (m :: * -> *).
Monad m =>
(forall d. Data d => d -> m d) -> RiddersStep -> m RiddersStep
forall (m :: * -> *).
MonadPlus m =>
(forall d. Data d => d -> m d) -> RiddersStep -> m RiddersStep
forall (c :: * -> *).
(forall b r. Data b => c (b -> r) -> c r)
-> (forall r. r -> c r) -> Constr -> c RiddersStep
forall (c :: * -> *).
(forall d b. Data d => c (d -> b) -> d -> c b)
-> (forall g. g -> c g) -> RiddersStep -> c RiddersStep
forall (t :: * -> *) (c :: * -> *).
Typeable t =>
(forall d. Data d => c (t d)) -> Maybe (c RiddersStep)
forall (t :: * -> * -> *) (c :: * -> *).
Typeable t =>
(forall d e. (Data d, Data e) => c (t d e))
-> Maybe (c RiddersStep)
$cgfoldl :: forall (c :: * -> *).
(forall d b. Data d => c (d -> b) -> d -> c b)
-> (forall g. g -> c g) -> RiddersStep -> c RiddersStep
gfoldl :: forall (c :: * -> *).
(forall d b. Data d => c (d -> b) -> d -> c b)
-> (forall g. g -> c g) -> RiddersStep -> c RiddersStep
$cgunfold :: forall (c :: * -> *).
(forall b r. Data b => c (b -> r) -> c r)
-> (forall r. r -> c r) -> Constr -> c RiddersStep
gunfold :: forall (c :: * -> *).
(forall b r. Data b => c (b -> r) -> c r)
-> (forall r. r -> c r) -> Constr -> c RiddersStep
$ctoConstr :: RiddersStep -> Constr
toConstr :: RiddersStep -> Constr
$cdataTypeOf :: RiddersStep -> DataType
dataTypeOf :: RiddersStep -> DataType
$cdataCast1 :: forall (t :: * -> *) (c :: * -> *).
Typeable t =>
(forall d. Data d => c (t d)) -> Maybe (c RiddersStep)
dataCast1 :: forall (t :: * -> *) (c :: * -> *).
Typeable t =>
(forall d. Data d => c (t d)) -> Maybe (c RiddersStep)
$cdataCast2 :: forall (t :: * -> * -> *) (c :: * -> *).
Typeable t =>
(forall d e. (Data d, Data e) => c (t d e))
-> Maybe (c RiddersStep)
dataCast2 :: forall (t :: * -> * -> *) (c :: * -> *).
Typeable t =>
(forall d e. (Data d, Data e) => c (t d e))
-> Maybe (c RiddersStep)
$cgmapT :: (forall b. Data b => b -> b) -> RiddersStep -> RiddersStep
gmapT :: (forall b. Data b => b -> b) -> RiddersStep -> RiddersStep
$cgmapQl :: forall r r'.
(r -> r' -> r)
-> r -> (forall d. Data d => d -> r') -> RiddersStep -> r
gmapQl :: forall r r'.
(r -> r' -> r)
-> r -> (forall d. Data d => d -> r') -> RiddersStep -> r
$cgmapQr :: forall r r'.
(r' -> r -> r)
-> r -> (forall d. Data d => d -> r') -> RiddersStep -> r
gmapQr :: forall r r'.
(r' -> r -> r)
-> r -> (forall d. Data d => d -> r') -> RiddersStep -> r
$cgmapQ :: forall u. (forall d. Data d => d -> u) -> RiddersStep -> [u]
gmapQ :: forall u. (forall d. Data d => d -> u) -> RiddersStep -> [u]
$cgmapQi :: forall u. Int -> (forall d. Data d => d -> u) -> RiddersStep -> u
gmapQi :: forall u. Int -> (forall d. Data d => d -> u) -> RiddersStep -> u
$cgmapM :: forall (m :: * -> *).
Monad m =>
(forall d. Data d => d -> m d) -> RiddersStep -> m RiddersStep
gmapM :: forall (m :: * -> *).
Monad m =>
(forall d. Data d => d -> m d) -> RiddersStep -> m RiddersStep
$cgmapMp :: forall (m :: * -> *).
MonadPlus m =>
(forall d. Data d => d -> m d) -> RiddersStep -> m RiddersStep
gmapMp :: forall (m :: * -> *).
MonadPlus m =>
(forall d. Data d => d -> m d) -> RiddersStep -> m RiddersStep
$cgmapMo :: forall (m :: * -> *).
MonadPlus m =>
(forall d. Data d => d -> m d) -> RiddersStep -> m RiddersStep
gmapMo :: forall (m :: * -> *).
MonadPlus m =>
(forall d. Data d => d -> m d) -> RiddersStep -> m RiddersStep
Data, (forall x. RiddersStep -> Rep RiddersStep x)
-> (forall x. Rep RiddersStep x -> RiddersStep)
-> Generic RiddersStep
forall x. Rep RiddersStep x -> RiddersStep
forall x. RiddersStep -> Rep RiddersStep x
forall a.
(forall x. a -> Rep a x) -> (forall x. Rep a x -> a) -> Generic a
$cfrom :: forall x. RiddersStep -> Rep RiddersStep x
from :: forall x. RiddersStep -> Rep RiddersStep x
$cto :: forall x. Rep RiddersStep x -> RiddersStep
to :: forall x. Rep RiddersStep x -> RiddersStep
Generic)

instance NFData RiddersStep where
  rnf :: RiddersStep -> ()
rnf RiddersStep
x = RiddersStep
x RiddersStep -> () -> ()
forall a b. a -> b -> b
`seq` ()

instance IterationStep RiddersStep where
  matchRoot :: Tolerance -> RiddersStep -> Maybe (Root Double)
matchRoot Tolerance
tol RiddersStep
r = case RiddersStep
r of
    RiddersRoot Double
x               -> Root Double -> Maybe (Root Double)
forall a. a -> Maybe a
Just (Root Double -> Maybe (Root Double))
-> Root Double -> Maybe (Root Double)
forall a b. (a -> b) -> a -> b
$ Double -> Root Double
forall a. a -> Root a
Root Double
x
    RiddersStep
RiddersNoBracket            -> Root Double -> Maybe (Root Double)
forall a. a -> Maybe a
Just Root Double
forall a. Root a
NotBracketed
    RiddersStep Double
a Double
b
      | Tolerance -> Double -> Double -> Bool
withinTolerance Tolerance
tol Double
a Double
b -> Root Double -> Maybe (Root Double)
forall a. a -> Maybe a
Just (Root Double -> Maybe (Root Double))
-> Root Double -> Maybe (Root Double)
forall a b. (a -> b) -> a -> b
$ Double -> Root Double
forall a. a -> Root a
Root ((Double
a Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
b) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
2)
      | Bool
otherwise               -> Maybe (Root Double)
forall a. Maybe a
Nothing
    RiddersBisect Double
a Double
b
      | Tolerance -> Double -> Double -> Bool
withinTolerance Tolerance
tol Double
a Double
b -> Root Double -> Maybe (Root Double)
forall a. a -> Maybe a
Just (Root Double -> Maybe (Root Double))
-> Root Double -> Maybe (Root Double)
forall a b. (a -> b) -> a -> b
$ Double -> Root Double
forall a. a -> Root a
Root ((Double
a Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
b) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
2)
      | Bool
otherwise               -> Maybe (Root Double)
forall a. Maybe a
Nothing


-- | Use the method of Ridders[Ridders1979] to compute a root of a
--   function. It doesn't require derivative and provide quadratic
--   convergence (number of significant digits grows quadratically
--   with number of iterations).
--
--   The function must have opposite signs when evaluated at the lower
--   and upper bounds of the search (i.e. the root must be
--   bracketed). If there's more that one root in the bracket
--   iteration will converge to some root in the bracket.
ridders
  :: RiddersParam               -- ^ Parameters for algorithms. @def@
                                --   provides reasonable defaults
  -> (Double,Double)            -- ^ Bracket for root
  -> (Double -> Double)         -- ^ Function to find roots
  -> Root Double
ridders :: RiddersParam
-> (Double, Double) -> (Double -> Double) -> Root Double
ridders RiddersParam
p (Double, Double)
bracket Double -> Double
fun
  = Int -> Tolerance -> [RiddersStep] -> Root Double
forall a. IterationStep a => Int -> Tolerance -> [a] -> Root Double
findRoot (RiddersParam -> Int
riddersMaxIter RiddersParam
p) (RiddersParam -> Tolerance
riddersTol RiddersParam
p)
  ([RiddersStep] -> Root Double) -> [RiddersStep] -> Root Double
forall a b. (a -> b) -> a -> b
$ (Double, Double) -> (Double -> Double) -> [RiddersStep]
riddersIterations (Double, Double)
bracket Double -> Double
fun

-- | List of iterations for Ridders methods. See 'ridders' for
--   documentation of parameters
riddersIterations :: (Double,Double) -> (Double -> Double) -> [RiddersStep]
riddersIterations :: (Double, Double) -> (Double -> Double) -> [RiddersStep]
riddersIterations (Double
lo,Double
hi) Double -> Double
f
  | Double
flo Double -> Double -> Bool
forall a. Eq a => a -> a -> Bool
== Double
0    = [Double -> RiddersStep
RiddersRoot Double
lo]
  | Double
fhi Double -> Double -> Bool
forall a. Eq a => a -> a -> Bool
== Double
0    = [Double -> RiddersStep
RiddersRoot Double
hi]
    -- root is not bracketed
  | Double
floDouble -> Double -> Double
forall a. Num a => a -> a -> a
*Double
fhi Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
> Double
0 = [RiddersStep
RiddersNoBracket]
    -- Ensure that a<b in iterations
  | Double
lo Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
< Double
hi     = Double -> Double -> RiddersStep
RiddersStep Double
lo Double
hi RiddersStep -> [RiddersStep] -> [RiddersStep]
forall a. a -> [a] -> [a]
: Double -> Double -> Double -> Double -> [RiddersStep]
go Double
lo Double
flo Double
hi Double
fhi
  | Bool
otherwise   = Double -> Double -> RiddersStep
RiddersStep Double
lo Double
hi RiddersStep -> [RiddersStep] -> [RiddersStep]
forall a. a -> [a] -> [a]
: Double -> Double -> Double -> Double -> [RiddersStep]
go Double
hi Double
fhi Double
lo Double
flo
  where
    flo :: Double
flo = Double -> Double
f Double
lo
    fhi :: Double
fhi = Double -> Double
f Double
hi
    --
    go :: Double -> Double -> Double -> Double -> [RiddersStep]
go !Double
a !Double
fa !Double
b !Double
fb
      | Double
fm Double -> Double -> Bool
forall a. Eq a => a -> a -> Bool
== Double
0       = [Double -> RiddersStep
RiddersRoot Double
m]
      | Double
fn Double -> Double -> Bool
forall a. Eq a => a -> a -> Bool
== Double
0       = [Double -> RiddersStep
RiddersRoot Double
n]
      -- Ridder's approximation coincide with one of old bounds or
      -- went out of (a,b) range due to numerical problems. Revert
      -- to bisection
      | Double
n Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
<= Double
a Bool -> Bool -> Bool
|| Double
n Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
>= Double
b   = case () of
          ()
_| Double
fmDouble -> Double -> Double
forall a. Num a => a -> a -> a
*Double
fa Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
< Double
0 -> Double -> Double -> Double -> Double -> [RiddersStep]
recBisect Double
a Double
fa Double
m Double
fm
           | Bool
otherwise -> Double -> Double -> Double -> Double -> [RiddersStep]
recBisect Double
m Double
fm Double
b Double
fb
      | Double
fnDouble -> Double -> Double
forall a. Num a => a -> a -> a
*Double
fm Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
< Double
0          = Double -> Double -> Double -> Double -> [RiddersStep]
recRidders Double
n Double
fn Double
m Double
fm
      | Double
fnDouble -> Double -> Double
forall a. Num a => a -> a -> a
*Double
fa Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
< Double
0          = Double -> Double -> Double -> Double -> [RiddersStep]
recRidders Double
a Double
fa Double
n Double
fn
      | Bool
otherwise          = Double -> Double -> Double -> Double -> [RiddersStep]
recRidders Double
n Double
fn Double
b Double
fb
      where
        recBisect :: Double -> Double -> Double -> Double -> [RiddersStep]
recBisect  Double
x Double
fx Double
y Double
fy = Double -> Double -> RiddersStep
RiddersBisect Double
x Double
y RiddersStep -> [RiddersStep] -> [RiddersStep]
forall a. a -> [a] -> [a]
: Double -> Double -> Double -> Double -> [RiddersStep]
go Double
x Double
fx Double
y Double
fy
        recRidders :: Double -> Double -> Double -> Double -> [RiddersStep]
recRidders Double
x Double
fx Double
y Double
fy = Double -> Double -> RiddersStep
RiddersStep   Double
x Double
y RiddersStep -> [RiddersStep] -> [RiddersStep]
forall a. a -> [a] -> [a]
: Double -> Double -> Double -> Double -> [RiddersStep]
go Double
x Double
fx Double
y Double
fy
        --
        dm :: Double
dm  = (Double
b Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
a) Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
0.5
        -- Mean point
        m :: Double
m   = (Double
a Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
b) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
2
        fm :: Double
fm  = Double -> Double
f Double
m
        -- Ridders update
        n :: Double
n   = Double
m Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double -> Double
forall a. Num a => a -> a
signum (Double
fb Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
fa) Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
dm Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
fm Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double -> Double
forall a. Floating a => a -> a
sqrt(Double
fmDouble -> Double -> Double
forall a. Num a => a -> a -> a
*Double
fm Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
faDouble -> Double -> Double
forall a. Num a => a -> a -> a
*Double
fb)
        fn :: Double
fn  = Double -> Double
f Double
n



----------------------------------------------------------------
-- Newton-Raphson algorithm
----------------------------------------------------------------

-- | Parameters for 'ridders' root finding
data NewtonParam = NewtonParam
  { NewtonParam -> Int
newtonMaxIter :: !Int
    -- ^ Maximum number of iterations. Default = 50
  , NewtonParam -> Tolerance
newtonTol     :: !Tolerance
    -- ^ Error tolerance for root approximation. Default is relative
    --   error 4·ε, where ε is machine precision
  }
  deriving (NewtonParam -> NewtonParam -> Bool
(NewtonParam -> NewtonParam -> Bool)
-> (NewtonParam -> NewtonParam -> Bool) -> Eq NewtonParam
forall a. (a -> a -> Bool) -> (a -> a -> Bool) -> Eq a
$c== :: NewtonParam -> NewtonParam -> Bool
== :: NewtonParam -> NewtonParam -> Bool
$c/= :: NewtonParam -> NewtonParam -> Bool
/= :: NewtonParam -> NewtonParam -> Bool
Eq, ReadPrec [NewtonParam]
ReadPrec NewtonParam
Int -> ReadS NewtonParam
ReadS [NewtonParam]
(Int -> ReadS NewtonParam)
-> ReadS [NewtonParam]
-> ReadPrec NewtonParam
-> ReadPrec [NewtonParam]
-> Read NewtonParam
forall a.
(Int -> ReadS a)
-> ReadS [a] -> ReadPrec a -> ReadPrec [a] -> Read a
$creadsPrec :: Int -> ReadS NewtonParam
readsPrec :: Int -> ReadS NewtonParam
$creadList :: ReadS [NewtonParam]
readList :: ReadS [NewtonParam]
$creadPrec :: ReadPrec NewtonParam
readPrec :: ReadPrec NewtonParam
$creadListPrec :: ReadPrec [NewtonParam]
readListPrec :: ReadPrec [NewtonParam]
Read, Int -> NewtonParam -> ShowS
[NewtonParam] -> ShowS
NewtonParam -> String
(Int -> NewtonParam -> ShowS)
-> (NewtonParam -> String)
-> ([NewtonParam] -> ShowS)
-> Show NewtonParam
forall a.
(Int -> a -> ShowS) -> (a -> String) -> ([a] -> ShowS) -> Show a
$cshowsPrec :: Int -> NewtonParam -> ShowS
showsPrec :: Int -> NewtonParam -> ShowS
$cshow :: NewtonParam -> String
show :: NewtonParam -> String
$cshowList :: [NewtonParam] -> ShowS
showList :: [NewtonParam] -> ShowS
Show, Typeable, Typeable NewtonParam
Typeable NewtonParam =>
(forall (c :: * -> *).
 (forall d b. Data d => c (d -> b) -> d -> c b)
 -> (forall g. g -> c g) -> NewtonParam -> c NewtonParam)
-> (forall (c :: * -> *).
    (forall b r. Data b => c (b -> r) -> c r)
    -> (forall r. r -> c r) -> Constr -> c NewtonParam)
-> (NewtonParam -> Constr)
-> (NewtonParam -> DataType)
-> (forall (t :: * -> *) (c :: * -> *).
    Typeable t =>
    (forall d. Data d => c (t d)) -> Maybe (c NewtonParam))
-> (forall (t :: * -> * -> *) (c :: * -> *).
    Typeable t =>
    (forall d e. (Data d, Data e) => c (t d e))
    -> Maybe (c NewtonParam))
-> ((forall b. Data b => b -> b) -> NewtonParam -> NewtonParam)
-> (forall r r'.
    (r -> r' -> r)
    -> r -> (forall d. Data d => d -> r') -> NewtonParam -> r)
-> (forall r r'.
    (r' -> r -> r)
    -> r -> (forall d. Data d => d -> r') -> NewtonParam -> r)
-> (forall u. (forall d. Data d => d -> u) -> NewtonParam -> [u])
-> (forall u.
    Int -> (forall d. Data d => d -> u) -> NewtonParam -> u)
-> (forall (m :: * -> *).
    Monad m =>
    (forall d. Data d => d -> m d) -> NewtonParam -> m NewtonParam)
-> (forall (m :: * -> *).
    MonadPlus m =>
    (forall d. Data d => d -> m d) -> NewtonParam -> m NewtonParam)
-> (forall (m :: * -> *).
    MonadPlus m =>
    (forall d. Data d => d -> m d) -> NewtonParam -> m NewtonParam)
-> Data NewtonParam
NewtonParam -> Constr
NewtonParam -> DataType
(forall b. Data b => b -> b) -> NewtonParam -> NewtonParam
forall a.
Typeable a =>
(forall (c :: * -> *).
 (forall d b. Data d => c (d -> b) -> d -> c b)
 -> (forall g. g -> c g) -> a -> c a)
-> (forall (c :: * -> *).
    (forall b r. Data b => c (b -> r) -> c r)
    -> (forall r. r -> c r) -> Constr -> c a)
-> (a -> Constr)
-> (a -> DataType)
-> (forall (t :: * -> *) (c :: * -> *).
    Typeable t =>
    (forall d. Data d => c (t d)) -> Maybe (c a))
-> (forall (t :: * -> * -> *) (c :: * -> *).
    Typeable t =>
    (forall d e. (Data d, Data e) => c (t d e)) -> Maybe (c a))
-> ((forall b. Data b => b -> b) -> a -> a)
-> (forall r r'.
    (r -> r' -> r) -> r -> (forall d. Data d => d -> r') -> a -> r)
-> (forall r r'.
    (r' -> r -> r) -> r -> (forall d. Data d => d -> r') -> a -> r)
-> (forall u. (forall d. Data d => d -> u) -> a -> [u])
-> (forall u. Int -> (forall d. Data d => d -> u) -> a -> u)
-> (forall (m :: * -> *).
    Monad m =>
    (forall d. Data d => d -> m d) -> a -> m a)
-> (forall (m :: * -> *).
    MonadPlus m =>
    (forall d. Data d => d -> m d) -> a -> m a)
-> (forall (m :: * -> *).
    MonadPlus m =>
    (forall d. Data d => d -> m d) -> a -> m a)
-> Data a
forall u. Int -> (forall d. Data d => d -> u) -> NewtonParam -> u
forall u. (forall d. Data d => d -> u) -> NewtonParam -> [u]
forall r r'.
(r -> r' -> r)
-> r -> (forall d. Data d => d -> r') -> NewtonParam -> r
forall r r'.
(r' -> r -> r)
-> r -> (forall d. Data d => d -> r') -> NewtonParam -> r
forall (m :: * -> *).
Monad m =>
(forall d. Data d => d -> m d) -> NewtonParam -> m NewtonParam
forall (m :: * -> *).
MonadPlus m =>
(forall d. Data d => d -> m d) -> NewtonParam -> m NewtonParam
forall (c :: * -> *).
(forall b r. Data b => c (b -> r) -> c r)
-> (forall r. r -> c r) -> Constr -> c NewtonParam
forall (c :: * -> *).
(forall d b. Data d => c (d -> b) -> d -> c b)
-> (forall g. g -> c g) -> NewtonParam -> c NewtonParam
forall (t :: * -> *) (c :: * -> *).
Typeable t =>
(forall d. Data d => c (t d)) -> Maybe (c NewtonParam)
forall (t :: * -> * -> *) (c :: * -> *).
Typeable t =>
(forall d e. (Data d, Data e) => c (t d e))
-> Maybe (c NewtonParam)
$cgfoldl :: forall (c :: * -> *).
(forall d b. Data d => c (d -> b) -> d -> c b)
-> (forall g. g -> c g) -> NewtonParam -> c NewtonParam
gfoldl :: forall (c :: * -> *).
(forall d b. Data d => c (d -> b) -> d -> c b)
-> (forall g. g -> c g) -> NewtonParam -> c NewtonParam
$cgunfold :: forall (c :: * -> *).
(forall b r. Data b => c (b -> r) -> c r)
-> (forall r. r -> c r) -> Constr -> c NewtonParam
gunfold :: forall (c :: * -> *).
(forall b r. Data b => c (b -> r) -> c r)
-> (forall r. r -> c r) -> Constr -> c NewtonParam
$ctoConstr :: NewtonParam -> Constr
toConstr :: NewtonParam -> Constr
$cdataTypeOf :: NewtonParam -> DataType
dataTypeOf :: NewtonParam -> DataType
$cdataCast1 :: forall (t :: * -> *) (c :: * -> *).
Typeable t =>
(forall d. Data d => c (t d)) -> Maybe (c NewtonParam)
dataCast1 :: forall (t :: * -> *) (c :: * -> *).
Typeable t =>
(forall d. Data d => c (t d)) -> Maybe (c NewtonParam)
$cdataCast2 :: forall (t :: * -> * -> *) (c :: * -> *).
Typeable t =>
(forall d e. (Data d, Data e) => c (t d e))
-> Maybe (c NewtonParam)
dataCast2 :: forall (t :: * -> * -> *) (c :: * -> *).
Typeable t =>
(forall d e. (Data d, Data e) => c (t d e))
-> Maybe (c NewtonParam)
$cgmapT :: (forall b. Data b => b -> b) -> NewtonParam -> NewtonParam
gmapT :: (forall b. Data b => b -> b) -> NewtonParam -> NewtonParam
$cgmapQl :: forall r r'.
(r -> r' -> r)
-> r -> (forall d. Data d => d -> r') -> NewtonParam -> r
gmapQl :: forall r r'.
(r -> r' -> r)
-> r -> (forall d. Data d => d -> r') -> NewtonParam -> r
$cgmapQr :: forall r r'.
(r' -> r -> r)
-> r -> (forall d. Data d => d -> r') -> NewtonParam -> r
gmapQr :: forall r r'.
(r' -> r -> r)
-> r -> (forall d. Data d => d -> r') -> NewtonParam -> r
$cgmapQ :: forall u. (forall d. Data d => d -> u) -> NewtonParam -> [u]
gmapQ :: forall u. (forall d. Data d => d -> u) -> NewtonParam -> [u]
$cgmapQi :: forall u. Int -> (forall d. Data d => d -> u) -> NewtonParam -> u
gmapQi :: forall u. Int -> (forall d. Data d => d -> u) -> NewtonParam -> u
$cgmapM :: forall (m :: * -> *).
Monad m =>
(forall d. Data d => d -> m d) -> NewtonParam -> m NewtonParam
gmapM :: forall (m :: * -> *).
Monad m =>
(forall d. Data d => d -> m d) -> NewtonParam -> m NewtonParam
$cgmapMp :: forall (m :: * -> *).
MonadPlus m =>
(forall d. Data d => d -> m d) -> NewtonParam -> m NewtonParam
gmapMp :: forall (m :: * -> *).
MonadPlus m =>
(forall d. Data d => d -> m d) -> NewtonParam -> m NewtonParam
$cgmapMo :: forall (m :: * -> *).
MonadPlus m =>
(forall d. Data d => d -> m d) -> NewtonParam -> m NewtonParam
gmapMo :: forall (m :: * -> *).
MonadPlus m =>
(forall d. Data d => d -> m d) -> NewtonParam -> m NewtonParam
Data, (forall x. NewtonParam -> Rep NewtonParam x)
-> (forall x. Rep NewtonParam x -> NewtonParam)
-> Generic NewtonParam
forall x. Rep NewtonParam x -> NewtonParam
forall x. NewtonParam -> Rep NewtonParam x
forall a.
(forall x. a -> Rep a x) -> (forall x. Rep a x -> a) -> Generic a
$cfrom :: forall x. NewtonParam -> Rep NewtonParam x
from :: forall x. NewtonParam -> Rep NewtonParam x
$cto :: forall x. Rep NewtonParam x -> NewtonParam
to :: forall x. Rep NewtonParam x -> NewtonParam
Generic)

instance Default NewtonParam where
  def :: NewtonParam
def = NewtonParam
        { newtonMaxIter :: Int
newtonMaxIter = Int
50
        , newtonTol :: Tolerance
newtonTol     = Double -> Tolerance
RelTol (Double
4 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
m_epsilon)
        }

-- | Steps for Newton iterations
data NewtonStep
  = NewtonStep         !Double !Double
  -- ^ Normal Newton-Raphson update. Parameters are: old guess, new guess
  | NewtonBisection    !Double !Double
  -- ^ Bisection fallback when Newton-Raphson iteration doesn't
  --   work. Parameters are bracket on root
  | NewtonRoot         !Double
  -- ^ Root is found
  | NewtonNoBracket
  -- ^ Root is not bracketed
  deriving (NewtonStep -> NewtonStep -> Bool
(NewtonStep -> NewtonStep -> Bool)
-> (NewtonStep -> NewtonStep -> Bool) -> Eq NewtonStep
forall a. (a -> a -> Bool) -> (a -> a -> Bool) -> Eq a
$c== :: NewtonStep -> NewtonStep -> Bool
== :: NewtonStep -> NewtonStep -> Bool
$c/= :: NewtonStep -> NewtonStep -> Bool
/= :: NewtonStep -> NewtonStep -> Bool
Eq, ReadPrec [NewtonStep]
ReadPrec NewtonStep
Int -> ReadS NewtonStep
ReadS [NewtonStep]
(Int -> ReadS NewtonStep)
-> ReadS [NewtonStep]
-> ReadPrec NewtonStep
-> ReadPrec [NewtonStep]
-> Read NewtonStep
forall a.
(Int -> ReadS a)
-> ReadS [a] -> ReadPrec a -> ReadPrec [a] -> Read a
$creadsPrec :: Int -> ReadS NewtonStep
readsPrec :: Int -> ReadS NewtonStep
$creadList :: ReadS [NewtonStep]
readList :: ReadS [NewtonStep]
$creadPrec :: ReadPrec NewtonStep
readPrec :: ReadPrec NewtonStep
$creadListPrec :: ReadPrec [NewtonStep]
readListPrec :: ReadPrec [NewtonStep]
Read, Int -> NewtonStep -> ShowS
[NewtonStep] -> ShowS
NewtonStep -> String
(Int -> NewtonStep -> ShowS)
-> (NewtonStep -> String)
-> ([NewtonStep] -> ShowS)
-> Show NewtonStep
forall a.
(Int -> a -> ShowS) -> (a -> String) -> ([a] -> ShowS) -> Show a
$cshowsPrec :: Int -> NewtonStep -> ShowS
showsPrec :: Int -> NewtonStep -> ShowS
$cshow :: NewtonStep -> String
show :: NewtonStep -> String
$cshowList :: [NewtonStep] -> ShowS
showList :: [NewtonStep] -> ShowS
Show, Typeable, Typeable NewtonStep
Typeable NewtonStep =>
(forall (c :: * -> *).
 (forall d b. Data d => c (d -> b) -> d -> c b)
 -> (forall g. g -> c g) -> NewtonStep -> c NewtonStep)
-> (forall (c :: * -> *).
    (forall b r. Data b => c (b -> r) -> c r)
    -> (forall r. r -> c r) -> Constr -> c NewtonStep)
-> (NewtonStep -> Constr)
-> (NewtonStep -> DataType)
-> (forall (t :: * -> *) (c :: * -> *).
    Typeable t =>
    (forall d. Data d => c (t d)) -> Maybe (c NewtonStep))
-> (forall (t :: * -> * -> *) (c :: * -> *).
    Typeable t =>
    (forall d e. (Data d, Data e) => c (t d e))
    -> Maybe (c NewtonStep))
-> ((forall b. Data b => b -> b) -> NewtonStep -> NewtonStep)
-> (forall r r'.
    (r -> r' -> r)
    -> r -> (forall d. Data d => d -> r') -> NewtonStep -> r)
-> (forall r r'.
    (r' -> r -> r)
    -> r -> (forall d. Data d => d -> r') -> NewtonStep -> r)
-> (forall u. (forall d. Data d => d -> u) -> NewtonStep -> [u])
-> (forall u.
    Int -> (forall d. Data d => d -> u) -> NewtonStep -> u)
-> (forall (m :: * -> *).
    Monad m =>
    (forall d. Data d => d -> m d) -> NewtonStep -> m NewtonStep)
-> (forall (m :: * -> *).
    MonadPlus m =>
    (forall d. Data d => d -> m d) -> NewtonStep -> m NewtonStep)
-> (forall (m :: * -> *).
    MonadPlus m =>
    (forall d. Data d => d -> m d) -> NewtonStep -> m NewtonStep)
-> Data NewtonStep
NewtonStep -> Constr
NewtonStep -> DataType
(forall b. Data b => b -> b) -> NewtonStep -> NewtonStep
forall a.
Typeable a =>
(forall (c :: * -> *).
 (forall d b. Data d => c (d -> b) -> d -> c b)
 -> (forall g. g -> c g) -> a -> c a)
-> (forall (c :: * -> *).
    (forall b r. Data b => c (b -> r) -> c r)
    -> (forall r. r -> c r) -> Constr -> c a)
-> (a -> Constr)
-> (a -> DataType)
-> (forall (t :: * -> *) (c :: * -> *).
    Typeable t =>
    (forall d. Data d => c (t d)) -> Maybe (c a))
-> (forall (t :: * -> * -> *) (c :: * -> *).
    Typeable t =>
    (forall d e. (Data d, Data e) => c (t d e)) -> Maybe (c a))
-> ((forall b. Data b => b -> b) -> a -> a)
-> (forall r r'.
    (r -> r' -> r) -> r -> (forall d. Data d => d -> r') -> a -> r)
-> (forall r r'.
    (r' -> r -> r) -> r -> (forall d. Data d => d -> r') -> a -> r)
-> (forall u. (forall d. Data d => d -> u) -> a -> [u])
-> (forall u. Int -> (forall d. Data d => d -> u) -> a -> u)
-> (forall (m :: * -> *).
    Monad m =>
    (forall d. Data d => d -> m d) -> a -> m a)
-> (forall (m :: * -> *).
    MonadPlus m =>
    (forall d. Data d => d -> m d) -> a -> m a)
-> (forall (m :: * -> *).
    MonadPlus m =>
    (forall d. Data d => d -> m d) -> a -> m a)
-> Data a
forall u. Int -> (forall d. Data d => d -> u) -> NewtonStep -> u
forall u. (forall d. Data d => d -> u) -> NewtonStep -> [u]
forall r r'.
(r -> r' -> r)
-> r -> (forall d. Data d => d -> r') -> NewtonStep -> r
forall r r'.
(r' -> r -> r)
-> r -> (forall d. Data d => d -> r') -> NewtonStep -> r
forall (m :: * -> *).
Monad m =>
(forall d. Data d => d -> m d) -> NewtonStep -> m NewtonStep
forall (m :: * -> *).
MonadPlus m =>
(forall d. Data d => d -> m d) -> NewtonStep -> m NewtonStep
forall (c :: * -> *).
(forall b r. Data b => c (b -> r) -> c r)
-> (forall r. r -> c r) -> Constr -> c NewtonStep
forall (c :: * -> *).
(forall d b. Data d => c (d -> b) -> d -> c b)
-> (forall g. g -> c g) -> NewtonStep -> c NewtonStep
forall (t :: * -> *) (c :: * -> *).
Typeable t =>
(forall d. Data d => c (t d)) -> Maybe (c NewtonStep)
forall (t :: * -> * -> *) (c :: * -> *).
Typeable t =>
(forall d e. (Data d, Data e) => c (t d e)) -> Maybe (c NewtonStep)
$cgfoldl :: forall (c :: * -> *).
(forall d b. Data d => c (d -> b) -> d -> c b)
-> (forall g. g -> c g) -> NewtonStep -> c NewtonStep
gfoldl :: forall (c :: * -> *).
(forall d b. Data d => c (d -> b) -> d -> c b)
-> (forall g. g -> c g) -> NewtonStep -> c NewtonStep
$cgunfold :: forall (c :: * -> *).
(forall b r. Data b => c (b -> r) -> c r)
-> (forall r. r -> c r) -> Constr -> c NewtonStep
gunfold :: forall (c :: * -> *).
(forall b r. Data b => c (b -> r) -> c r)
-> (forall r. r -> c r) -> Constr -> c NewtonStep
$ctoConstr :: NewtonStep -> Constr
toConstr :: NewtonStep -> Constr
$cdataTypeOf :: NewtonStep -> DataType
dataTypeOf :: NewtonStep -> DataType
$cdataCast1 :: forall (t :: * -> *) (c :: * -> *).
Typeable t =>
(forall d. Data d => c (t d)) -> Maybe (c NewtonStep)
dataCast1 :: forall (t :: * -> *) (c :: * -> *).
Typeable t =>
(forall d. Data d => c (t d)) -> Maybe (c NewtonStep)
$cdataCast2 :: forall (t :: * -> * -> *) (c :: * -> *).
Typeable t =>
(forall d e. (Data d, Data e) => c (t d e)) -> Maybe (c NewtonStep)
dataCast2 :: forall (t :: * -> * -> *) (c :: * -> *).
Typeable t =>
(forall d e. (Data d, Data e) => c (t d e)) -> Maybe (c NewtonStep)
$cgmapT :: (forall b. Data b => b -> b) -> NewtonStep -> NewtonStep
gmapT :: (forall b. Data b => b -> b) -> NewtonStep -> NewtonStep
$cgmapQl :: forall r r'.
(r -> r' -> r)
-> r -> (forall d. Data d => d -> r') -> NewtonStep -> r
gmapQl :: forall r r'.
(r -> r' -> r)
-> r -> (forall d. Data d => d -> r') -> NewtonStep -> r
$cgmapQr :: forall r r'.
(r' -> r -> r)
-> r -> (forall d. Data d => d -> r') -> NewtonStep -> r
gmapQr :: forall r r'.
(r' -> r -> r)
-> r -> (forall d. Data d => d -> r') -> NewtonStep -> r
$cgmapQ :: forall u. (forall d. Data d => d -> u) -> NewtonStep -> [u]
gmapQ :: forall u. (forall d. Data d => d -> u) -> NewtonStep -> [u]
$cgmapQi :: forall u. Int -> (forall d. Data d => d -> u) -> NewtonStep -> u
gmapQi :: forall u. Int -> (forall d. Data d => d -> u) -> NewtonStep -> u
$cgmapM :: forall (m :: * -> *).
Monad m =>
(forall d. Data d => d -> m d) -> NewtonStep -> m NewtonStep
gmapM :: forall (m :: * -> *).
Monad m =>
(forall d. Data d => d -> m d) -> NewtonStep -> m NewtonStep
$cgmapMp :: forall (m :: * -> *).
MonadPlus m =>
(forall d. Data d => d -> m d) -> NewtonStep -> m NewtonStep
gmapMp :: forall (m :: * -> *).
MonadPlus m =>
(forall d. Data d => d -> m d) -> NewtonStep -> m NewtonStep
$cgmapMo :: forall (m :: * -> *).
MonadPlus m =>
(forall d. Data d => d -> m d) -> NewtonStep -> m NewtonStep
gmapMo :: forall (m :: * -> *).
MonadPlus m =>
(forall d. Data d => d -> m d) -> NewtonStep -> m NewtonStep
Data, (forall x. NewtonStep -> Rep NewtonStep x)
-> (forall x. Rep NewtonStep x -> NewtonStep) -> Generic NewtonStep
forall x. Rep NewtonStep x -> NewtonStep
forall x. NewtonStep -> Rep NewtonStep x
forall a.
(forall x. a -> Rep a x) -> (forall x. Rep a x -> a) -> Generic a
$cfrom :: forall x. NewtonStep -> Rep NewtonStep x
from :: forall x. NewtonStep -> Rep NewtonStep x
$cto :: forall x. Rep NewtonStep x -> NewtonStep
to :: forall x. Rep NewtonStep x -> NewtonStep
Generic)

instance NFData NewtonStep where
  rnf :: NewtonStep -> ()
rnf NewtonStep
x = NewtonStep
x NewtonStep -> () -> ()
forall a b. a -> b -> b
`seq` ()

instance IterationStep NewtonStep where
  matchRoot :: Tolerance -> NewtonStep -> Maybe (Root Double)
matchRoot Tolerance
tol NewtonStep
r = case NewtonStep
r of
    NewtonRoot Double
x                 -> Root Double -> Maybe (Root Double)
forall a. a -> Maybe a
Just (Double -> Root Double
forall a. a -> Root a
Root Double
x)
    NewtonStep
NewtonNoBracket              -> Root Double -> Maybe (Root Double)
forall a. a -> Maybe a
Just Root Double
forall a. Root a
NotBracketed
    NewtonStep Double
x Double
x'
      | Tolerance -> Double -> Double -> Bool
withinTolerance Tolerance
tol Double
x Double
x' -> Root Double -> Maybe (Root Double)
forall a. a -> Maybe a
Just (Double -> Root Double
forall a. a -> Root a
Root Double
x')
      | Bool
otherwise                -> Maybe (Root Double)
forall a. Maybe a
Nothing
    NewtonBisection Double
a Double
b
      | Tolerance -> Double -> Double -> Bool
withinTolerance Tolerance
tol Double
a Double
b  -> Root Double -> Maybe (Root Double)
forall a. a -> Maybe a
Just (Double -> Root Double
forall a. a -> Root a
Root ((Double
a Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
b) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
2))
      | Bool
otherwise                -> Maybe (Root Double)
forall a. Maybe a
Nothing
  {-# INLINE matchRoot #-}


-- | Solve equation using Newton-Raphson iterations.
--
--   This method require both initial guess and bounds for root. If
--   Newton step takes us out of bounds on root function reverts to
--   bisection.
newtonRaphson
  :: NewtonParam                 -- ^ Parameters for algorithm. @def@
                                 --   provide reasonable defaults.
  -> (Double,Double,Double)      -- ^ Triple of @(low bound, initial
                                 --   guess, upper bound)@. If initial
                                 --   guess if out of bracket middle
                                 --   of bracket is taken as
                                 --   approximation
  -> (Double -> (Double,Double)) -- ^ Function to find root of. It
                                 --   returns pair of function value and
                                 --   its first derivative
  -> Root Double
newtonRaphson :: NewtonParam
-> (Double, Double, Double)
-> (Double -> (Double, Double))
-> Root Double
newtonRaphson NewtonParam
p (Double, Double, Double)
guess Double -> (Double, Double)
fun
  = Int -> Tolerance -> [NewtonStep] -> Root Double
forall a. IterationStep a => Int -> Tolerance -> [a] -> Root Double
findRoot (NewtonParam -> Int
newtonMaxIter NewtonParam
p) (NewtonParam -> Tolerance
newtonTol NewtonParam
p)
  ([NewtonStep] -> Root Double) -> [NewtonStep] -> Root Double
forall a b. (a -> b) -> a -> b
$ (Double, Double, Double)
-> (Double -> (Double, Double)) -> [NewtonStep]
newtonRaphsonIterations (Double, Double, Double)
guess Double -> (Double, Double)
fun

-- | List of iteration for Newton-Raphson algorithm. See documentation
--   for 'newtonRaphson' for meaning of parameters.
newtonRaphsonIterations :: (Double,Double,Double) -> (Double -> (Double,Double)) -> [NewtonStep]
newtonRaphsonIterations :: (Double, Double, Double)
-> (Double -> (Double, Double)) -> [NewtonStep]
newtonRaphsonIterations (Double
lo,Double
guess,Double
hi) Double -> (Double, Double)
function
  | Double
flo Double -> Double -> Bool
forall a. Eq a => a -> a -> Bool
== Double
0    = [Double -> NewtonStep
NewtonRoot Double
lo]
  | Double
fhi Double -> Double -> Bool
forall a. Eq a => a -> a -> Bool
== Double
0    = [Double -> NewtonStep
NewtonRoot Double
hi]
  | Double
floDouble -> Double -> Double
forall a. Num a => a -> a -> a
*Double
fhi Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
> Double
0 = [NewtonStep
NewtonNoBracket]
    -- Ensure that function value on low bound is negative
  | Double
flo Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
> Double
0     = Double -> Double -> Double -> [NewtonStep]
go Double
hi Double
guess' Double
lo
  | Bool
otherwise   = Double -> Double -> Double -> [NewtonStep]
go Double
lo Double
guess Double
hi
  where
    (Double
flo,Double
_) = Double -> (Double, Double)
function Double
lo
    (Double
fhi,Double
_) = Double -> (Double, Double)
function Double
hi
    -- Ensure that initial guess is within bracket
    guess' :: Double
guess'
      | Double
guess Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
>= Double
lo Bool -> Bool -> Bool
&& Double
guess Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
<= Double
hi = Double
guess
      | Double
guess Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
>= Double
hi Bool -> Bool -> Bool
&& Double
guess Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
<= Double
lo = Double
guess
      | Bool
otherwise                  = (Double
lo Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
hi) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
2
    -- Newton iterations. Invariant:
    --   > f xA < 0
    --   > f xB > 0
    go :: Double -> Double -> Double -> [NewtonStep]
go Double
xA Double
x Double
xB
      | Double
f  Double -> Double -> Bool
forall a. Eq a => a -> a -> Bool
== Double
0   = [Double -> NewtonStep
NewtonRoot Double
x]
      | Double
f' Double -> Double -> Bool
forall a. Eq a => a -> a -> Bool
== Double
0   = [NewtonStep]
bisectionStep
      -- Accept Newton step since it stays within bracket.
      | (Double
x' Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
xA) Double -> Double -> Double
forall a. Num a => a -> a -> a
* (Double
x' Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
xB) Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
< Double
0 = [NewtonStep]
newtonStep
      -- Otherwise bracket root and pick new approximation as
      -- midpoint.
      | Bool
otherwise                 = [NewtonStep]
bisectionStep
      where
        -- Calculate Newton step
        (Double
f,Double
f') = Double -> (Double, Double)
function Double
x
        x' :: Double
x'   = Double
x Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
f Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
f'
        -- Newton step
        newtonStep :: [NewtonStep]
newtonStep
          | Double
f Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
> Double
0     = Double -> Double -> NewtonStep
NewtonStep Double
x Double
x' NewtonStep -> [NewtonStep] -> [NewtonStep]
forall a. a -> [a] -> [a]
: Double -> Double -> Double -> [NewtonStep]
go Double
xA Double
x' Double
x
          | Bool
otherwise = Double -> Double -> NewtonStep
NewtonStep Double
x Double
x' NewtonStep -> [NewtonStep] -> [NewtonStep]
forall a. a -> [a] -> [a]
: Double -> Double -> Double -> [NewtonStep]
go Double
x  Double
x' Double
xB
        -- Fallback bisection step
        bisectionStep :: [NewtonStep]
bisectionStep
          | Double
f Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
> Double
0     = Double -> Double -> NewtonStep
NewtonBisection Double
xA Double
x NewtonStep -> [NewtonStep] -> [NewtonStep]
forall a. a -> [a] -> [a]
: Double -> Double -> Double -> [NewtonStep]
go Double
xA ((Double
xA Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
x) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
2) Double
x
          | Bool
otherwise = Double -> Double -> NewtonStep
NewtonBisection Double
x Double
xB NewtonStep -> [NewtonStep] -> [NewtonStep]
forall a. a -> [a] -> [a]
: Double -> Double -> Double -> [NewtonStep]
go Double
x  ((Double
x Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
xB) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
2) Double
xB



----------------------------------------------------------------
-- Internal functions
----------------------------------------------------------------

-- $references
--
-- * Ridders, C.F.J. (1979) A new algorithm for computing a single
--   root of a real continuous function.
--   /IEEE Transactions on Circuits and Systems/ 26:979&#8211;980.
--
-- * Press W.H.; Teukolsky S.A.; Vetterling W.T.; Flannery B.P.
--   (2007). \"Section 9.2.1. Ridders' Method\". /Numerical Recipes: The
--   Art of Scientific Computing (3rd ed.)./ New York: Cambridge
--   University Press. ISBN 978-0-521-88068-8.