X-Git-Url: http://gitweb.michael.orlitzky.com/?a=blobdiff_plain;f=src%2FRoots%2FFast.hs;h=8b69786379218b12448b114d4183ca9c3ef64c5d;hb=613b25028368b3651bea0c801e05d9962ad9b604;hp=f7cde8a82a70f3f98605ab87fc94fae5a59147dc;hpb=04f56a8882bb0c574b603f8c3fed9481ea934f7f;p=numerical-analysis.git diff --git a/src/Roots/Fast.hs b/src/Roots/Fast.hs index f7cde8a..8b69786 100644 --- a/src/Roots/Fast.hs +++ b/src/Roots/Fast.hs @@ -9,6 +9,7 @@ module Roots.Fast where import Data.List (find) +import Data.Maybe (fromMaybe) import Normed @@ -30,33 +31,20 @@ has_root :: (RealField.C a, -> 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 b epsilon f_of_a f_of_b + | (signum (f_of_a')) * (signum (f_of_b')) /= 1 = True + | (b - a) <= epsilon' = False + | otherwise = (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 + epsilon' = fromMaybe (b-a) epsilon -- 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 + f_of_a' = fromMaybe (f a) f_of_a + f_of_b' = fromMaybe (f b) f_of_b c = (a + b)/2 @@ -87,13 +75,8 @@ bisect f a b epsilon f_of_a f_of_b 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 + f_of_a' = fromMaybe (f a) f_of_a + f_of_b' = fromMaybe (f b) f_of_b c = (a + b) / 2 @@ -118,26 +101,21 @@ trisect f a b epsilon f_of_a f_of_b | f_of_b' == 0 = Just b | otherwise = -- Use a 'prime' just for consistency. - let (a', b', fa', fb') = - if (has_root f d b (Just epsilon) (Just f_of_d') (Just f_of_b')) - then (d, b, f_of_d', f_of_b') - else - if (has_root f c d (Just epsilon) (Just f_of_c') (Just f_of_d')) - then (c, d, f_of_c', f_of_d') - else (a, c, f_of_a', f_of_c') + let (a', b', fa', fb') + | has_root f d b (Just epsilon) (Just f_of_d') (Just f_of_b') = + (d, b, f_of_d', f_of_b') + | has_root f c d (Just epsilon) (Just f_of_c') (Just f_of_d') = + (c, d, f_of_c', f_of_d') + | otherwise = + (a, c, f_of_a', f_of_c') in if (b-a) < 2*epsilon then Just ((b+a)/2) else trisect f a' b' epsilon (Just fa') (Just fb') 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 + f_of_a' = fromMaybe (f a) f_of_a + f_of_b' = fromMaybe (f b) f_of_b c = (2*a + b) / 3 @@ -153,8 +131,8 @@ trisect f a b epsilon f_of_a f_of_b fixed_point_iterations :: (a -> a) -- ^ The function @f@ to iterate. -> a -- ^ The initial value @x0@. -> [a] -- ^ The resulting sequence of x_{n}. -fixed_point_iterations f x0 = - iterate f x0 +fixed_point_iterations = + iterate -- | Find a fixed point of the function @f@ with the search starting