Update numeric-prelude and fixed-vector.

Use the (!) function from the new fixed-vector.
Add the Fast/Simple trisection method and a test.

diff --git a/numerical-analysis.cabal b/numerical-analysis.cabal
index e459791a107b894a9475d99e0e4477b8e48959ae..0c2c30e6d3d5d96db0125a6c9e4ea88651fc22c2 100644 (file)
--- a/numerical-analysis.cabal
+++ b/numerical-analysis.cabal
@@ -23,9 +23,9 @@ library
 
   build-depends:
     base              >= 3 && < 5,
-    fixed-vector      == 0.2.*,
+    fixed-vector      == 0.4.*,
     numbers           == 3000.1.*,
-    numeric-prelude   == 0.3.*
+    numeric-prelude   == 0.4.*
 
   hs-source-dirs:
     src/

diff --git a/src/Linear/Matrix.hs b/src/Linear/Matrix.hs
index dbd321fefc6e076ade9344495ddfd1c6bd26858c..b17d0eb391484767cfed99d70f9b6ee7268a2394 100644 (file)
--- a/src/Linear/Matrix.hs
+++ b/src/Linear/Matrix.hs
@@ -17,6 +17,7 @@ where
 import Data.List (intercalate)
 
 import Data.Vector.Fixed (
+  (!),
   N1,
   N2,
   N3,
@@ -42,7 +43,7 @@ import qualified Data.Vector.Fixed as V (
   zipWith
   )
 import Data.Vector.Fixed.Boxed (Vec)
-import Data.Vector.Fixed.Internal (Arity, arity)
+import Data.Vector.Fixed.Internal.Arity (Arity, arity)
 import Linear.Vector
 import Normed
 

diff --git a/src/Linear/Vector.hs b/src/Linear/Vector.hs
index aa1568bec6c2ff3b8ed0baa724177349fb04e3bb..42d47a39e58670136b5300dd1f83f007ce85fc2e 100644 (file)
--- a/src/Linear/Vector.hs
+++ b/src/Linear/Vector.hs
@@ -18,6 +18,7 @@ import Data.Vector.Fixed (
   toList,
   )
 import qualified Data.Vector.Fixed as V (
+  (!),
   length,
   )
 import Data.Vector.Fixed.Boxed
@@ -27,17 +28,6 @@ type Vec4 = Vec N4
 type Vec5 = Vec N5
 
 
--- | Unsafe indexing.
---
---   Examples:
---
---   >>> import Data.Vector.Fixed (mk2)
---   >>> let v1 = mk2 1 2 :: Vec2 Int
---   >>> v1 ! 1
---   2
---
-(!) :: (Vector v a) => v a -> Int -> a
-(!) v1 idx = (toList v1) !! idx
 
 -- | Safe indexing.
 --
@@ -53,7 +43,7 @@ type Vec5 = Vec N5
 (!?) :: (Vector v a) => v a -> Int -> Maybe a
 (!?) v1 idx
   | idx < 0 || idx >= V.length v1 = Nothing
-  | otherwise                     = Just $ v1 ! idx
+  | otherwise                     = Just $ v1 V.! idx
 
 
 -- | Remove an element of the given vector.

diff --git a/src/Roots/Fast.hs b/src/Roots/Fast.hs
index b3b5818782fede0db3e6b465dc9186d8e29a1ca5..f7cde8a82a70f3f98605ab87fc94fae5a59147dc 100644 (file)
--- a/src/Roots/Fast.hs
+++ b/src/Roots/Fast.hs
@@ -99,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.

diff --git a/src/Roots/Simple.hs b/src/Roots/Simple.hs
index 414b77787262e608ceef4a66981a5dc5ab82cc50..a6aa09e497ba841d2c3a8e0be2749aab27c72662 100644 (file)
--- a/src/Roots/Simple.hs
+++ b/src/Roots/Simple.hs
@@ -56,7 +56,7 @@ has_root f a b epsilon =
 
 
 -- | We are given a function @f@ and an interval [a,b]. The bisection
---   method checks finds a root by splitting [a,b] in half repeatedly.
+--   method finds a root by splitting [a,b] in half repeatedly.
 --
 --   If one is found within some prescribed tolerance @epsilon@, it is
 --   returned. Otherwise, the interval [a,b] is split into two
@@ -68,8 +68,12 @@ has_root f a b epsilon =
 --
 --   Examples:
 --
---   >>> bisect cos 1 2 0.001
---   Just 1.5712890625
+--   >>> let actual = 1.5707963267948966
+--   >>> let Just root = bisect cos 1 2 0.001
+--   >>> root
+--   1.5712890625
+--   >>> abs (root - actual) < 0.001
+--   True
 --
 --   >>> bisect sin (-1) 1 0.001
 --   Just 0.0
@@ -84,6 +88,40 @@ bisect f a b epsilon =
   F.bisect f a b epsilon Nothing Nothing
 
 
+-- | We are given a function @f@ and an interval [a,b]. The trisection
+--   method finds a root by splitting [a,b] into three
+--   subintervals repeatedly.
+--
+--   If one is found within some prescribed tolerance @epsilon@, it is
+--   returned. Otherwise, the interval [a,b] is split into two
+--   subintervals [a,c] and [c,b] of equal length which are then both
+--   checked via the same process.
+--
+--   Returns 'Just' the value x for which f(x) == 0 if one is found,
+--   or Nothing if one of the preconditions is violated.
+--
+--   Examples:
+--
+--   >>> let actual = 1.5707963267948966
+--   >>> let Just root = trisect cos 1 2 0.001
+--   >>> root
+--   1.5713305898491083
+--   >>> abs (root - actual) < 0.001
+--   True
+--
+--   >>> trisect sin (-1) 1 0.001
+--   Just 0.0
+--
+trisect :: (RealField.C a, RealRing.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 a
+trisect f a b epsilon =
+  F.trisect f a b epsilon Nothing Nothing
+
+
 -- | Find a fixed point of the function @f@ with the search starting
 --   at x0. We delegate to the version that returns the number of
 --   iterations and simply discard the number of iterations.