From: Michael Orlitzky Date: Tue, 11 Sep 2012 12:40:18 +0000 (-0400) Subject: Split Roots into Roots.Simple and Roots.Fast. X-Git-Url: https://gitweb.michael.orlitzky.com/?a=commitdiff_plain;h=61dcb661ee35ff693a9d64d4d21a65791a480a52;p=numerical-analysis.git Split Roots into Roots.Simple and Roots.Fast. Allow the makefile to find sources in subdirectories. Update module exports in the cabal file. --- diff --git a/makefile b/makefile index 8b9c250..0ffe8d3 100644 --- a/makefile +++ b/makefile @@ -1,23 +1,27 @@ -BIN = dist/build/numerical-analysis/numerical-analysis +# There's only one '$' in the awk script, but we have to double-money +# it for make. +PN = $(shell grep 'name:' *.cabal | awk '{ print $$2 }') +BIN = dist/build/$(PN)/$(PN) DOCTESTS_BIN = dist/build/doctests/doctests +SRCS = $(shell find src/ -name '*.hs') .PHONY : test publish_doc doc dist hlint -$(BIN): src/*.hs +$(BIN): $(SRCS) runghc Setup.hs clean runghc Setup.hs configure --user --flags=${FLAGS} runghc Setup.hs build -$(DOCTESTS_BIN): src/*.hs test/Doctests.hs +$(DOCTESTS_BIN): $(SRCS) test/Doctests.hs runghc Setup.hs configure --user --flags=${FLAGS} --enable-tests runghc Setup.hs build -profile: src/*.hs +profile: $(SRCS) runghc Setup.hs configure --user --enable-executable-profiling runghc Setup.hs build -hpc: src/*.hs +hpc: $(SRCS) FLAGS="hpc" make clean: @@ -34,8 +38,8 @@ dist: runghc Setup.hs configure runghc Setup.hs sdist -# Neither 'haddock' nor 'hscolour' seem to work properly. -doc: src/*.hs + +doc: $(SRCS) runghc Setup.hs configure --user --flags=${FLAGS} runghc Setup.hs hscolour --executables runghc Setup.hs haddock --internal \ diff --git a/numerical-analysis.cabal b/numerical-analysis.cabal index 708644f..7bd0fb4 100644 --- a/numerical-analysis.cabal +++ b/numerical-analysis.cabal @@ -16,7 +16,7 @@ description: data-files: makefile library - exposed-modules: Roots + exposed-modules: Roots.Simple, Roots.Fast build-depends: base == 4.5.*, diff --git a/src/Roots/Fast.hs b/src/Roots/Fast.hs new file mode 100644 index 0000000..8e49750 --- /dev/null +++ b/src/Roots/Fast.hs @@ -0,0 +1,82 @@ +-- | The Roots.Fast module contains faster implementations of the +-- 'Roots.Simple' algorithms. Generally, we will pass precomputed +-- values to the next iteration of a function rather than passing +-- the function and the points at which to (re)evaluate it. + +module Roots.Fast +where + +has_root :: (Fractional a, Ord a, Ord b, Num b) + => (a -> b) -- ^ The function @f@ + -> a -- ^ The \"left\" endpoint, @a@ + -> a -- ^ The \"right\" endpoint, @b@ + -> Maybe a -- ^ The size of the smallest subinterval + -- we'll examine, @epsilon@ + -> Maybe b -- ^ Precoumpted f(a) + -> Maybe b -- ^ Precoumpted f(b) + -> Bool +has_root f a b epsilon f_of_a f_of_b = + if not ((signum (f_of_a')) * (signum (f_of_b')) == 1) then + -- We don't care about epsilon here, there's definitely a root! + True + else + if (b - a) <= epsilon' then + -- Give up, return false. + False + else + -- If either [a,c] or [c,b] have roots, we do too. + (has_root f a c (Just epsilon') (Just f_of_a') Nothing) || + (has_root f c b (Just epsilon') Nothing (Just f_of_b')) + where + -- If the size of the smallest subinterval is not specified, + -- assume we just want to check once on all of [a,b]. + epsilon' = case epsilon of + Nothing -> (b-a) + Just eps -> eps + + -- Compute f(a) and f(b) only if needed. + f_of_a' = case f_of_a of + Nothing -> f a + Just v -> v + + f_of_b' = case f_of_b of + Nothing -> f b + Just v -> v + + c = (a + b)/2 + + + +bisect :: (Fractional a, Ord a, Num b, Ord b) + => (a -> b) -- ^ The function @f@ whose root we seek + -> a -- ^ The \"left\" endpoint of the interval, @a@ + -> a -- ^ The \"right\" endpoint of the interval, @b@ + -> a -- ^ The tolerance, @epsilon@ + -> Maybe b -- ^ Precomputed f(a) + -> Maybe b -- ^ Precomputed f(b) + -> Maybe a +bisect f a b epsilon f_of_a f_of_b + -- We pass @epsilon@ to the 'has_root' function because if we want a + -- result within epsilon of the true root, we need to know that + -- there *is* a root within an interval of length epsilon. + | not (has_root f a b (Just epsilon) (Just f_of_a') (Just f_of_b')) = Nothing + | f_of_a' == 0 = Just a + | f_of_b' == 0 = Just b + | (b - c) < epsilon = Just c + | otherwise = + -- Use a 'prime' just for consistency. + let f_of_c' = f c in + if (has_root f a c (Just epsilon) (Just f_of_a') (Just f_of_c')) + then bisect f a c epsilon (Just f_of_a') (Just f_of_c') + else bisect f c b epsilon (Just f_of_c') (Just f_of_b') + where + -- Compute f(a) and f(b) only if needed. + f_of_a' = case f_of_a of + Nothing -> f a + Just v -> v + + f_of_b' = case f_of_b of + Nothing -> f b + Just v -> v + + c = (a + b) / 2 diff --git a/src/Roots.hs b/src/Roots/Simple.hs similarity index 61% rename from src/Roots.hs rename to src/Roots/Simple.hs index a8c3c76..64124a8 100644 --- a/src/Roots.hs +++ b/src/Roots/Simple.hs @@ -1,14 +1,16 @@ --- | The Roots module contains root-finding algorithms. That is, --- procedures to (numerically) find solutions to the equation, +-- | The Roots.Simple module contains root-finding algorithms. That +-- is, procedures to (numerically) find solutions to the equation, -- -- > f(x) = 0 -- -- where f is assumed to be continuous on the interval of interest. -- -module Roots +module Roots.Simple where +import qualified Roots.Fast as F + -- | Does the (continuous) function @f@ have a root on the interval -- [a,b]? If f(a) <] 0 and f(b) ]> 0, we know that there's a root in @@ -37,23 +39,8 @@ has_root :: (Fractional a, Ord a, Ord b, Num b) -- we'll examine, @epsilon@ -> Bool has_root f a b epsilon = - if not ((signum (f a)) * (signum (f b)) == 1) then - -- We don't care about epsilon here, there's definitely a root! - True - else - if (b - a) <= epsilon' then - -- Give up, return false. - False - else - -- If either [a,c] or [c,b] have roots, we do too. - (has_root f a c (Just epsilon')) || (has_root f c b (Just epsilon')) - where - -- If the size of the smallest subinterval is not specified, - -- assume we just want to check once on all of [a,b]. - epsilon' = case epsilon of - Nothing -> (b-a) - Just eps -> eps - c = (a + b)/2 + F.has_root f a b epsilon Nothing Nothing + @@ -82,16 +69,5 @@ bisect :: (Fractional a, Ord a, Num b, Ord b) -> a -- ^ The \"right\" endpoint of the interval, @b@ -> a -- ^ The tolerance, @epsilon@ -> Maybe a -bisect f a b epsilon - -- We pass @epsilon@ to the 'has_root' function because if we want a - -- result within epsilon of the true root, we need to know that - -- there *is* a root within an interval of length epsilon. - | not (has_root f a b (Just epsilon)) = Nothing - | f a == 0 = Just a - | f b == 0 = Just b - | (b - c) < epsilon = Just c - | otherwise = - if (has_root f a c (Just epsilon)) then bisect f a c epsilon - else bisect f c b epsilon - where - c = (a + b) / 2 +bisect f a b epsilon = + F.bisect f a b epsilon Nothing Nothing diff --git a/test/Doctests.hs b/test/Doctests.hs index f9bacc5..3ce1b4a 100644 --- a/test/Doctests.hs +++ b/test/Doctests.hs @@ -4,4 +4,6 @@ where import Test.DocTest main :: IO () -main = doctest ["--optghc=-isrc", "src/Root.hs"] +main = doctest ["--optghc=-isrc", + "src/Roots/Simple.hs", + "src/Roots/Fast.hs"]