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 import Data.Maybe (fromMaybe)
  13
  14 import Normed
  15
  16 import NumericPrelude hiding (abs)
  17 import qualified Algebra.Absolute as Absolute
  18 import qualified Algebra.Additive as Additive
  19 import qualified Algebra.Algebraic as Algebraic
  20 import qualified Algebra.RealRing as RealRing
  21 import qualified Algebra.RealField as RealField
  22
  23 has_root :: (RealField.C a,
  24              RealRing.C b,
  25              Absolute.C b)
  26          => (a -> b) -- ^ The function @f@
  27          -> a       -- ^ The \"left\" endpoint, @a@
  28          -> a       -- ^ The \"right\" endpoint, @b@
  29          -> Maybe a -- ^ The size of the smallest subinterval
  30                    --   we'll examine, @epsilon@
  31          -> Maybe b -- ^ Precoumpted f(a)
  32          -> Maybe b -- ^ Precoumpted f(b)
  33          -> Bool
  34 has_root f a b epsilon f_of_a f_of_b
  35   | (signum (f_of_a')) * (signum (f_of_b')) /= 1 = True
  36   | (b - a) <= epsilon' = False
  37   | otherwise =
  38       (has_root f a c (Just epsilon') (Just f_of_a') Nothing) ||
  39         (has_root f c b (Just epsilon') Nothing (Just f_of_b'))
  40   where
  41     -- If the size of the smallest subinterval is not specified,
  42     -- assume we just want to check once on all of [a,b].
  43     epsilon' = fromMaybe (b-a) epsilon
  44
  45     -- Compute f(a) and f(b) only if needed.
  46     f_of_a'  = fromMaybe (f a) f_of_a
  47     f_of_b'  = fromMaybe (f b) f_of_b
  48
  49     c = (a + b)/2
  50
  51
  52 bisect :: (RealField.C a,
  53            RealRing.C b,
  54            Absolute.C b)
  55        => (a -> b) -- ^ The function @f@ whose root we seek
  56        -> a       -- ^ The \"left\" endpoint of the interval, @a@
  57        -> a       -- ^ The \"right\" endpoint of the interval, @b@
  58        -> a       -- ^ The tolerance, @epsilon@
  59        -> Maybe b -- ^ Precomputed f(a)
  60        -> Maybe b -- ^ Precomputed f(b)
  61        -> Maybe a
  62 bisect f a b epsilon f_of_a f_of_b
  63   -- We pass @epsilon@ to the 'has_root' function because if we want a
  64   -- result within epsilon of the true root, we need to know that
  65   -- there *is* a root within an interval of length epsilon.
  66   | not (has_root f a b (Just epsilon) (Just f_of_a') (Just f_of_b')) = Nothing
  67   | f_of_a' == 0 = Just a
  68   | f_of_b' == 0 = Just b
  69   | (b - c) < epsilon = Just c
  70   | otherwise =
  71       -- Use a 'prime' just for consistency.
  72       let f_of_c' = f c in
  73         if (has_root f a c (Just epsilon) (Just f_of_a') (Just f_of_c'))
  74         then bisect f a c epsilon (Just f_of_a') (Just f_of_c')
  75         else bisect f c b epsilon (Just f_of_c') (Just f_of_b')
  76   where
  77     -- Compute f(a) and f(b) only if needed.
  78     f_of_a'  = fromMaybe (f a) f_of_a
  79     f_of_b'  = fromMaybe (f b) f_of_b
  80
  81     c = (a + b) / 2
  82
  83
  84
  85 trisect :: (RealField.C a,
  86             RealRing.C b,
  87             Absolute.C b)
  88         => (a -> b) -- ^ The function @f@ whose root we seek
  89         -> a       -- ^ The \"left\" endpoint of the interval, @a@
  90         -> a       -- ^ The \"right\" endpoint of the interval, @b@
  91         -> a       -- ^ The tolerance, @epsilon@
  92         -> Maybe b -- ^ Precomputed f(a)
  93         -> Maybe b -- ^ Precomputed f(b)
  94         -> Maybe a
  95 trisect f a b epsilon f_of_a f_of_b
  96   -- We pass @epsilon@ to the 'has_root' function because if we want a
  97   -- result within epsilon of the true root, we need to know that
  98   -- there *is* a root within an interval of length epsilon.
  99   | not (has_root f a b (Just epsilon) (Just f_of_a') (Just f_of_b')) = Nothing
 100   | f_of_a' == 0 = Just a
 101   | f_of_b' == 0 = Just b
 102   | otherwise =
 103       -- Use a 'prime' just for consistency.
 104     let (a', b', fa', fb')
 105           | has_root f d b (Just epsilon) (Just f_of_d') (Just f_of_b') =
 106               (d, b, f_of_d', f_of_b')
 107           | has_root f c d (Just epsilon) (Just f_of_c') (Just f_of_d') =
 108               (c, d, f_of_c', f_of_d')
 109           | otherwise =
 110               (a, c, f_of_a', f_of_c')
 111     in
 112       if (b-a) < 2*epsilon
 113       then Just ((b+a)/2)
 114       else trisect f a' b' epsilon (Just fa') (Just fb')
 115   where
 116     -- Compute f(a) and f(b) only if needed.
 117     f_of_a'  = fromMaybe (f a) f_of_a
 118     f_of_b'  = fromMaybe (f b) f_of_b
 119
 120     c = (2*a + b) / 3
 121
 122     d = (a + 2*b) / 3
 123
 124     f_of_c' = f c
 125     f_of_d' = f d
 126
 127
 128
 129 -- | Iterate the function @f@ with the initial guess @x0@ in hopes of
 130 --   finding a fixed point.
 131 fixed_point_iterations :: (a -> a) -- ^ The function @f@ to iterate.
 132                        -> a       -- ^ The initial value @x0@.
 133                        -> [a]     -- ^ The resulting sequence of x_{n}.
 134 fixed_point_iterations =
 135   iterate
 136
 137
 138 -- | Find a fixed point of the function @f@ with the search starting
 139 --   at x0. This will find the first element in the chain f(x0),
 140 --   f(f(x0)),... such that the magnitude of the difference between it
 141 --   and the next element is less than epsilon.
 142 --
 143 --   We also return the number of iterations required.
 144 --
 145 fixed_point_with_iterations :: (Normed a,
 146                                 Additive.C a,
 147                                 RealField.C b,
 148                                 Algebraic.C b)
 149                             => (a -> a)  -- ^ The function @f@ to iterate.
 150                             -> b        -- ^ The tolerance, @epsilon@.
 151                             -> a        -- ^ The initial value @x0@.
 152                             -> (Int, a) -- ^ The (iterations, fixed point) pair
 153 fixed_point_with_iterations f epsilon x0 =
 154   (fst winning_pair)
 155   where
 156     xn = fixed_point_iterations f x0
 157     xn_plus_one = tail xn
 158
 159     abs_diff v w = norm (v - w)
 160
 161     -- The nth entry in this list is the absolute value of x_{n} -
 162     -- x_{n+1}.
 163     differences = zipWith abs_diff xn xn_plus_one
 164
 165     -- This produces the list [(n, xn)] so that we can determine
 166     -- the number of iterations required.
 167     numbered_xn = zip [0..] xn
 168
 169     -- A list of pairs, (xn, |x_{n} - x_{n+1}|).
 170     pairs = zip numbered_xn differences
 171
 172     -- The pair (xn, |x_{n} - x_{n+1}|) with
 173     -- |x_{n} - x_{n+1}| < epsilon. The pattern match on 'Just' is
 174     -- "safe" since the list is infinite. We'll succeed or loop
 175     -- forever.
 176     Just winning_pair = find (\(_, diff) -> diff < epsilon) pairs