From 04f56a8882bb0c574b603f8c3fed9481ea934f7f Mon Sep 17 00:00:00 2001 From: Michael Orlitzky Date: Sat, 8 Jun 2013 17:19:36 -0400 Subject: [PATCH] Update numeric-prelude and fixed-vector. Use the (!) function from the new fixed-vector. Add the Fast/Simple trisection method and a test. --- numerical-analysis.cabal | 4 ++-- src/Linear/Matrix.hs | 3 ++- src/Linear/Vector.hs | 14 ++---------- src/Roots/Fast.hs | 48 ++++++++++++++++++++++++++++++++++++++++ src/Roots/Simple.hs | 44 +++++++++++++++++++++++++++++++++--- 5 files changed, 95 insertions(+), 18 deletions(-) diff --git a/numerical-analysis.cabal b/numerical-analysis.cabal index e459791..0c2c30e 100644 --- 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 dbd321f..b17d0eb 100644 --- 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 aa1568b..42d47a3 100644 --- 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 b3b5818..f7cde8a 100644 --- 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 414b777..a6aa09e 100644 --- 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. -- 2.43.2