summaryrefslogtreecommitdiff
diff options
context:
space:
mode:
authorJamesCook <>2011-03-20 18:54:56 (GMT)
committerLuite Stegeman <luite@luite.com>2011-03-20 18:54:56 (GMT)
commit826f600284576a29fe1d68635ddadb12668aa261 (patch)
treeb67d3478bbd7fcc66f0ed54b975cdbddbabe41a5
parentb245cd924549f5780a829de89c999da9ddce8b41 (diff)
version 0.1.1.00.1.1.0
-rw-r--r--roots.cabal6
-rw-r--r--src/Math/Root/Finder.hs125
-rw-r--r--src/Math/Root/Finder/Bisection.hs17
-rw-r--r--src/Math/Root/Finder/Brent.hs6
4 files changed, 123 insertions, 31 deletions
diff --git a/roots.cabal b/roots.cabal
index b9cf1a8..8fb3846 100644
--- 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 @@ tested-with: GHC == 6.8.3,
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
index ed5cdf6..5d6c83d 100644
--- 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 @@ class RootFinder r a b where
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 @@ class RootFinder r a b where
-- 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 @@ traceRoot f a b xacc = go nSteps start (stepRootFinder f start)
-- 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
+ 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
- Tagged nSteps = (const :: Tagged a b -> a -> Tagged a b) defaultNSteps start
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 :: RealFloat a => a
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))
+ ]
+ \ No newline at end of file
diff --git a/src/Math/Root/Finder/Bisection.hs b/src/Math/Root/Finder/Bisection.hs
index a511ad1..37a2940 100644
--- a/src/Math/Root/Finder/Bisection.hs
+++ b/src/Math/Root/Finder/Bisection.hs
@@ -17,18 +17,19 @@ instance (Fractional a, Ord b, Num b) => RootFinder Bisect a b where
| 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
index 0acd5a8..51720cd 100644
--- a/src/Math/Root/Finder/Brent.hs
+++ b/src/Math/Root/Finder/Brent.hs
@@ -38,7 +38,9 @@ instance (RealFloat a, Real b, Fractional b) => RootFinder Brent a b where
| 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 @@ instance (RealFloat a, Real b, Fractional b) => RootFinder Brent a b where
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