From: Michael Orlitzky Date: Wed, 20 Feb 2013 23:35:05 +0000 (-0500) Subject: Rename Aliases.hs to BigFloat.hs, now containing numeric-prelude instances for BigFloats. X-Git-Url: http://gitweb.michael.orlitzky.com/?p=numerical-analysis.git;a=commitdiff_plain;h=fe73028041fe3becce6ce1ff268181d55d54a011 Rename Aliases.hs to BigFloat.hs, now containing numeric-prelude instances for BigFloats. Rewrite Normed.hs, Roots/Fast.hs, and Roots/Simple.hs for numeric-prelude. --- diff --git a/src/Aliases.hs b/src/Aliases.hs deleted file mode 100644 index 2534e35..0000000 --- a/src/Aliases.hs +++ /dev/null @@ -1,6 +0,0 @@ -module Aliases -where - -import Data.Number.BigFloat - -type R = BigFloat Prec50 diff --git a/src/BigFloat.hs b/src/BigFloat.hs new file mode 100644 index 0000000..e23ebe6 --- /dev/null +++ b/src/BigFloat.hs @@ -0,0 +1,40 @@ +{-# LANGUAGE RebindableSyntax #-} + +module BigFloat + (module Data.Number.BigFloat) +where + +import Data.Number.BigFloat + +import NumericPrelude hiding (abs) +import qualified Algebra.Absolute as Absolute +import qualified Algebra.Additive as Additive +import qualified Algebra.Field as Field +import qualified Algebra.Ring as Ring +import qualified Algebra.ToRational as ToRational +import qualified Algebra.ZeroTestable as ZeroTestable +import qualified Prelude as P + +type R = BigFloat Prec50 + +instance Epsilon e => Additive.C (BigFloat e) where + (+) = (P.+) + zero = (P.fromInteger 0) + negate = (P.negate) + +instance Epsilon e => Ring.C (BigFloat e) where + (*) = (P.*) + fromInteger = P.fromInteger + +instance Epsilon e => Absolute.C (BigFloat e) where + abs = P.abs + signum = P.signum + +instance Epsilon e => Field.C (BigFloat e) where + recip = P.recip + +instance Epsilon e => ZeroTestable.C (BigFloat e) where + isZero = ZeroTestable.defltIsZero + +instance Epsilon e => ToRational.C (BigFloat e) where + toRational = fromRational . P.toRational diff --git a/src/Normed.hs b/src/Normed.hs index 9bef763..8554bd9 100644 --- a/src/Normed.hs +++ b/src/Normed.hs @@ -1,41 +1,53 @@ {-# LANGUAGE FlexibleInstances #-} - +{-# LANGUAGE RebindableSyntax #-} -- | The 'Normed' class represents elements of a normed vector -- space. We define instances for all common numeric types. module Normed where -import Data.Number.BigFloat +import BigFloat + +import NumericPrelude hiding (abs) +import Algebra.Absolute +import Algebra.Field +import Algebra.Ring +import Algebra.ToInteger -- Since the norm is defined on a vector space, we should be able to -- add and subtract anything on which a norm is defined. Of course -- 'Num' is a bad choice here, but we really prefer to use the normal -- addition and subtraction operators. -class (Num a) => Normed a where - norm_p :: (Integral c, RealFrac b) => c -> a -> b - norm_infty :: RealFrac b => a -> b +class (Algebra.Ring.C a, Algebra.Absolute.C a) => Normed a where + norm_p :: (Algebra.ToInteger.C c, + Algebra.Field.C b, + Algebra.Absolute.C b) + => c -> a -> b + + norm_infty :: (Algebra.Field.C b, + Algebra.Absolute.C b) + => a -> b -- | The "usual" norm. Defaults to the Euclidean norm. - norm :: RealFrac b => a -> b + norm :: (Algebra.Field.C b, Algebra.Absolute.C b) => a -> b norm = norm_p (2 :: Integer) -- Define instances for common numeric types. instance Normed Integer where - norm_p _ = fromInteger - norm_infty = fromInteger + norm_p _ = abs . fromInteger + norm_infty = abs . fromInteger instance Normed Rational where - norm_p _ = realToFrac - norm_infty = realToFrac + norm_p _ = abs . fromRational' + norm_infty = abs . fromRational' instance Epsilon e => Normed (BigFloat e) where - norm_p _ = realToFrac - norm_infty = realToFrac + norm_p _ = abs . fromRational' . toRational + norm_infty = abs . fromRational' . toRational instance Normed Float where - norm_p _ = realToFrac - norm_infty = realToFrac + norm_p _ = abs . fromRational' . toRational + norm_infty = abs . fromRational' . toRational instance Normed Double where - norm_p _ = realToFrac - norm_infty = realToFrac + norm_p _ = abs . fromRational' . toRational + norm_infty = abs . fromRational' . toRational diff --git a/src/Roots/Fast.hs b/src/Roots/Fast.hs index 5efdf3b..47fa512 100644 --- a/src/Roots/Fast.hs +++ b/src/Roots/Fast.hs @@ -1,3 +1,5 @@ +{-# LANGUAGE RebindableSyntax #-} + -- | The Roots.Fast module contains faster implementations of the -- 'Roots.Simple' algorithms. Generally, we will pass precomputed -- values to the next iteration of a function rather than passing @@ -10,8 +12,16 @@ import Data.List (find) import Normed +import NumericPrelude hiding (abs) +import Algebra.Absolute +import Algebra.Field +import Algebra.Ring -has_root :: (Fractional a, Ord a, Ord b, Num b) +has_root :: (Algebra.Field.C a, + Ord a, + Algebra.Ring.C b, + Ord b, + Algebra.Absolute.C b) => (a -> b) -- ^ The function @f@ -> a -- ^ The \"left\" endpoint, @a@ -> a -- ^ The \"right\" endpoint, @b@ @@ -51,8 +61,11 @@ has_root f a b epsilon f_of_a f_of_b = c = (a + b)/2 - -bisect :: (Fractional a, Ord a, Num b, Ord b) +bisect :: (Algebra.Field.C a, + Ord a, + Algebra.Ring.C b, + Ord b, + Algebra.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@ @@ -88,6 +101,7 @@ bisect f a b epsilon f_of_a f_of_b + -- | Iterate the function @f@ with the initial guess @x0@ in hopes of -- finding a fixed point. fixed_point_iterations :: (a -> a) -- ^ The function @f@ to iterate. @@ -104,7 +118,10 @@ fixed_point_iterations f x0 = -- -- We also return the number of iterations required. -- -fixed_point_with_iterations :: (Normed a, RealFrac b) +fixed_point_with_iterations :: (Normed a, + Algebra.Field.C b, + Algebra.Absolute.C b, + Ord b) => (a -> a) -- ^ The function @f@ to iterate. -> b -- ^ The tolerance, @epsilon@. -> a -- ^ The initial value @x0@. @@ -133,4 +150,3 @@ fixed_point_with_iterations f epsilon x0 = -- "safe" since the list is infinite. We'll succeed or loop -- forever. Just winning_pair = find (\(_, diff) -> diff < epsilon) pairs - diff --git a/src/Roots/Simple.hs b/src/Roots/Simple.hs index 6e3ff51..44d3d62 100644 --- a/src/Roots/Simple.hs +++ b/src/Roots/Simple.hs @@ -1,3 +1,5 @@ +{-# LANGUAGE RebindableSyntax #-} + -- | The Roots.Simple module contains root-finding algorithms. That -- is, procedures to (numerically) find solutions to the equation, -- @@ -15,6 +17,11 @@ import Normed import qualified Roots.Fast as F +import NumericPrelude hiding (abs) +import Algebra.Absolute +import Algebra.Field +import Algebra.Ring + -- | Does the (continuous) function @f@ have a root on the interval -- [a,b]? If f(a) <] 0 and f(b) ]> 0, we know that there's a root in -- [a,b] by the intermediate value theorem. Likewise when f(a) >= 0 @@ -34,7 +41,11 @@ import qualified Roots.Fast as F -- >>> has_root cos (-2) 2 (Just 0.001) -- True -- -has_root :: (Fractional a, Ord a, Ord b, Num b) +has_root :: (Algebra.Field.C a, + Ord a, + Algebra.Ring.C b, + Algebra.Absolute.C b, + Ord b) => (a -> b) -- ^ The function @f@ -> a -- ^ The \"left\" endpoint, @a@ -> a -- ^ The \"right\" endpoint, @b@ @@ -45,8 +56,6 @@ has_root f a b epsilon = F.has_root f a b epsilon Nothing Nothing - - -- | 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. -- @@ -66,7 +75,11 @@ has_root f a b epsilon = -- >>> bisect sin (-1) 1 0.001 -- Just 0.0 -- -bisect :: (Fractional a, Ord a, Num b, Ord b) +bisect :: (Algebra.Field.C a, + Ord a, + Algebra.Ring.C b, + Algebra.Absolute.C b, + Ord 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@ @@ -76,6 +89,65 @@ bisect f a b epsilon = F.bisect 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. +-- +fixed_point :: (Normed a, + Algebra.Field.C b, + Algebra.Absolute.C b, + Ord b) + => (a -> a) -- ^ The function @f@ to iterate. + -> b -- ^ The tolerance, @epsilon@. + -> a -- ^ The initial value @x0@. + -> a -- ^ The fixed point. +fixed_point f epsilon x0 = + snd $ F.fixed_point_with_iterations f epsilon x0 + + +-- | Return the number of iterations required to find a fixed point of +-- the function @f@ with the search starting at x0 and tolerance +-- @epsilon@. We delegate to the version that returns the number of +-- iterations and simply discard the fixed point. +fixed_point_iteration_count :: (Normed a, + Algebra.Field.C b, + Algebra.Absolute.C b, + Ord b) + => (a -> a) -- ^ The function @f@ to iterate. + -> b -- ^ The tolerance, @epsilon@. + -> a -- ^ The initial value @x0@. + -> Int -- ^ The fixed point. +fixed_point_iteration_count f epsilon x0 = + fst $ F.fixed_point_with_iterations f epsilon x0 + + +-- | Returns a list of ratios, +-- +-- ||x^{*} - x_{n+1}|| / ||x^{*} - x_{n}||^{p} +-- +-- of fixed point iterations for the function @f@ with initial guess +-- @x0@ and @p@ some positive power. +-- +-- This is used to determine the rate of convergence. +-- +fixed_point_error_ratios :: (Normed a, + Algebra.Field.C b, + Algebra.Absolute.C b, + Ord b) + => (a -> a) -- ^ The function @f@ to iterate. + -> a -- ^ The initial value @x0@. + -> a -- ^ The true solution, @x_star@. + -> Integer -- ^ The power @p@. + -> [b] -- ^ The resulting sequence of x_{n}. +fixed_point_error_ratios f x0 x_star p = + zipWith (/) en_plus_one en_exp + where + xn = F.fixed_point_iterations f x0 + en = map (\x -> norm (x_star - x)) xn + en_plus_one = tail en + en_exp = map (^p) en + + -- | The sequence x_{n} of values obtained by applying Newton's method -- on the function @f@ and initial guess @x0@. @@ -88,7 +160,7 @@ bisect f a b epsilon = -- >>> tail $ take 4 $ newton_iterations f f' 2 -- [1.6806282722513088,1.4307389882390624,1.2549709561094362] -- -newton_iterations :: (Fractional a, Ord a) +newton_iterations :: (Algebra.Field.C a) => (a -> a) -- ^ The function @f@ whose root we seek -> (a -> a) -- ^ The derivative of @f@ -> a -- ^ Initial guess, x-naught @@ -100,7 +172,6 @@ newton_iterations f f' x0 = xn - ( (f xn) / (f' xn) ) - -- | Use Newton's method to find a root of @f@ near the initial guess -- @x0@. If your guess is bad, this will recurse forever! -- @@ -124,7 +195,7 @@ newton_iterations f f' x0 = -- >>> abs (f root) < eps -- True -- -newtons_method :: (Fractional a, Ord a) +newtons_method :: (Algebra.Field.C a, Algebra.Absolute.C a, Ord a) => (a -> a) -- ^ The function @f@ whose root we seek -> (a -> a) -- ^ The derivative of @f@ -> a -- ^ The tolerance epsilon @@ -136,7 +207,6 @@ newtons_method f f' epsilon x0 = x_n = newton_iterations f f' x0 - -- | Takes a function @f@ of two arguments and repeatedly applies @f@ -- to the previous two values. Returns a list containing all -- generated values, f(x0, x1), f(x1, x2), f(x2, x3)... @@ -158,6 +228,7 @@ iterate2 f x0 x1 = let next = f prev2 prev1 in next : go prev1 next + -- | The sequence x_{n} of values obtained by applying the secant -- method on the function @f@ and initial guesses @x0@, @x1@. -- @@ -174,7 +245,7 @@ iterate2 f x0 x1 = -- >>> take 4 $ secant_iterations f 2 1 -- [2.0,1.0,1.0161290322580645,1.190577768676638] -- -secant_iterations :: (Fractional a, Ord a) +secant_iterations :: (Algebra.Field.C a) => (a -> a) -- ^ The function @f@ whose root we seek -> a -- ^ Initial guess, x-naught -> a -- ^ Second initial guess, x-one @@ -202,7 +273,7 @@ secant_iterations f x0 x1 = -- >>> abs (f root) < (1/10^9) -- True -- -secant_method :: (Fractional a, Ord a) +secant_method :: (Algebra.Field.C a, Algebra.Absolute.C a, Ord a) => (a -> a) -- ^ The function @f@ whose root we seek -> a -- ^ The tolerance epsilon -> a -- ^ Initial guess, x-naught @@ -212,54 +283,3 @@ secant_method f epsilon x0 x1 = find (\x -> abs (f x) < epsilon) x_n where x_n = secant_iterations f x0 x1 - - - --- | 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. --- -fixed_point :: (Normed a, RealFrac b) - => (a -> a) -- ^ The function @f@ to iterate. - -> b -- ^ The tolerance, @epsilon@. - -> a -- ^ The initial value @x0@. - -> a -- ^ The fixed point. -fixed_point f epsilon x0 = - snd $ F.fixed_point_with_iterations f epsilon x0 - - --- | Return the number of iterations required to find a fixed point of --- the function @f@ with the search starting at x0 and tolerance --- @epsilon@. We delegate to the version that returns the number of --- iterations and simply discard the fixed point. -fixed_point_iteration_count :: (Normed a, RealFrac b) - => (a -> a) -- ^ The function @f@ to iterate. - -> b -- ^ The tolerance, @epsilon@. - -> a -- ^ The initial value @x0@. - -> Int -- ^ The fixed point. -fixed_point_iteration_count f epsilon x0 = - fst $ F.fixed_point_with_iterations f epsilon x0 - - --- | Returns a list of ratios, --- --- ||x^{*} - x_{n+1}|| / ||x^{*} - x_{n}||^{p} --- --- of fixed point iterations for the function @f@ with initial guess --- @x0@ and @p@ some positive power. --- --- This is used to determine the rate of convergence. --- -fixed_point_error_ratios :: (Normed a, RealFrac b) - => (a -> a) -- ^ The function @f@ to iterate. - -> a -- ^ The initial value @x0@. - -> a -- ^ The true solution, @x_star@. - -> Integer -- ^ The power @p@. - -> [b] -- ^ The resulting sequence of x_{n}. -fixed_point_error_ratios f x0 x_star p = - zipWith (/) en_plus_one en_exp - where - xn = F.fixed_point_iterations f x0 - en = map (\x -> norm (x_star - x)) xn - en_plus_one = tail en - en_exp = map (^p) en