1 {-# LANGUAGE RebindableSyntax #-}
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.
11 import Data.List (find)
12 import Data.Maybe (fromMaybe)
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
23 has_root :: (RealField.C a,
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)
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
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'))
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
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
52 bisect :: (RealField.C a,
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)
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
71 -- Use a 'prime' just for consistency.
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')
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
85 trisect :: (RealField.C a,
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)
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
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')
110 (a, c, f_of_a', f_of_c')
114 else trisect f a' b' epsilon (Just fa') (Just fb')
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
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 =
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.
143 -- We also return the number of iterations required.
145 fixed_point_with_iterations :: (Normed a,
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 =
156 xn = fixed_point_iterations f x0
157 xn_plus_one = tail xn
159 abs_diff v w = norm (v - w)
161 -- The nth entry in this list is the absolute value of x_{n} -
163 differences = zipWith abs_diff xn xn_plus_one
165 -- This produces the list [(n, xn)] so that we can determine
166 -- the number of iterations required.
167 numbered_xn = zip [0..] xn
169 -- A list of pairs, (xn, |x_{n} - x_{n+1}|).
170 pairs = zip numbered_xn differences
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
176 Just winning_pair = find (\(_, diff) -> diff < epsilon) pairs