roots 0.1 → 0.1.1.0
raw patch · 4 files changed
+123/−31 lines, 4 files
Files
- roots.cabal +3/−3
- src/Math/Root/Finder.hs +106/−19
- src/Math/Root/Finder/Bisection.hs +9/−8
- src/Math/Root/Finder/Brent.hs +5/−1
roots.cabal view
@@ -1,5 +1,5 @@ name: roots-version: 0.1+version: 0.1.1.0 stability: experimental cabal-version: >= 1.6@@ -20,8 +20,8 @@ GHC == 6.12.1, GHC == 6.12.3 source-repository head- type: darcs- location: http://code.haskell.org/~mokus/roots+ type: git+ location: git://github.com/mokus0/roots.git Library ghc-options: -Wall
src/Math/Root/Finder.hs view
@@ -1,5 +1,13 @@ {-# LANGUAGE MultiParamTypeClasses, ScopedTypeVariables, FlexibleContexts #-}-module Math.Root.Finder where+module Math.Root.Finder+ ( RootFinder(..)+ , getDefaultNSteps+ , runRootFinder+ , traceRoot+ , findRoot, findRootN+ , eps+ , realFloatDefaultNSteps+ ) where import Control.Monad.Instances () import Data.Tagged@@ -37,9 +45,61 @@ defaultNSteps :: Tagged (r a b) Int defaultNSteps = Tagged 250 +-- |Convenience function to access 'defaultNSteps' for a root finder, +-- which requires a little bit of type-gymnastics.+-- +-- This function does not evaluate its argument.+getDefaultNSteps :: RootFinder r a b => r a b -> Int+getDefaultNSteps rf = nSteps+ where+ Tagged nSteps = + (const :: Tagged a b -> a -> Tagged a b)+ defaultNSteps rf++-- |General-purpose driver for stepping a root finder. Given a \"control\"+-- function, the function being searched, and an initial 'RootFinder' state,+-- @runRootFinder step f state@ repeatedly steps the root-finder and passes+-- each intermediate state, along with a count of steps taken, to @step@.+-- +-- The @step@ funtion will be called with the following arguments:+--+-- [@ n :: 'Int' @] +-- The number of steps taken thus far+-- +-- [@ currentState :: r a b @]+-- The current state of the root finder+--+-- [@ continue :: c @]+-- The result of the \"rest\" of the iteration+--+-- For example, the following function simply iterates a root finder+-- and returns every intermediate state (similar to 'traceRoot'):+-- +-- > iterateRoot :: RootFinder r a b => (a -> b) -> a -> a -> [r a b]+-- > iterateRoot f a b = runRootFinder (const (:)) f (initRootFinder f a b)+--+-- And the following function simply iterates the root finder to +-- convergence or throws an error after a given number of steps:+--+-- > solve :: (RootFinder r a b, RealFloat a)+-- > => Int -> (a -> b) -> a -> a -> r a b+-- > solve maxN f a b = runRootFinder step f (initRootFinder f a b)+-- > where+-- > step n x continue+-- > | converged eps x = x+-- > | n > maxN = error "solve: step limit exceeded"+-- > | otherwise = continue+-- +runRootFinder :: (RootFinder r a b) =>+ (Int -> r a b -> c -> c) -> (a -> b) -> r a b -> c+runRootFinder cons f = go 0+ where+ go n x = n `seq` cons n x (go (n+1) (stepRootFinder f x))+ -- |@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.+-- steps it, returning each step of the process in a list. No step limit+-- is imposed.+-- -- 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@@ -48,20 +108,20 @@ -- 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)+traceRoot f a b mbEps = runRootFinder cons 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)+ cons _n x rest = x : if done x rest then [] else rest+ + -- 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.+ done = case mbEps of+ Nothing -> \x (next:_) -> x == next+ Just xacc -> \x _rest -> converged xacc x -- |@findRoot f x0 x1 eps@ initializes a root finder and repeatedly -- steps it. When the algorithm converges to @eps@ or the 'defaultNSteps'@@ -70,15 +130,23 @@ -- 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+findRoot f a b xacc = result where- Tagged nSteps = (const :: Tagged a b -> a -> Tagged a b) defaultNSteps start+ result = findRootN nSteps f a b xacc+ nSteps = getDefaultNSteps (either id id result)++-- |Like 'findRoot' but with a specified limit on the number of steps (rather+-- than using 'defaultNSteps').+findRootN :: (RootFinder r a b, Num a, Ord a) =>+ Int -> (a -> b) -> a -> a -> a -> Either (r a b) (r a b)+findRootN nSteps f a b xacc = runRootFinder step f start+ where start = initRootFinder f a b - go n x+ step n x continue | converged xacc x = Right x- | n <= 0 = Left x- | otherwise = go (n-1) (stepRootFinder f x)+ | n > nSteps = Left x+ | otherwise = continue -- |A useful constant: 'eps' is (for most 'RealFloat' types) the smallest -- positive number such that @1 + eps /= 1@.@@ -86,3 +154,22 @@ eps = eps' where eps' = encodeFloat 1 (1 - floatDigits eps')++-- |For 'RealFloat' types, computes a suitable default step limit based+-- on the precision of the type and a margin of error.+realFloatDefaultNSteps :: RealFloat a => Float -> Tagged (r a b) Int+realFloatDefaultNSteps margin = nSteps+ where+ f :: (Int -> Tagged (r a b) Int) -> (a -> Int) -> a -> Tagged (r a b) Int+ f = (.)+ + nSteps :: RealFloat a => Tagged (r a b) Int+ nSteps = f Tagged n 0+ + n :: RealFloat a => a -> Int+ n x = round $ product+ [ margin+ , realToFrac (floatDigits x)+ , logBase 2 (realToFrac (floatRadix x))+ ]+
src/Math/Root/Finder/Bisection.hs view
@@ -17,18 +17,19 @@ | 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+ stepRootFinder _ orig@(Bisect _ _ 0) = orig+ stepRootFinder f orig@(Bisect x fx dx)+ | x == xMid = orig+ | otherwise = case fMid `compare` 0 of LT -> Bisect xMid fMid dx2- EQ -> orig+ EQ -> Bisect xMid fMid dx2 GT -> Bisect x fx dx2 - where- dx2 = dx * 0.5- xMid = x + dx2- fMid = f xMid+ 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.
src/Math/Root/Finder/Brent.hs view
@@ -38,7 +38,9 @@ | otherwise = advance m (abs (b - a)) where -- Minimum step size to continue with inverse-quadratic interpolation- tol1 = eps * (abs b + 0.5)+ -- This should not be too low; if it is, convergence can be+ -- spectacularly slow+ tol1 = 1e-3 * (abs b + 0.5) abs_s = abs s -- midpoint for bisection step@@ -73,6 +75,8 @@ converged _ Brent{brFB = 0} = True converged tol Brent{brB = b, brE = e} = abs e <= 4 * eps * abs b + tol+ + defaultNSteps = realFloatDefaultNSteps 5 -- |@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