diff --git a/roots.cabal b/roots.cabal --- a/roots.cabal +++ b/roots.cabal @@ -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 diff --git a/src/Math/Root/Finder.hs b/src/Math/Root/Finder.hs --- a/src/Math/Root/Finder.hs +++ b/src/Math/Root/Finder.hs @@ -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)) + ] + diff --git a/src/Math/Root/Finder/Bisection.hs b/src/Math/Root/Finder/Bisection.hs --- a/src/Math/Root/Finder/Bisection.hs +++ b/src/Math/Root/Finder/Bisection.hs @@ -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. diff --git a/src/Math/Root/Finder/Brent.hs b/src/Math/Root/Finder/Brent.hs --- a/src/Math/Root/Finder/Brent.hs +++ b/src/Math/Root/Finder/Brent.hs @@ -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