Add 'diagonal' function.

[numerical-analysis.git] / src / Roots / Fast.hs

diff --git a/src/Roots/Fast.hs b/src/Roots/Fast.hs
index 0deb1fd6237a5909ea7d15a252f09b524093d1bf..f7cde8a82a70f3f98605ab87fc94fae5a59147dc 100644 (file)
--- a/src/Roots/Fast.hs
+++ b/src/Roots/Fast.hs
@@ -16,7 +16,6 @@ import NumericPrelude hiding (abs)
 import qualified Algebra.Absolute as Absolute
 import qualified Algebra.Additive as Additive
 import qualified Algebra.Algebraic as Algebraic
 import qualified Algebra.Absolute as Absolute
 import qualified Algebra.Additive as Additive
 import qualified Algebra.Algebraic as Algebraic

-import qualified Algebra.Field as Field

 import qualified Algebra.RealRing as RealRing
 import qualified Algebra.RealField as RealField
 
 import qualified Algebra.RealRing as RealRing
 import qualified Algebra.RealField as RealField
 


@@ -100,6 +99,54 @@ bisect f a b epsilon f_of_a f_of_b
 
 
 
 
 
 

+trisect :: (RealField.C a,
+            RealRing.C b,
+            Absolute.C 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
+trisect 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
+  | 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')
+    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
+
+    c = (2*a + b) / 3
+
+    d = (a + 2*b) / 3
+
+    f_of_c' = f c
+    f_of_d' = f d
+
+

 
 -- | Iterate the function @f@ with the initial guess @x0@ in hopes of
 --   finding a fixed point.
 
 -- | Iterate the function @f@ with the initial guess @x0@ in hopes of
 --   finding a fixed point.


@@ -118,7 +165,7 @@ fixed_point_iterations f x0 =
 --   We also return the number of iterations required.
 --
 fixed_point_with_iterations :: (Normed a,
 --   We also return the number of iterations required.
 --
 fixed_point_with_iterations :: (Normed a,

-                                Algebraic.C a,
+                                Additive.C a,

                                 RealField.C b,
                                 Algebraic.C b)
                             => (a -> a)  -- ^ The function @f@ to iterate.
                                 RealField.C b,
                                 Algebraic.C b)
                             => (a -> a)  -- ^ The function @f@ to iterate.