src/Roots/Fast.hs

   1 {-# LANGUAGE RebindableSyntax #-}
   2
   3 -- | The Roots.Fast module contains faster implementations of the
   4 --   'Roots.Simple' algorithms. Generally, we will pass precomputed
   5 --   values to the next iteration of a function rather than passing
   6 --   the function and the points at which to (re)evaluate it.
   7
   8 module Roots.Fast
   9 where
  10
  11 import Data.List (find)
  12
  13 import Normed
  14
  15 import NumericPrelude hiding (abs)
  16 import qualified Algebra.Absolute as Absolute
  17 import qualified Algebra.Field as Field
  18 import qualified Algebra.RealRing as RealRing
  19 import qualified Algebra.RealField as RealField
  20
  21 has_root :: (RealField.C a,
  22              RealRing.C b,
  23              Absolute.C b)
  24          => (a -> b) -- ^ The function @f@
  25          -> a       -- ^ The \"left\" endpoint, @a@
  26          -> a       -- ^ The \"right\" endpoint, @b@
  27          -> Maybe a -- ^ The size of the smallest subinterval
  28                    --   we'll examine, @epsilon@
  29          -> Maybe b -- ^ Precoumpted f(a)
  30          -> Maybe b -- ^ Precoumpted f(b)
  31          -> Bool
  32 has_root f a b epsilon f_of_a f_of_b =
  33   if not ((signum (f_of_a')) * (signum (f_of_b')) == 1) then
  34     -- We don't care about epsilon here, there's definitely a root!
  35     True
  36   else
  37     if (b - a) <= epsilon' then
  38       -- Give up, return false.
  39       False
  40     else
  41       -- If either [a,c] or [c,b] have roots, we do too.
  42       (has_root f a c (Just epsilon') (Just f_of_a') Nothing) ||
  43         (has_root f c b (Just epsilon') Nothing (Just f_of_b'))
  44   where
  45     -- If the size of the smallest subinterval is not specified,
  46     -- assume we just want to check once on all of [a,b].
  47     epsilon' = case epsilon of
  48                  Nothing -> (b-a)
  49                  Just eps -> eps
  50
  51     -- Compute f(a) and f(b) only if needed.
  52     f_of_a'  = case f_of_a of
  53                  Nothing -> f a
  54                  Just v -> v
  55
  56     f_of_b'  = case f_of_b of
  57                  Nothing -> f b
  58                  Just v -> v
  59
  60     c = (a + b)/2
  61
  62
  63 bisect :: (RealField.C a,
  64            RealRing.C b,
  65            Absolute.C b)
  66        => (a -> b) -- ^ The function @f@ whose root we seek
  67        -> a       -- ^ The \"left\" endpoint of the interval, @a@
  68        -> a       -- ^ The \"right\" endpoint of the interval, @b@
  69        -> a       -- ^ The tolerance, @epsilon@
  70        -> Maybe b -- ^ Precomputed f(a)
  71        -> Maybe b -- ^ Precomputed f(b)
  72        -> Maybe a
  73 bisect f a b epsilon f_of_a f_of_b
  74   -- We pass @epsilon@ to the 'has_root' function because if we want a
  75   -- result within epsilon of the true root, we need to know that
  76   -- there *is* a root within an interval of length epsilon.
  77   | not (has_root f a b (Just epsilon) (Just f_of_a') (Just f_of_b')) = Nothing
  78   | f_of_a' == 0 = Just a
  79   | f_of_b' == 0 = Just b
  80   | (b - c) < epsilon = Just c
  81   | otherwise =
  82       -- Use a 'prime' just for consistency.
  83       let f_of_c' = f c in
  84         if (has_root f a c (Just epsilon) (Just f_of_a') (Just f_of_c'))
  85         then bisect f a c epsilon (Just f_of_a') (Just f_of_c')
  86         else bisect f c b epsilon (Just f_of_c') (Just f_of_b')
  87   where
  88     -- Compute f(a) and f(b) only if needed.
  89     f_of_a'  = case f_of_a of
  90                  Nothing -> f a
  91                  Just v -> v
  92
  93     f_of_b'  = case f_of_b of
  94                  Nothing -> f b
  95                  Just v -> v
  96
  97     c = (a + b) / 2
  98
  99
 100
 101
 102 -- | Iterate the function @f@ with the initial guess @x0@ in hopes of
 103 --   finding a fixed point.
 104 fixed_point_iterations :: (a -> a) -- ^ The function @f@ to iterate.
 105                        -> a       -- ^ The initial value @x0@.
 106                        -> [a]     -- ^ The resulting sequence of x_{n}.
 107 fixed_point_iterations f x0 =
 108   iterate f x0
 109
 110
 111 -- | Find a fixed point of the function @f@ with the search starting
 112 --   at x0. This will find the first element in the chain f(x0),
 113 --   f(f(x0)),... such that the magnitude of the difference between it
 114 --   and the next element is less than epsilon.
 115 --
 116 --   We also return the number of iterations required.
 117 --
 118 fixed_point_with_iterations :: (Normed a,
 119                                 Field.C b,
 120                                 Absolute.C b,
 121                                 Ord b)
 122                             => (a -> a)  -- ^ The function @f@ to iterate.
 123                             -> b        -- ^ The tolerance, @epsilon@.
 124                             -> a        -- ^ The initial value @x0@.
 125                             -> (Int, a) -- ^ The (iterations, fixed point) pair
 126 fixed_point_with_iterations f epsilon x0 =
 127   (fst winning_pair)
 128   where
 129     xn = fixed_point_iterations f x0
 130     xn_plus_one = tail xn
 131
 132     abs_diff v w = norm (v - w)
 133
 134     -- The nth entry in this list is the absolute value of x_{n} -
 135     -- x_{n+1}.
 136     differences = zipWith abs_diff xn xn_plus_one
 137
 138     -- This produces the list [(n, xn)] so that we can determine
 139     -- the number of iterations required.
 140     numbered_xn = zip [0..] xn
 141
 142     -- A list of pairs, (xn, |x_{n} - x_{n+1}|).
 143     pairs = zip numbered_xn differences
 144
 145     -- The pair (xn, |x_{n} - x_{n+1}|) with
 146     -- |x_{n} - x_{n+1}| < epsilon. The pattern match on 'Just' is
 147     -- "safe" since the list is infinite. We'll succeed or loop
 148     -- forever.
 149     Just winning_pair = find (\(_, diff) -> diff < epsilon) pairs