diff --git a/Data/Number/BigFloat.hs b/Data/Number/BigFloat.hs --- a/Data/Number/BigFloat.hs +++ b/Data/Number/BigFloat.hs @@ -1,110 +1,110 @@ --- | A simple implementation of floating point numbers with a selectable --- precision. The number of digits in the mantissa is selected by the --- 'Epsilon' type class from the "Fixed" module. --- --- The numbers are stored in base 10. -module Data.Number.BigFloat( - BigFloat, - Epsilon, Eps1, EpsDiv10, Prec10, Prec50, PrecPlus20 - ) where - -import Numeric(showSigned) -import Data.Number.Fixed -import qualified Data.Number.FixedFunctions as F - -base :: (Num a) => a -base = 10 - --- This representation is stupid, two Integers makes more sense, --- but is more work. --- | Floating point number where the precision is determined by the type /e/. -data BigFloat e = BF (Fixed e) Integer - deriving (Eq) - -instance (Epsilon e) => Show (BigFloat e) where - showsPrec = showSigned showBF - -- Assumes base is 10 - where showBF (BF m e) = showsPrec 0 m . showString "e" . showsPrec 0 e - -instance (Epsilon e) => Num (BigFloat e) where - BF m1 e1 + BF m2 e2 = bf (m1' + m2') e - where (m1', m2') = if e == e1 then (m1, m2 / base^(e-e2)) - else (m1 / base^(e-e1), m2) - e = e1 `max` e2 - -- Do - via negate - BF m1 e1 * BF m2 e2 = bf (m1 * m2) (e1 + e2) - negate (BF m e) = BF (-m) e - abs (BF m e) = BF (abs m) e - signum (BF m _) = bf (signum m) 0 - fromInteger i = bf (fromInteger i) 0 - -instance (Epsilon e) => Real (BigFloat e) where - toRational (BF e m) = toRational e * base^^m - -instance (Epsilon e) => Ord (BigFloat e) where - compare x y = compare (toRational x) (toRational y) - -instance (Epsilon e) => Fractional (BigFloat e) where - recip (BF m e) = bf (base / m) (-(e + 1)) - -- Take care not to lose precision for small numbers - fromRational x - | x == 0 || abs x >= 1 = bf (fromRational x) 0 - | otherwise = recip $ bf (fromRational (recip x)) 0 - - --- normalizing constructor --- XXX The scaling is very inefficient -bf :: (Epsilon e) => Fixed e -> Integer -> BigFloat e -bf m e | m == 0 = BF 0 0 - | m < 0 = - bf (-m) e - | m >= base = bf (m / base) (e + 1) - | m < 1 = bf (m * base) (e - 1) - | otherwise = BF m e - -instance (Epsilon e) => RealFrac (BigFloat e) where - properFraction x@(BF m e) = - if e < 0 then (0, x) - else let (i, f) = properFraction (m * base^^e) - in (i, bf f 0) - -instance (Epsilon e) => Floating (BigFloat e) where - pi = bf pi 0 - sqrt = toFloat1 F.sqrt - exp = toFloat1 F.exp - log = toFloat1 F.log - sin = toFloat1 F.sin - cos = toFloat1 F.cos - tan = toFloat1 F.tan - asin = toFloat1 F.asin - acos = toFloat1 F.acos - atan = toFloat1 F.atan - sinh = toFloat1 F.sinh - cosh = toFloat1 F.cosh - tanh = toFloat1 F.tanh - asinh = toFloat1 F.asinh - acosh = toFloat1 F.acosh - atanh = toFloat1 F.atanh - -instance (Epsilon e) => RealFloat (BigFloat e) where - floatRadix _ = base - floatDigits (BF m _) = - floor $ logBase base $ recip $ fromRational $ precision m - floatRange _ = (minBound, maxBound) - decodeFloat x@(BF m e) = - let d = floatDigits x - in (round $ m * base^d, fromInteger e - d) - encodeFloat m e = bf (fromInteger m) (toInteger e) - exponent (BF _ e) = fromInteger e - significand (BF m _) = BF m 0 - scaleFloat n (BF m e) = BF m (e + toInteger n) - isNaN _ = False - isInfinite _ = False - isDenormalized _ = False - isNegativeZero _ = False - isIEEE _ = False - -toFloat1 :: (Epsilon e) => (Rational -> Rational -> Rational) -> - BigFloat e -> BigFloat e -toFloat1 f x@(BF m e) = - fromRational $ f (precision m * scl) (toRational m * scl) - where scl = base^^e +-- | A simple implementation of floating point numbers with a selectable +-- precision. The number of digits in the mantissa is selected by the +-- 'Epsilon' type class from the "Fixed" module. +-- +-- The numbers are stored in base 10. +module Data.Number.BigFloat( + BigFloat, + Epsilon, Eps1, EpsDiv10, Prec10, Prec50, PrecPlus20 + ) where + +import Numeric(showSigned) +import Data.Number.Fixed +import qualified Data.Number.FixedFunctions as F + +base :: (Num a) => a +base = 10 + +-- This representation is stupid, two Integers makes more sense, +-- but is more work. +-- | Floating point number where the precision is determined by the type /e/. +data BigFloat e = BF (Fixed e) Integer + deriving (Eq) + +instance (Epsilon e) => Show (BigFloat e) where + showsPrec = showSigned showBF + -- Assumes base is 10 + where showBF (BF m e) = showsPrec 0 m . showString "e" . showsPrec 0 e + +instance (Epsilon e) => Num (BigFloat e) where + BF m1 e1 + BF m2 e2 = bf (m1' + m2') e + where (m1', m2') = if e == e1 then (m1, m2 / base^(e-e2)) + else (m1 / base^(e-e1), m2) + e = e1 `max` e2 + -- Do - via negate + BF m1 e1 * BF m2 e2 = bf (m1 * m2) (e1 + e2) + negate (BF m e) = BF (-m) e + abs (BF m e) = BF (abs m) e + signum (BF m _) = bf (signum m) 0 + fromInteger i = bf (fromInteger i) 0 + +instance (Epsilon e) => Real (BigFloat e) where + toRational (BF e m) = toRational e * base^^m + +instance (Epsilon e) => Ord (BigFloat e) where + compare x y = compare (toRational x) (toRational y) + +instance (Epsilon e) => Fractional (BigFloat e) where + recip (BF m e) = bf (base / m) (-(e + 1)) + -- Take care not to lose precision for small numbers + fromRational x + | x == 0 || abs x >= 1 = bf (fromRational x) 0 + | otherwise = recip $ bf (fromRational (recip x)) 0 + + +-- normalizing constructor +-- XXX The scaling is very inefficient +bf :: (Epsilon e) => Fixed e -> Integer -> BigFloat e +bf m e | m == 0 = BF 0 0 + | m < 0 = - bf (-m) e + | m >= base = bf (m / base) (e + 1) + | m < 1 = bf (m * base) (e - 1) + | otherwise = BF m e + +instance (Epsilon e) => RealFrac (BigFloat e) where + properFraction x@(BF m e) = + if e < 0 then (0, x) + else let (i, f) = properFraction (m * base^^e) + in (i, bf f 0) + +instance (Epsilon e) => Floating (BigFloat e) where + pi = bf pi 0 + sqrt = toFloat1 F.sqrt + exp = toFloat1 F.exp + log = toFloat1 F.log + sin = toFloat1 F.sin + cos = toFloat1 F.cos + tan = toFloat1 F.tan + asin = toFloat1 F.asin + acos = toFloat1 F.acos + atan = toFloat1 F.atan + sinh = toFloat1 F.sinh + cosh = toFloat1 F.cosh + tanh = toFloat1 F.tanh + asinh = toFloat1 F.asinh + acosh = toFloat1 F.acosh + atanh = toFloat1 F.atanh + +instance (Epsilon e) => RealFloat (BigFloat e) where + floatRadix _ = base + floatDigits (BF m _) = + floor $ logBase base $ recip $ fromRational $ precision m + floatRange _ = (minBound, maxBound) + decodeFloat x@(BF m e) = + let d = floatDigits x + in (round $ m * base^d, fromInteger e - d) + encodeFloat m e = bf (fromInteger m) (toInteger e) + exponent (BF _ e) = fromInteger e + significand (BF m _) = BF m 0 + scaleFloat n (BF m e) = BF m (e + toInteger n) + isNaN _ = False + isInfinite _ = False + isDenormalized _ = False + isNegativeZero _ = False + isIEEE _ = False + +toFloat1 :: (Epsilon e) => (Rational -> Rational -> Rational) -> + BigFloat e -> BigFloat e +toFloat1 f x@(BF m e) = + fromRational $ f (precision m * scl) (toRational m * scl) + where scl = base^^e diff --git a/Data/Number/CReal.hs b/Data/Number/CReal.hs --- a/Data/Number/CReal.hs +++ b/Data/Number/CReal.hs @@ -177,7 +177,7 @@ properFraction x@(CR x') = (fromInteger n, x - fromInteger n) where n = x' 0 instance RealFloat CReal where - floatRadix _ = error "CCeal.floatRadix" + floatRadix _ = error "CReal.floatRadix" floatDigits _ = error "CReal.floatDigits" floatRange _ = error "CReal.floatRange" decodeFloat _ = error "CReal.decodeFloat" diff --git a/Data/Number/Dif.hs b/Data/Number/Dif.hs --- a/Data/Number/Dif.hs +++ b/Data/Number/Dif.hs @@ -1,183 +1,183 @@ --- | The 'Data.Number.Dif' module contains a data type, 'Dif', that allows for --- automatic forward differentiation. --- --- All the ideas are from Jerzy Karczmarczuk\'s work, --- see . --- --- A simple example, if we define --- --- > foo x = x*x --- --- then the function --- --- > foo' = deriv foo --- --- will behave as if its body was 2*x. --- -module Data.Number.Dif(Dif, val, df, mkDif, dCon, dVar, deriv, unDif) where - --- |The 'Dif' type is the type of differentiable numbers. --- It's an instance of all the usual numeric classes. --- The computed derivative of a function is is correct --- except where the function is discontinuous, at these points --- the derivative should be a Dirac pulse, but it isn\'t. --- --- The 'Dif' numbers are printed with a trailing ~~ to --- indicate that there is a \"tail\" of derivatives. -data Dif a = D !a (Dif a) | C !a - --- |The 'dCon' function turns a normal number into a 'Dif' --- number with the same value. Not that numeric literals --- do not need an explicit conversion due to the normal --- Haskell overloading of literals. -dCon :: (Num a) => a -> Dif a -dCon x = C x - --- |The 'dVar' function turns a number into a variable --- number. This is the number with with respect to which --- the derivaticve is computed. -dVar :: (Num a, Eq a) => a -> Dif a -dVar x = D x 1 - --- |The 'df' takes a 'Dif' number and returns its first --- derivative. The function can be iterated to to get --- higher derivaties. -df :: (Num a, Eq a) => Dif a -> Dif a -df (D _ x') = x' -df (C _ ) = 0 - --- |The 'val' function takes a 'Dif' number back to a normal --- number, thus forgetting about all the derivatives. -val :: Dif a -> a -val (D x _) = x -val (C x ) = x - --- |The 'mkDif' takes a value and 'Dif' value and makes --- a 'Dif' number that has the given value as its normal --- value, and the 'Dif' number as its derivatives. -mkDif :: a -> Dif a -> Dif a -mkDif = D - --- |The 'deriv' function is a simple utility to take the --- derivative of a (single argument) function. --- It is simply defined as --- --- > deriv f = val . df . f . dVar --- -deriv :: (Num a, Num b, Eq a, Eq b) => (Dif a -> Dif b) -> (a -> b) -deriv f = val . df . f . dVar - --- |Convert a 'Dif' function to an ordinary function. -unDif :: (Num a, Eq a) => (Dif a -> Dif b) -> (a -> b) -unDif f = val . f . dVar - -instance (Show a) => Show (Dif a) where - show x = show (val x) ++ "~~" - -instance (Read a) => Read (Dif a) where - readsPrec p s = [(C x, s') | (x, s') <- readsPrec p s] - -instance (Eq a) => Eq (Dif a) where - x == y = val x == val y - -instance (Ord a) => Ord (Dif a) where - x `compare` y = val x `compare` val y - -instance (Num a, Eq a) => Num (Dif a) where - (C x) + (C y) = C (x + y) - (C x) + (D y y') = D (x + y) y' - (D x x') + (C y) = D (x + y) x' - (D x x') + (D y y') = D (x + y) (x' + y') - - (C x) - (C y) = C (x - y) - (C x) - (D y y') = D (x - y) (-y') - (D x x') - (C y) = D (x - y) x' - (D x x') - (D y y') = D (x - y) (x' - y') - - (C 0) * _ = C 0 - _ * (C 0) = C 0 - (C x) * (C y) = C (x * y) - p@(C x) * (D y y') = D (x * y) (p * y') - (D x x') * q@(C y) = D (x * y) (x' * q) - p@(D x x') * q@(D y y') = D (x * y) (x' * q + p * y') - - negate (C x) = C (negate x) - negate (D x x') = D (negate x) (negate x') - - fromInteger i = C (fromInteger i) - - abs (C x) = C (abs x) - abs p@(D x x') = D (abs x) (signum p * x') - - -- The derivative of the signum function is (2*) the Dirac impulse, - -- but there's not really any good way to encode this. - -- We could do it by +Infinity (1/0) at 0. - signum (C x) = C (signum x) - signum (D x _) = C (signum x) - -instance (Fractional a, Eq a) => Fractional (Dif a) where - recip (C x) = C (recip x) - recip (D x x') = ip - where ip = D (recip x) (-x' * ip * ip) - fromRational r = C (fromRational r) - -lift :: (Num a, Eq a) => [a -> a] -> Dif a -> Dif a -lift (f : _) (C x) = C (f x) -lift (f : f') p@(D x x') = D (f x) (x' * lift f' p) -lift _ _ = error "lift" - -instance (Floating a, Eq a) => Floating (Dif a) where - pi = C pi - - exp (C x) = C (exp x) - exp (D x x') = r where r = D (exp x) (x' * r) - - log (C x) = C (log x) - log p@(D x x') = D (log x) (x' / p) - - sqrt (C x) = C (sqrt x) - sqrt (D x x') = r where r = D (sqrt x) (x' / (2 * r)) - - sin = lift (cycle [sin, cos, negate . sin, negate . cos]) - cos = lift (cycle [cos, negate . sin, negate . cos, sin]) - - acos (C x) = C (acos x) - acos p@(D x x') = D (acos x) (-x' / sqrt(1 - p*p)) - asin (C x) = C (asin x) - asin p@(D x x') = D (asin x) ( x' / sqrt(1 - p*p)) - atan (C x) = C (atan x) - atan p@(D x x') = D (atan x) ( x' / (p*p - 1)) - - sinh x = (exp x - exp (-x)) / 2 - cosh x = (exp x + exp (-x)) / 2 - asinh x = log (x + sqrt (x*x + 1)) - acosh x = log (x + sqrt (x*x - 1)) - atanh x = (log (1 + x) - log (1 - x)) / 2 - -instance (Real a) => Real (Dif a) where - toRational = toRational . val - -instance (RealFrac a) => RealFrac (Dif a) where - -- Second component should have an impulse derivative. - properFraction x = (i, x - fromIntegral i) where (i, _) = properFraction (val x) - truncate = truncate . val - round = round . val - ceiling = ceiling . val - floor = floor . val - --- Partial definition on purpose, more could be defined. -instance (RealFloat a) => RealFloat (Dif a) where - floatRadix = floatRadix . val - floatDigits = floatDigits . val - floatRange = floatRange . val - exponent _ = 0 - scaleFloat 0 x = x - isNaN = isNaN . val - isInfinite = isInfinite . val - isDenormalized = isDenormalized . val - isNegativeZero = isNegativeZero . val - isIEEE = isIEEE . val - -- Set these to undefined rather than omit them to avoid compiler - -- warnings. - decodeFloat = undefined - encodeFloat = undefined +-- | The 'Data.Number.Dif' module contains a data type, 'Dif', that allows for +-- automatic forward differentiation. +-- +-- All the ideas are from Jerzy Karczmarczuk\'s work, +-- see . +-- +-- A simple example, if we define +-- +-- > foo x = x*x +-- +-- then the function +-- +-- > foo' = deriv foo +-- +-- will behave as if its body was 2*x. +-- +module Data.Number.Dif(Dif, val, df, mkDif, dCon, dVar, deriv, unDif) where + +-- |The 'Dif' type is the type of differentiable numbers. +-- It's an instance of all the usual numeric classes. +-- The computed derivative of a function is is correct +-- except where the function is discontinuous, at these points +-- the derivative should be a Dirac pulse, but it isn\'t. +-- +-- The 'Dif' numbers are printed with a trailing ~~ to +-- indicate that there is a \"tail\" of derivatives. +data Dif a = D !a (Dif a) | C !a + +-- |The 'dCon' function turns a normal number into a 'Dif' +-- number with the same value. Not that numeric literals +-- do not need an explicit conversion due to the normal +-- Haskell overloading of literals. +dCon :: (Num a) => a -> Dif a +dCon x = C x + +-- |The 'dVar' function turns a number into a variable +-- number. This is the number with with respect to which +-- the derivaticve is computed. +dVar :: (Num a, Eq a) => a -> Dif a +dVar x = D x 1 + +-- |The 'df' takes a 'Dif' number and returns its first +-- derivative. The function can be iterated to to get +-- higher derivaties. +df :: (Num a, Eq a) => Dif a -> Dif a +df (D _ x') = x' +df (C _ ) = 0 + +-- |The 'val' function takes a 'Dif' number back to a normal +-- number, thus forgetting about all the derivatives. +val :: Dif a -> a +val (D x _) = x +val (C x ) = x + +-- |The 'mkDif' takes a value and 'Dif' value and makes +-- a 'Dif' number that has the given value as its normal +-- value, and the 'Dif' number as its derivatives. +mkDif :: a -> Dif a -> Dif a +mkDif = D + +-- |The 'deriv' function is a simple utility to take the +-- derivative of a (single argument) function. +-- It is simply defined as +-- +-- > deriv f = val . df . f . dVar +-- +deriv :: (Num a, Num b, Eq a, Eq b) => (Dif a -> Dif b) -> (a -> b) +deriv f = val . df . f . dVar + +-- |Convert a 'Dif' function to an ordinary function. +unDif :: (Num a, Eq a) => (Dif a -> Dif b) -> (a -> b) +unDif f = val . f . dVar + +instance (Show a) => Show (Dif a) where + show x = show (val x) ++ "~~" + +instance (Read a) => Read (Dif a) where + readsPrec p s = [(C x, s') | (x, s') <- readsPrec p s] + +instance (Eq a) => Eq (Dif a) where + x == y = val x == val y + +instance (Ord a) => Ord (Dif a) where + x `compare` y = val x `compare` val y + +instance (Num a, Eq a) => Num (Dif a) where + (C x) + (C y) = C (x + y) + (C x) + (D y y') = D (x + y) y' + (D x x') + (C y) = D (x + y) x' + (D x x') + (D y y') = D (x + y) (x' + y') + + (C x) - (C y) = C (x - y) + (C x) - (D y y') = D (x - y) (-y') + (D x x') - (C y) = D (x - y) x' + (D x x') - (D y y') = D (x - y) (x' - y') + + (C 0) * _ = C 0 + _ * (C 0) = C 0 + (C x) * (C y) = C (x * y) + p@(C x) * (D y y') = D (x * y) (p * y') + (D x x') * q@(C y) = D (x * y) (x' * q) + p@(D x x') * q@(D y y') = D (x * y) (x' * q + p * y') + + negate (C x) = C (negate x) + negate (D x x') = D (negate x) (negate x') + + fromInteger i = C (fromInteger i) + + abs (C x) = C (abs x) + abs p@(D x x') = D (abs x) (signum p * x') + + -- The derivative of the signum function is (2*) the Dirac impulse, + -- but there's not really any good way to encode this. + -- We could do it by +Infinity (1/0) at 0. + signum (C x) = C (signum x) + signum (D x _) = C (signum x) + +instance (Fractional a, Eq a) => Fractional (Dif a) where + recip (C x) = C (recip x) + recip (D x x') = ip + where ip = D (recip x) (-x' * ip * ip) + fromRational r = C (fromRational r) + +lift :: (Num a, Eq a) => [a -> a] -> Dif a -> Dif a +lift (f : _) (C x) = C (f x) +lift (f : f') p@(D x x') = D (f x) (x' * lift f' p) +lift _ _ = error "lift" + +instance (Floating a, Eq a) => Floating (Dif a) where + pi = C pi + + exp (C x) = C (exp x) + exp (D x x') = r where r = D (exp x) (x' * r) + + log (C x) = C (log x) + log p@(D x x') = D (log x) (x' / p) + + sqrt (C x) = C (sqrt x) + sqrt (D x x') = r where r = D (sqrt x) (x' / (2 * r)) + + sin = lift (cycle [sin, cos, negate . sin, negate . cos]) + cos = lift (cycle [cos, negate . sin, negate . cos, sin]) + + acos (C x) = C (acos x) + acos p@(D x x') = D (acos x) (-x' / sqrt(1 - p*p)) + asin (C x) = C (asin x) + asin p@(D x x') = D (asin x) ( x' / sqrt(1 - p*p)) + atan (C x) = C (atan x) + atan p@(D x x') = D (atan x) ( x' / (p*p - 1)) + + sinh x = (exp x - exp (-x)) / 2 + cosh x = (exp x + exp (-x)) / 2 + asinh x = log (x + sqrt (x*x + 1)) + acosh x = log (x + sqrt (x*x - 1)) + atanh x = (log (1 + x) - log (1 - x)) / 2 + +instance (Real a) => Real (Dif a) where + toRational = toRational . val + +instance (RealFrac a) => RealFrac (Dif a) where + -- Second component should have an impulse derivative. + properFraction x = (i, x - fromIntegral i) where (i, _) = properFraction (val x) + truncate = truncate . val + round = round . val + ceiling = ceiling . val + floor = floor . val + +-- Partial definition on purpose, more could be defined. +instance (RealFloat a) => RealFloat (Dif a) where + floatRadix = floatRadix . val + floatDigits = floatDigits . val + floatRange = floatRange . val + exponent _ = 0 + scaleFloat 0 x = x + isNaN = isNaN . val + isInfinite = isInfinite . val + isDenormalized = isDenormalized . val + isNegativeZero = isNegativeZero . val + isIEEE = isIEEE . val + -- Set these to undefined rather than omit them to avoid compiler + -- warnings. + decodeFloat = undefined + encodeFloat = undefined diff --git a/Data/Number/Fixed.hs b/Data/Number/Fixed.hs --- a/Data/Number/Fixed.hs +++ b/Data/Number/Fixed.hs @@ -1,158 +1,158 @@ -{-# LANGUAGE - EmptyDataDecls, - GeneralizedNewtypeDeriving, - ScopedTypeVariables, - Rank2Types #-} - --- | Numbers with a fixed number of decimals. -module Data.Number.Fixed( - Fixed, - Epsilon, Eps1, EpsDiv10, Prec10, Prec50, PrecPlus20, - convertFixed, dynamicEps, precision) where -import Numeric -import Data.Char -import Data.Ratio -import qualified Data.Number.FixedFunctions as F - --- | The 'Epsilon' class contains the types that can be used to determine the --- precision of a 'Fixed' number. -class Epsilon e where - eps :: e -> Rational - --- | An epsilon of 1, i.e., no decimals. -data Eps1 -instance Epsilon Eps1 where - eps _ = 1 - --- | A type construct that gives one more decimals than the argument. -data EpsDiv10 p -instance (Epsilon e) => Epsilon (EpsDiv10 e) where - eps e = eps (un e) / 10 - where un :: EpsDiv10 e -> e - un = undefined - --- | Ten decimals. -data Prec10 -instance Epsilon Prec10 where - eps _ = 1e-10 - --- | 50 decimals. -data Prec50 -instance Epsilon Prec50 where - eps _ = 1e-50 - --- | 500 decimals. -data Prec500 -instance Epsilon Prec500 where - eps _ = 1e-500 - --- A type that gives 20 more decimals than the argument. -data PrecPlus20 e -instance (Epsilon e) => Epsilon (PrecPlus20 e) where - eps e = 1e-20 * eps (un e) - where un :: PrecPlus20 e -> e - un = undefined - ------------ - --- The type of fixed precision numbers. The type /e/ determines the precision. -newtype Fixed e = F Rational deriving (Eq, Ord, Enum, Real, RealFrac) - --- Get the accuracy (the epsilon) of the type. -precision :: (Epsilon e) => Fixed e -> Rational -precision = getEps - -instance (Epsilon e) => Num (Fixed e) where - (+) = lift2 (+) - (-) = lift2 (-) - (*) = lift2 (*) - negate (F x) = F (negate x) - abs (F x) = F (abs x) - signum (F x) = F (signum x) - fromInteger = F . fromInteger - -instance (Epsilon e) => Fractional (Fixed e) where - (/) = lift2 (/) - fromRational x = r - where r = F $ approx x (getEps r) - -lift2 :: (Epsilon e) => (Rational -> Rational -> Rational) -> Fixed e -> Fixed e -> Fixed e -lift2 op fx@(F x) (F y) = F $ approx (x `op` y) (getEps fx) - -approx :: Rational -> Rational -> Rational -approx x eps = approxRational x (eps/2) - --- | Convert between two arbitrary fixed precision types. -convertFixed :: (Epsilon e, Epsilon f) => Fixed e -> Fixed f -convertFixed e@(F x) = f - where f = F $ if feps > eeps then approx x feps else x - feps = getEps f - eeps = getEps e - -getEps :: (Epsilon e) => Fixed e -> Rational -getEps = eps . un - where un :: Fixed e -> e - un = undefined - -instance (Epsilon e) => Show (Fixed e) where - showsPrec = showSigned showFixed - where showFixed f@(F x) = showString $ show q ++ "." ++ decimals r e - where q :: Integer - (q, r) = properFraction (x + e/2) - e = getEps f - decimals a e | e >= 1 = "" - | otherwise = intToDigit b : decimals c (10 * e) - where (b, c) = properFraction (10 * a) - -instance (Epsilon e) => Read (Fixed e) where - readsPrec _ = readSigned readFixed - where readFixed s = [ (toFixed0 (approxRational x), s') | (x, s') <- readFloat s ] - -instance (Epsilon e) => Floating (Fixed e) where - pi = toFixed0 F.pi - sqrt = toFixed1 F.sqrt - exp = toFixed1 F.exp - log = toFixed1 F.log - sin = toFixed1 F.sin - cos = toFixed1 F.cos - tan = toFixed1 F.tan - asin = toFixed1 F.asin - acos = toFixed1 F.acos - atan = toFixed1 F.atan - sinh = toFixed1 F.sinh - cosh = toFixed1 F.cosh - tanh = toFixed1 F.tanh - asinh = toFixed1 F.asinh - acosh = toFixed1 F.acosh - atanh = toFixed1 F.atanh - -toFixed0 :: (Epsilon e) => (Rational -> Rational) -> Fixed e -toFixed0 f = r - where r = F $ f $ getEps r - -toFixed1 :: (Epsilon e) => (Rational -> Rational -> Rational) -> Fixed e -> Fixed e -toFixed1 f x@(F r) = F $ f (getEps x) r - -instance (Epsilon e) => RealFloat (Fixed e) where - exponent _ = 0 - scaleFloat 0 x = x - isNaN _ = False - isInfinite _ = False - isDenormalized _ = False - isNegativeZero _ = False - isIEEE _ = False - -- Explicitly undefine these rather than omitting them; this - -- prevents a compiler warning at least. - floatRadix = undefined - floatDigits = undefined - floatRange = undefined - decodeFloat = undefined - encodeFloat = undefined - ------------ - --- The call @dynmicEps r f v@ evaluates @f v@ to a precsion of @r@. -dynamicEps :: forall a . Rational -> (forall e . Epsilon e => Fixed e -> a) -> Rational -> a -dynamicEps r f v = loop (undefined :: Eps1) - where loop :: forall x . (Epsilon x) => x -> a - loop e = if eps e <= r then f (fromRational v :: Fixed x) else loop (undefined :: EpsDiv10 x) +{-# LANGUAGE + EmptyDataDecls, + GeneralizedNewtypeDeriving, + ScopedTypeVariables, + Rank2Types #-} + +-- | Numbers with a fixed number of decimals. +module Data.Number.Fixed( + Fixed, + Epsilon, Eps1, EpsDiv10, Prec10, Prec50, PrecPlus20, + convertFixed, dynamicEps, precision) where +import Numeric +import Data.Char +import Data.Ratio +import qualified Data.Number.FixedFunctions as F + +-- | The 'Epsilon' class contains the types that can be used to determine the +-- precision of a 'Fixed' number. +class Epsilon e where + eps :: e -> Rational + +-- | An epsilon of 1, i.e., no decimals. +data Eps1 +instance Epsilon Eps1 where + eps _ = 1 + +-- | A type construct that gives one more decimals than the argument. +data EpsDiv10 p +instance (Epsilon e) => Epsilon (EpsDiv10 e) where + eps e = eps (un e) / 10 + where un :: EpsDiv10 e -> e + un = undefined + +-- | Ten decimals. +data Prec10 +instance Epsilon Prec10 where + eps _ = 1e-10 + +-- | 50 decimals. +data Prec50 +instance Epsilon Prec50 where + eps _ = 1e-50 + +-- | 500 decimals. +data Prec500 +instance Epsilon Prec500 where + eps _ = 1e-500 + +-- A type that gives 20 more decimals than the argument. +data PrecPlus20 e +instance (Epsilon e) => Epsilon (PrecPlus20 e) where + eps e = 1e-20 * eps (un e) + where un :: PrecPlus20 e -> e + un = undefined + +----------- + +-- The type of fixed precision numbers. The type /e/ determines the precision. +newtype Fixed e = F Rational deriving (Eq, Ord, Enum, Real, RealFrac) + +-- Get the accuracy (the epsilon) of the type. +precision :: (Epsilon e) => Fixed e -> Rational +precision = getEps + +instance (Epsilon e) => Num (Fixed e) where + (+) = lift2 (+) + (-) = lift2 (-) + (*) = lift2 (*) + negate (F x) = F (negate x) + abs (F x) = F (abs x) + signum (F x) = F (signum x) + fromInteger = F . fromInteger + +instance (Epsilon e) => Fractional (Fixed e) where + (/) = lift2 (/) + fromRational x = r + where r = F $ approx x (getEps r) + +lift2 :: (Epsilon e) => (Rational -> Rational -> Rational) -> Fixed e -> Fixed e -> Fixed e +lift2 op fx@(F x) (F y) = F $ approx (x `op` y) (getEps fx) + +approx :: Rational -> Rational -> Rational +approx x eps = approxRational x (eps/2) + +-- | Convert between two arbitrary fixed precision types. +convertFixed :: (Epsilon e, Epsilon f) => Fixed e -> Fixed f +convertFixed e@(F x) = f + where f = F $ if feps > eeps then approx x feps else x + feps = getEps f + eeps = getEps e + +getEps :: (Epsilon e) => Fixed e -> Rational +getEps = eps . un + where un :: Fixed e -> e + un = undefined + +instance (Epsilon e) => Show (Fixed e) where + showsPrec = showSigned showFixed + where showFixed f@(F x) = showString $ show q ++ "." ++ decimals r e + where q :: Integer + (q, r) = properFraction (x + e/2) + e = getEps f + decimals a e | e >= 1 = "" + | otherwise = intToDigit b : decimals c (10 * e) + where (b, c) = properFraction (10 * a) + +instance (Epsilon e) => Read (Fixed e) where + readsPrec _ = readSigned readFixed + where readFixed s = [ (toFixed0 (approxRational x), s') | (x, s') <- readFloat s ] + +instance (Epsilon e) => Floating (Fixed e) where + pi = toFixed0 F.pi + sqrt = toFixed1 F.sqrt + exp = toFixed1 F.exp + log = toFixed1 F.log + sin = toFixed1 F.sin + cos = toFixed1 F.cos + tan = toFixed1 F.tan + asin = toFixed1 F.asin + acos = toFixed1 F.acos + atan = toFixed1 F.atan + sinh = toFixed1 F.sinh + cosh = toFixed1 F.cosh + tanh = toFixed1 F.tanh + asinh = toFixed1 F.asinh + acosh = toFixed1 F.acosh + atanh = toFixed1 F.atanh + +toFixed0 :: (Epsilon e) => (Rational -> Rational) -> Fixed e +toFixed0 f = r + where r = F $ f $ getEps r + +toFixed1 :: (Epsilon e) => (Rational -> Rational -> Rational) -> Fixed e -> Fixed e +toFixed1 f x@(F r) = F $ f (getEps x) r + +instance (Epsilon e) => RealFloat (Fixed e) where + exponent _ = 0 + scaleFloat 0 x = x + isNaN _ = False + isInfinite _ = False + isDenormalized _ = False + isNegativeZero _ = False + isIEEE _ = False + -- Explicitly undefine these rather than omitting them; this + -- prevents a compiler warning at least. + floatRadix = undefined + floatDigits = undefined + floatRange = undefined + decodeFloat = undefined + encodeFloat = undefined + +----------- + +-- The call @dynmicEps r f v@ evaluates @f v@ to a precsion of @r@. +dynamicEps :: forall a . Rational -> (forall e . Epsilon e => Fixed e -> a) -> Rational -> a +dynamicEps r f v = loop (undefined :: Eps1) + where loop :: forall x . (Epsilon x) => x -> a + loop e = if eps e <= r then f (fromRational v :: Fixed x) else loop (undefined :: EpsDiv10 x) diff --git a/Data/Number/FixedFunctions.hs b/Data/Number/FixedFunctions.hs --- a/Data/Number/FixedFunctions.hs +++ b/Data/Number/FixedFunctions.hs @@ -1,471 +1,471 @@ --- Modified by Lennart Augustsson to fit into Haskell numerical hierarchy. --- --- Module: --- --- Fraction.hs --- --- Language: --- --- Haskell --- --- Description: Rational with transcendental functionalities --- --- --- This is a generalized Rational in disguise. Rational, as a type --- synonim, could not be directly made an instance of any new class --- at all. --- But we would like it to be an instance of Transcendental, where --- trigonometry, hyperbolics, logarithms, etc. are defined. --- So here we are tiptoe-ing around, re-defining everything from --- scratch, before designing the transcendental functions -- which --- is the main motivation for this module. --- --- Aside from its ability to compute transcendentals, Fraction --- allows for denominators zero. Unlike Rational, Fraction does --- not produce run-time errors for zero denominators, but use such --- entities as indicators of invalid results -- plus or minus --- infinities. Operations on fractions never fail in principle. --- --- However, some function may compute slowly when both numerators --- and denominators of their arguments are chosen to be huge. --- For example, periodicity relations are utilized with large --- arguments in trigonometric functions to reduce the arguments --- to smaller values and thus improve on the convergence --- of continued fractions. Yet, if pi number is chosen to --- be extremely accurate then the reduced argument would --- become a fraction with huge numerator and denominator --- -- thus slowing down the entire computation of a trigonometric --- function. --- --- Usage: --- --- When computation speed is not an issue and accuracy is important --- this module replaces some of the functionalities typically handled --- by the floating point numbers: trigonometry, hyperbolics, roots --- and some special functions. All computations, including definitions --- of the basic constants pi and e, can be carried with any desired --- accuracy. One suggested usage is for mathematical servers, where --- safety might be more important than speed. See also the module --- Numerus, which supports mixed arithmetic between Integer, --- Fraction and Cofra (Complex fraction), and returns complex --- legal answers in some cases where Fraction would produce --- infinities: log (-5), sqrt (-1), etc. --- --- --- Required: --- --- Haskell Prelude --- --- Author: --- --- Jan Skibinski, Numeric Quest Inc. --- --- Date: --- --- 1998.08.16, last modified 2000.05.31 --- --- See also bottom of the page for description of the format used --- for continued fractions, references, etc. -------------------------------------------------------------------- - -module Data.Number.FixedFunctions where -import Prelude hiding (pi, sqrt, tan, atan, exp, log) -import Data.Ratio - -approx :: Rational -> Rational -> Rational -approx eps x = approxRational x eps - ------------------------------------------------------------------- --- Category: Conversion --- from continued fraction to fraction and vice versa, --- from Taylor series to continued fraction. -------------------------------------------------------------------- -type CF = [(Rational, Rational)] - -fromCF :: CF -> Rational -fromCF x = - -- - -- Convert finite continued fraction to fraction - -- evaluating from right to left. This is used - -- mainly for testing in conjunction with "toCF". - -- - foldr g 1 x - where - g :: (Rational, Rational) -> Rational -> Rational - g u v = (fst u) + (snd u) / v - -toCF :: Rational -> CF -toCF x = - -- - -- Convert fraction to finite continued fraction - -- - toCF' x [] - where - toCF' u lst = - case r of - 0 -> reverse (((q%1),(0%1)):lst) - _ -> toCF' (b%r) (((q%1),(1%1)):lst) - where - a = numerator u - b = denominator u - (q,r) = quotRem a b - - -approxCF :: Rational -> CF -> Rational -approxCF eps [] = 0 -approxCF eps x - -- - -- Approximate infinite continued fraction x by fraction, - -- evaluating from left to right, and stopping when - -- accuracy eps is achieved, or when a partial numerator - -- is zero -- as it indicates the end of CF. - -- - -- This recursive function relates continued fraction - -- to rational approximation. - -- - = approxCF' eps x 0 1 1 q' p' 1 - where - h = fst (x!!0) - (q', p') = x!!0 - approxCF' eps x v2 v1 u2 u1 a' n - | abs (1 - f1/f) < eps = approx eps f - | a == 0 = approx eps f - | otherwise = approxCF' eps x v1 v u1 u a (n+1) - where - (b, a) = x!!n - u = b*u1 + a'*u2 - v = b*v1 + a'*v2 - f = u/v - f1 = u1/v1 - - --- Type signature determined by GHC. -fromTaylorToCF :: Fractional a => [a] -> a -> [(a, a)] -fromTaylorToCF s x = - -- - -- Convert infinite number of terms of Taylor expansion of - -- a function f(x) to an infinite continued fraction, - -- where s = [s0,s1,s2,s3....] is a list of Taylor - -- series coefficients, such that f(x)=s0 + s1*x + s2*x^2.... - -- - -- Require: No Taylor coefficient is zero - -- - zero:one:[higher m | m <- [2..]] - where - zero = (s!!0, s!!1 * x) - one = (1, -s!!2/s!!1 * x) - higher m = (1 + s!!m/s!!(m-1) * x, -s!!(m+1)/s!!m * x) - - ------------------------------------------------------------------- --- Category: Auxiliaries ------------------------------------------------------------------- - -fac :: Integer -> Integer -fac = product . enumFromTo 1 - -integerRoot2 :: Integer -> Integer -integerRoot2 1 = 1 -integerRoot2 x = - -- - -- Biggest integer m, such that x - m^2 >= 0, - -- where x is a positive integer - -- - integerRoot2' 0 x (x `div` 2) x - where - integerRoot2' lo hi r y - | c > y = integerRoot2' lo r ((r + lo) `div` 2) y - | c == y = r - | otherwise = - if (r+1)^2 > y then - r - else - integerRoot2' r hi ((r + hi) `div` 2) y - where c = r^2 - -------------------------------------------------------------------- --- Everything below is the instantiation of class Transcendental --- for type Rational. See also modules Cofra and Numerus. --- --- Category: Constants -------------------------------------------------------------------- - -pi :: Rational -> Rational -pi eps = - -- - -- pi with accuracy eps - -- - -- Based on Ramanujan formula, as described in Ref. 3 - -- Accuracy: extremely good, 10^-19 for one term of continued - -- fraction - -- - (sqrt eps d) / (approxCF eps (fromTaylorToCF s x)) - where - x = 1%(640320^3)::Rational - s = [((-1)^k*(fac (6*k))%((fac k)^3*(fac (3*k))))*((a*k+b)%c) | k<-[0..]] - a = 545140134 - b = 13591409 - c = 426880 - d = 10005 - ---------------------------------------------------------------------- --- Category: Trigonometry ---------------------------------------------------------------------- - -tan :: Rational -> Rational -> Rational -tan eps 0 = 0 -tan eps x - -- - -- Tangent x computed with accuracy of eps. - -- - -- Trigonometric identities are used first to reduce - -- the value of x to a value from within the range of [-pi/2,pi/2] - -- - | x >= half_pi' = tan eps (x - ((1+m)%1)*xpi) - | x <= -half_pi' = tan eps (x + ((1+m)%1)*xpi) - --- | absx > 1 = 2 * t/(1 - t^2) - | otherwise = approxCF eps (cf x) - where - absx = abs x - t = tan eps (x/2) - m = floor ((absx - half_pi)/ xpi) - xpi = pi eps - half_pi'= 158%100 - half_pi = xpi * (1%2) - cf u = ((0%1,1%1):[((2*r + 1)/u, -1) | r <- [0..]]) - -sin :: Rational -> Rational -> Rational -sin eps 0 = 0 -sin eps x = 2*t/(1 + t*t) - where - t = tan eps (x/2) - -cos :: Rational -> Rational -> Rational -cos eps 0 = 1 -cos eps x = (1 - p)/(1 + p) - where - t = tan eps (x/2) - p = t*t - -atan :: Rational -> Rational -> Rational -atan eps x - -- - -- Inverse tangent of x with approximation eps - -- - | x == 0 = 0 - | x > 1 = (pi eps)/2 - atan eps (1/x) - | x < -1 = -(pi eps)/2 - atan eps (1/x) - | otherwise = approxCF eps ((0,x):[((2*m - 1),(m*x)^2) | m<- [1..]]) - - -asin :: Rational -> Rational -> Rational -asin eps x - -- - -- Inverse sine of x with approximation eps - -- - | x == 0 = 0 - | abs x > 1 = error "Fraction.asin" - | x == 1 = (pi eps) * (1%2) - | x == -1 = (pi eps) * (-1%2) - | otherwise = atan eps (x / (sqrt eps (1 - x^2))) - - -acos :: Rational -> Rational -> Rational -acos eps x - -- - -- Inverse cosine of x with approximation eps - -- - | x == 0 = (pi eps)*(1%2) - | abs x > 1 = error "Fraction.sin" - | x == 1 = 0 - | x == -1 = pi eps - | otherwise = atan eps ((sqrt eps (1 - x^2)) / x) - ---------------------------------------------------------------------- --- Category: Roots ---------------------------------------------------------------------- - -sqrt :: Rational -> Rational -> Rational -sqrt eps x - -- - -- Square root of x with approximation eps - -- - -- The CF pattern is: [(m,x-m^2),(2m,x-m^2),(2m,x-m^2)....] - -- where m is the biggest integer such that x-m^2 >= 0 - -- - | x < 0 = error "Fraction.sqrt" - | x == 0 = 0 - | x < 1 = 1/(sqrt eps (1/x)) - | otherwise = approxCF eps ((m,x-m^2):[(2*m,x-m^2) | r<-[0..]]) - where - m = (integerRoot2 (floor x))%1 - ---------------------------------------------------------------------- --- Category: Exponentials and hyperbolics ---------------------------------------------------------------------- - -exp :: Rational -> Rational -> Rational -exp eps x - -- - -- Exponent of x with approximation eps - -- - -- Based on Jacobi type continued fraction for exponential, - -- with fractional terms: - -- n == 0 ==> (1,x) - -- n == 1 ==> (1 -x/2, x^2/12) - -- n >= 2 ==> (1, x^2/(16*n^2 - 4)) - -- For x outside [-1,1] apply identity exp(x) = (exp(x/2))^2 - -- - | x == 0 = 1 - | x > 1 = (approxCF eps (f (x*(1%p))))^p - | x < (-1) = (approxCF eps (f (x*(1%q))))^q - | otherwise = approxCF eps (f x) - where - p = ceiling x - q = -(floor x) - f y = (1,y):(1-y/2,y^2/12):[(1,y^2/(16*n^2-4)) | n<-[2..]] - - -cosh :: Rational -> Rational -> Rational -cosh eps x = - -- - -- Hyperbolic cosine with approximation eps - -- - (a + b)*(1%2) - where - a = exp eps x - b = 1/a - -sinh :: Rational -> Rational -> Rational -sinh eps x = - -- - -- Hyperbolic sine with approximation eps - -- - (a - b)*(1%2) - where - a = exp eps x - b = 1/a - -tanh :: Rational -> Rational -> Rational -tanh eps x = - -- - -- Hyperbolic tangent with approximation eps - -- - (a - b)/ (a + b) - where - a = exp eps x - b = 1/a - -atanh :: Rational -> Rational -> Rational -atanh eps x - -- - -- Inverse hyperbolic tangent with approximation eps - -- - --- | x >= 1 = 1%0 --- | x <= -1 = -1%0 - | otherwise = (1%2) * (log eps ((1 + x) / (1 - x))) - -asinh :: Rational -> Rational -> Rational -asinh eps x - -- - -- Inverse hyperbolic sine - -- --- | x == 1%0 = 1%0 --- | x == -1%0 = -1%0 - | otherwise = log eps (x + (sqrt eps (x^2 + 1))) - -acosh :: Rational -> Rational -> Rational -acosh eps x - -- - -- Inverse hyperbolic cosine - -- --- | x == 1%0 = 1%0 --- | x < 1 = 1%0 - | otherwise = log eps (x + (sqrt eps (x^2 - 1))) - ---------------------------------------------------------------------- --- Category: Logarithms ---------------------------------------------------------------------- - -log :: Rational -> Rational -> Rational -log eps x - -- - -- Natural logarithm of strictly positive x - -- - -- Based on Stieltjes type continued fraction for log (1+y) - -- (0,y):(1,y/2):[(1,my/(4m+2)),(1,(m+1)y/(4m+2)),.... - -- (m >= 1, two elements per m) - -- Efficient only for x close to one. For larger x we recursively - -- apply the identity log(x) = log(x/2) + log(2) - -- - | x <= 0 = error "Fraction.log" - | x < 1 = -log eps (1/x) - | x == 1 = 0 - | otherwise = - case (scaled (x,0)) of - (1,s) -> (s%1) * approxCF eps (series 1) - (y,0) -> approxCF eps (series (y-1)) - (y,s) -> approxCF eps (series (y-1)) + (s%1)*approxCF eps (series 1) - where - series :: Rational -> CF - series u = (0,u):(1,u/2):[(1,u*((m+n)%(4*m + 2)))|m<-[1..],n<-[0,1]] - scaled :: (Rational,Integer) -> (Rational, Integer) - scaled (x, n) - | x == 2 = (1,n+1) - | x < 2 = (x, n) - | otherwise = scaled (x*(1%2), n+1) - - ---------------------------------------------------------------------------- --- References: --- --- 1. Classical Gosper notes on continued fraction arithmetic: --- http:%www.inwap.com/pdp10/hbaker/hakmem/cf.html --- 2. Pages on numerical constants represented as continued fractions: --- http:%www.mathsoft.com/asolve/constant/cntfrc/cntfrc.html --- 3. "Efficient on-line computation of real functions using exact floating --- point", by Peter John Potts, Imperial College --- http:%theory.doc.ic.ac.uk/~pjp/ieee.html --------------------------------------------------------------------------- - --------------------------------------------------------------------------- - --- The following representation of continued fractions is used: --- --- Continued fraction: CF representation: --- ================== ==================== --- b0 + a0 --- ------- ==> [(b0, a0), (b1, a1), (b2, a2).....] --- b1 + a1 --- ------- --- b2 + ... --- --- where "a's" and "b's" are Rationals. --- --- Many continued fractions could be represented by much simpler form --- [b1,b2,b3,b4..], where all coefficients "a" would have the same value 1 --- and would not need to be explicitely listed; and the coefficients "b" --- could be chosen as integers. --- However, there are some useful continued fractions that are --- given with fraction coefficients: "a", "b" or both. --- A fractional form can always be converted to an integer form, but --- a conversion process is not always simple and such an effort is not --- always worth of the achieved savings in the storage space or the --- computational efficiency. --- ----------------------------------------------------------------------------- --- --- Copyright: --- --- (C) 1998 Numeric Quest, All rights reserved --- --- --- --- http://www.numeric-quest.com --- --- License: --- --- GNU General Public License, GPL --- ------------------------------------------------------------------------------ +-- Modified by Lennart Augustsson to fit into Haskell numerical hierarchy. +-- +-- Module: +-- +-- Fraction.hs +-- +-- Language: +-- +-- Haskell +-- +-- Description: Rational with transcendental functionalities +-- +-- +-- This is a generalized Rational in disguise. Rational, as a type +-- synonim, could not be directly made an instance of any new class +-- at all. +-- But we would like it to be an instance of Transcendental, where +-- trigonometry, hyperbolics, logarithms, etc. are defined. +-- So here we are tiptoe-ing around, re-defining everything from +-- scratch, before designing the transcendental functions -- which +-- is the main motivation for this module. +-- +-- Aside from its ability to compute transcendentals, Fraction +-- allows for denominators zero. Unlike Rational, Fraction does +-- not produce run-time errors for zero denominators, but use such +-- entities as indicators of invalid results -- plus or minus +-- infinities. Operations on fractions never fail in principle. +-- +-- However, some function may compute slowly when both numerators +-- and denominators of their arguments are chosen to be huge. +-- For example, periodicity relations are utilized with large +-- arguments in trigonometric functions to reduce the arguments +-- to smaller values and thus improve on the convergence +-- of continued fractions. Yet, if pi number is chosen to +-- be extremely accurate then the reduced argument would +-- become a fraction with huge numerator and denominator +-- -- thus slowing down the entire computation of a trigonometric +-- function. +-- +-- Usage: +-- +-- When computation speed is not an issue and accuracy is important +-- this module replaces some of the functionalities typically handled +-- by the floating point numbers: trigonometry, hyperbolics, roots +-- and some special functions. All computations, including definitions +-- of the basic constants pi and e, can be carried with any desired +-- accuracy. One suggested usage is for mathematical servers, where +-- safety might be more important than speed. See also the module +-- Numerus, which supports mixed arithmetic between Integer, +-- Fraction and Cofra (Complex fraction), and returns complex +-- legal answers in some cases where Fraction would produce +-- infinities: log (-5), sqrt (-1), etc. +-- +-- +-- Required: +-- +-- Haskell Prelude +-- +-- Author: +-- +-- Jan Skibinski, Numeric Quest Inc. +-- +-- Date: +-- +-- 1998.08.16, last modified 2000.05.31 +-- +-- See also bottom of the page for description of the format used +-- for continued fractions, references, etc. +------------------------------------------------------------------- + +module Data.Number.FixedFunctions where +import Prelude hiding (pi, sqrt, tan, atan, exp, log) +import Data.Ratio + +approx :: Rational -> Rational -> Rational +approx eps x = approxRational x eps + +------------------------------------------------------------------ +-- Category: Conversion +-- from continued fraction to fraction and vice versa, +-- from Taylor series to continued fraction. +------------------------------------------------------------------- +type CF = [(Rational, Rational)] + +fromCF :: CF -> Rational +fromCF x = + -- + -- Convert finite continued fraction to fraction + -- evaluating from right to left. This is used + -- mainly for testing in conjunction with "toCF". + -- + foldr g 1 x + where + g :: (Rational, Rational) -> Rational -> Rational + g u v = (fst u) + (snd u) / v + +toCF :: Rational -> CF +toCF x = + -- + -- Convert fraction to finite continued fraction + -- + toCF' x [] + where + toCF' u lst = + case r of + 0 -> reverse (((q%1),(0%1)):lst) + _ -> toCF' (b%r) (((q%1),(1%1)):lst) + where + a = numerator u + b = denominator u + (q,r) = quotRem a b + + +approxCF :: Rational -> CF -> Rational +approxCF eps [] = 0 +approxCF eps x + -- + -- Approximate infinite continued fraction x by fraction, + -- evaluating from left to right, and stopping when + -- accuracy eps is achieved, or when a partial numerator + -- is zero -- as it indicates the end of CF. + -- + -- This recursive function relates continued fraction + -- to rational approximation. + -- + = approxCF' eps x 0 1 1 q' p' 1 + where + h = fst (x!!0) + (q', p') = x!!0 + approxCF' eps x v2 v1 u2 u1 a' n + | abs (1 - f1/f) < eps = approx eps f + | a == 0 = approx eps f + | otherwise = approxCF' eps x v1 v u1 u a (n+1) + where + (b, a) = x!!n + u = b*u1 + a'*u2 + v = b*v1 + a'*v2 + f = u/v + f1 = u1/v1 + + +-- Type signature determined by GHC. +fromTaylorToCF :: Fractional a => [a] -> a -> [(a, a)] +fromTaylorToCF s x = + -- + -- Convert infinite number of terms of Taylor expansion of + -- a function f(x) to an infinite continued fraction, + -- where s = [s0,s1,s2,s3....] is a list of Taylor + -- series coefficients, such that f(x)=s0 + s1*x + s2*x^2.... + -- + -- Require: No Taylor coefficient is zero + -- + zero:one:[higher m | m <- [2..]] + where + zero = (s!!0, s!!1 * x) + one = (1, -s!!2/s!!1 * x) + higher m = (1 + s!!m/s!!(m-1) * x, -s!!(m+1)/s!!m * x) + + +------------------------------------------------------------------ +-- Category: Auxiliaries +------------------------------------------------------------------ + +fac :: Integer -> Integer +fac = product . enumFromTo 1 + +integerRoot2 :: Integer -> Integer +integerRoot2 1 = 1 +integerRoot2 x = + -- + -- Biggest integer m, such that x - m^2 >= 0, + -- where x is a positive integer + -- + integerRoot2' 0 x (x `div` 2) x + where + integerRoot2' lo hi r y + | c > y = integerRoot2' lo r ((r + lo) `div` 2) y + | c == y = r + | otherwise = + if (r+1)^2 > y then + r + else + integerRoot2' r hi ((r + hi) `div` 2) y + where c = r^2 + +------------------------------------------------------------------- +-- Everything below is the instantiation of class Transcendental +-- for type Rational. See also modules Cofra and Numerus. +-- +-- Category: Constants +------------------------------------------------------------------- + +pi :: Rational -> Rational +pi eps = + -- + -- pi with accuracy eps + -- + -- Based on Ramanujan formula, as described in Ref. 3 + -- Accuracy: extremely good, 10^-19 for one term of continued + -- fraction + -- + (sqrt eps d) / (approxCF eps (fromTaylorToCF s x)) + where + x = 1%(640320^3)::Rational + s = [((-1)^k*(fac (6*k))%((fac k)^3*(fac (3*k))))*((a*k+b)%c) | k<-[0..]] + a = 545140134 + b = 13591409 + c = 426880 + d = 10005 + +--------------------------------------------------------------------- +-- Category: Trigonometry +--------------------------------------------------------------------- + +tan :: Rational -> Rational -> Rational +tan eps 0 = 0 +tan eps x + -- + -- Tangent x computed with accuracy of eps. + -- + -- Trigonometric identities are used first to reduce + -- the value of x to a value from within the range of [-pi/2,pi/2] + -- + | x >= half_pi' = tan eps (x - ((1+m)%1)*xpi) + | x <= -half_pi' = tan eps (x + ((1+m)%1)*xpi) + --- | absx > 1 = 2 * t/(1 - t^2) + | otherwise = approxCF eps (cf x) + where + absx = abs x + t = tan eps (x/2) + m = floor ((absx - half_pi)/ xpi) + xpi = pi eps + half_pi'= 158%100 + half_pi = xpi * (1%2) + cf u = ((0%1,1%1):[((2*r + 1)/u, -1) | r <- [0..]]) + +sin :: Rational -> Rational -> Rational +sin eps 0 = 0 +sin eps x = 2*t/(1 + t*t) + where + t = tan eps (x/2) + +cos :: Rational -> Rational -> Rational +cos eps 0 = 1 +cos eps x = (1 - p)/(1 + p) + where + t = tan eps (x/2) + p = t*t + +atan :: Rational -> Rational -> Rational +atan eps x + -- + -- Inverse tangent of x with approximation eps + -- + | x == 0 = 0 + | x > 1 = (pi eps)/2 - atan eps (1/x) + | x < -1 = -(pi eps)/2 - atan eps (1/x) + | otherwise = approxCF eps ((0,x):[((2*m - 1),(m*x)^2) | m<- [1..]]) + + +asin :: Rational -> Rational -> Rational +asin eps x + -- + -- Inverse sine of x with approximation eps + -- + | x == 0 = 0 + | abs x > 1 = error "Fraction.asin" + | x == 1 = (pi eps) * (1%2) + | x == -1 = (pi eps) * (-1%2) + | otherwise = atan eps (x / (sqrt eps (1 - x^2))) + + +acos :: Rational -> Rational -> Rational +acos eps x + -- + -- Inverse cosine of x with approximation eps + -- + | x == 0 = (pi eps)*(1%2) + | abs x > 1 = error "Fraction.sin" + | x == 1 = 0 + | x == -1 = pi eps + | otherwise = atan eps ((sqrt eps (1 - x^2)) / x) + +--------------------------------------------------------------------- +-- Category: Roots +--------------------------------------------------------------------- + +sqrt :: Rational -> Rational -> Rational +sqrt eps x + -- + -- Square root of x with approximation eps + -- + -- The CF pattern is: [(m,x-m^2),(2m,x-m^2),(2m,x-m^2)....] + -- where m is the biggest integer such that x-m^2 >= 0 + -- + | x < 0 = error "Fraction.sqrt" + | x == 0 = 0 + | x < 1 = 1/(sqrt eps (1/x)) + | otherwise = approxCF eps ((m,x-m^2):[(2*m,x-m^2) | r<-[0..]]) + where + m = (integerRoot2 (floor x))%1 + +--------------------------------------------------------------------- +-- Category: Exponentials and hyperbolics +--------------------------------------------------------------------- + +exp :: Rational -> Rational -> Rational +exp eps x + -- + -- Exponent of x with approximation eps + -- + -- Based on Jacobi type continued fraction for exponential, + -- with fractional terms: + -- n == 0 ==> (1,x) + -- n == 1 ==> (1 -x/2, x^2/12) + -- n >= 2 ==> (1, x^2/(16*n^2 - 4)) + -- For x outside [-1,1] apply identity exp(x) = (exp(x/2))^2 + -- + | x == 0 = 1 + | x > 1 = (approxCF eps (f (x*(1%p))))^p + | x < (-1) = (approxCF eps (f (x*(1%q))))^q + | otherwise = approxCF eps (f x) + where + p = ceiling x + q = -(floor x) + f y = (1,y):(1-y/2,y^2/12):[(1,y^2/(16*n^2-4)) | n<-[2..]] + + +cosh :: Rational -> Rational -> Rational +cosh eps x = + -- + -- Hyperbolic cosine with approximation eps + -- + (a + b)*(1%2) + where + a = exp eps x + b = 1/a + +sinh :: Rational -> Rational -> Rational +sinh eps x = + -- + -- Hyperbolic sine with approximation eps + -- + (a - b)*(1%2) + where + a = exp eps x + b = 1/a + +tanh :: Rational -> Rational -> Rational +tanh eps x = + -- + -- Hyperbolic tangent with approximation eps + -- + (a - b)/ (a + b) + where + a = exp eps x + b = 1/a + +atanh :: Rational -> Rational -> Rational +atanh eps x + -- + -- Inverse hyperbolic tangent with approximation eps + -- + +-- | x >= 1 = 1%0 +-- | x <= -1 = -1%0 + | otherwise = (1%2) * (log eps ((1 + x) / (1 - x))) + +asinh :: Rational -> Rational -> Rational +asinh eps x + -- + -- Inverse hyperbolic sine + -- +-- | x == 1%0 = 1%0 +-- | x == -1%0 = -1%0 + | otherwise = log eps (x + (sqrt eps (x^2 + 1))) + +acosh :: Rational -> Rational -> Rational +acosh eps x + -- + -- Inverse hyperbolic cosine + -- +-- | x == 1%0 = 1%0 +-- | x < 1 = 1%0 + | otherwise = log eps (x + (sqrt eps (x^2 - 1))) + +--------------------------------------------------------------------- +-- Category: Logarithms +--------------------------------------------------------------------- + +log :: Rational -> Rational -> Rational +log eps x + -- + -- Natural logarithm of strictly positive x + -- + -- Based on Stieltjes type continued fraction for log (1+y) + -- (0,y):(1,y/2):[(1,my/(4m+2)),(1,(m+1)y/(4m+2)),.... + -- (m >= 1, two elements per m) + -- Efficient only for x close to one. For larger x we recursively + -- apply the identity log(x) = log(x/2) + log(2) + -- + | x <= 0 = error "Fraction.log" + | x < 1 = -log eps (1/x) + | x == 1 = 0 + | otherwise = + case (scaled (x,0)) of + (1,s) -> (s%1) * approxCF eps (series 1) + (y,0) -> approxCF eps (series (y-1)) + (y,s) -> approxCF eps (series (y-1)) + (s%1)*approxCF eps (series 1) + where + series :: Rational -> CF + series u = (0,u):(1,u/2):[(1,u*((m+n)%(4*m + 2)))|m<-[1..],n<-[0,1]] + scaled :: (Rational,Integer) -> (Rational, Integer) + scaled (x, n) + | x == 2 = (1,n+1) + | x < 2 = (x, n) + | otherwise = scaled (x*(1%2), n+1) + + +--------------------------------------------------------------------------- +-- References: +-- +-- 1. Classical Gosper notes on continued fraction arithmetic: +-- http:%www.inwap.com/pdp10/hbaker/hakmem/cf.html +-- 2. Pages on numerical constants represented as continued fractions: +-- http:%www.mathsoft.com/asolve/constant/cntfrc/cntfrc.html +-- 3. "Efficient on-line computation of real functions using exact floating +-- point", by Peter John Potts, Imperial College +-- http:%theory.doc.ic.ac.uk/~pjp/ieee.html +-------------------------------------------------------------------------- + +-------------------------------------------------------------------------- + +-- The following representation of continued fractions is used: +-- +-- Continued fraction: CF representation: +-- ================== ==================== +-- b0 + a0 +-- ------- ==> [(b0, a0), (b1, a1), (b2, a2).....] +-- b1 + a1 +-- ------- +-- b2 + ... +-- +-- where "a's" and "b's" are Rationals. +-- +-- Many continued fractions could be represented by much simpler form +-- [b1,b2,b3,b4..], where all coefficients "a" would have the same value 1 +-- and would not need to be explicitely listed; and the coefficients "b" +-- could be chosen as integers. +-- However, there are some useful continued fractions that are +-- given with fraction coefficients: "a", "b" or both. +-- A fractional form can always be converted to an integer form, but +-- a conversion process is not always simple and such an effort is not +-- always worth of the achieved savings in the storage space or the +-- computational efficiency. +-- +---------------------------------------------------------------------------- +-- +-- Copyright: +-- +-- (C) 1998 Numeric Quest, All rights reserved +-- +-- +-- +-- http://www.numeric-quest.com +-- +-- License: +-- +-- GNU General Public License, GPL +-- +----------------------------------------------------------------------------- diff --git a/Data/Number/Interval.hs b/Data/Number/Interval.hs --- a/Data/Number/Interval.hs +++ b/Data/Number/Interval.hs @@ -1,45 +1,45 @@ --- | An incomplete implementation of interval aritrhmetic. -module Data.Number.Interval(Interval, ival, getIval) where - -data Interval a = I a a - -ival :: (Ord a) => a -> a -> Interval a -ival l h | l <= h = I l h - | otherwise = error "Interval.ival: low > high" - -getIval :: Interval a -> (a, a) -getIval (I l h) = (l, h) - -instance (Ord a) => Eq (Interval a) where - I l h == I l' h' = l == h' && h == l' - I l h /= I l' h' = h < l' || h' < l - -instance (Ord a) => Ord (Interval a) where - I l h < I l' h' = h < l' - I l h <= I l' h' = h <= l' - I l h > I l' h' = l > h' - I l h >= I l' h' = l >= h' - -- These funcions are partial, so we just leave them out. - compare _ _ = error "Interval compare" - max _ _ = error "Interval max" - min _ _ = error "Interval min" - -instance (Eq a, Show a) => Show (Interval a) where - showsPrec p (I l h) | l == h = showsPrec p l - | otherwise = showsPrec p l . showString ".." . showsPrec p h - -instance (Ord a, Num a) => Num (Interval a) where - I l h + I l' h' = I (l + l') (h + h') - I l h - I l' h' = I (l - h') (h - l') - I l h * I l' h' = I (minimum xs) (maximum xs) where xs = [l*l', l*h', h*l', h*h'] - negate (I l h) = I (-h) (-l) - -- leave out abs and signum - abs _ = error "Interval abs" - signum _ = error "Interval signum" - fromInteger i = I l l where l = fromInteger i - -instance (Ord a, Fractional a) => Fractional (Interval a) where - I l h / I l' h' | signum l' == signum h' && l' /= 0 = I (minimum xs) (maximum xs) - | otherwise = error "Interval: division by 0" - where xs = [l/l', l/h', h/l', h/h'] - fromRational r = I l l where l = fromRational r +-- | An incomplete implementation of interval aritrhmetic. +module Data.Number.Interval(Interval, ival, getIval) where + +data Interval a = I a a + +ival :: (Ord a) => a -> a -> Interval a +ival l h | l <= h = I l h + | otherwise = error "Interval.ival: low > high" + +getIval :: Interval a -> (a, a) +getIval (I l h) = (l, h) + +instance (Ord a) => Eq (Interval a) where + I l h == I l' h' = l == h' && h == l' + I l h /= I l' h' = h < l' || h' < l + +instance (Ord a) => Ord (Interval a) where + I l h < I l' h' = h < l' + I l h <= I l' h' = h <= l' + I l h > I l' h' = l > h' + I l h >= I l' h' = l >= h' + -- These funcions are partial, so we just leave them out. + compare _ _ = error "Interval compare" + max _ _ = error "Interval max" + min _ _ = error "Interval min" + +instance (Eq a, Show a) => Show (Interval a) where + showsPrec p (I l h) | l == h = showsPrec p l + | otherwise = showsPrec p l . showString ".." . showsPrec p h + +instance (Ord a, Num a) => Num (Interval a) where + I l h + I l' h' = I (l + l') (h + h') + I l h - I l' h' = I (l - h') (h - l') + I l h * I l' h' = I (minimum xs) (maximum xs) where xs = [l*l', l*h', h*l', h*h'] + negate (I l h) = I (-h) (-l) + -- leave out abs and signum + abs _ = error "Interval abs" + signum _ = error "Interval signum" + fromInteger i = I l l where l = fromInteger i + +instance (Ord a, Fractional a) => Fractional (Interval a) where + I l h / I l' h' | signum l' == signum h' && l' /= 0 = I (minimum xs) (maximum xs) + | otherwise = error "Interval: division by 0" + where xs = [l/l', l/h', h/l', h/h'] + fromRational r = I l l where l = fromRational r diff --git a/Data/Number/Natural.hs b/Data/Number/Natural.hs --- a/Data/Number/Natural.hs +++ b/Data/Number/Natural.hs @@ -1,97 +1,97 @@ --- | Lazy natural numbers. --- Addition and multiplication recurses over the first argument, i.e., --- @1 + n@ is the way to write the constant time successor function. --- --- Note that (+) and (*) are not commutative for lazy natural numbers --- when considering bottom. -module Data.Number.Natural(Natural, infinity) where - -import Data.Maybe - -data Natural = Z | S Natural - -instance Show Natural where - showsPrec p n = showsPrec p (toInteger n) - -instance Eq Natural where - x == y = x `compare` y == EQ - -instance Ord Natural where - Z `compare` Z = EQ - Z `compare` S _ = LT - S _ `compare` Z = GT - S x `compare` S y = x `compare` y - - -- (_|_) `compare` Z == _|_, but (_|_) >= Z = True - -- so for maximum laziness, we need a specialized version of (>=) and (<=) - _ >= Z = True - Z >= S _ = False - S a >= S b = a >= b - - (<=) = flip (>=) - - S x `max` S y = S (x `max` y) - x `max` y = x + y - - S x `min` S y = S (x `min` y) - _ `min` _ = Z - -maybeSubtract :: Natural -> Natural -> Maybe Natural -a `maybeSubtract` Z = Just a -S a `maybeSubtract` S b = a `maybeSubtract` b -_ `maybeSubtract` _ = Nothing - -instance Num Natural where - Z + y = y - S x + y = S (x + y) - - x - y = fromMaybe (error "Natural: (-)") (x `maybeSubtract` y) - - Z * y = Z - S x * y = y + x * y - - abs x = x - signum Z = Z - signum (S _) = S Z - - fromInteger x | x < 0 = error "Natural: fromInteger" - fromInteger 0 = Z - fromInteger x = S (fromInteger (x-1)) - -instance Integral Natural where - -- Not the most efficient version, but efficiency isn't the point of this module. :) - quotRem x y = - if x < y then - (0, x) - else - let (q, r) = quotRem (x-y) y - in (1+q, r) - div = quot - mod = rem - toInteger Z = 0 - toInteger (S x) = 1 + toInteger x - -instance Real Natural where - toRational = toRational . toInteger - -instance Enum Natural where - succ = S - pred Z = error "Natural: pred 0" - pred (S a) = a - toEnum = fromIntegral - fromEnum = fromIntegral - enumFromThenTo from thn to | from <= thn = go from (to `maybeSubtract` from) where - go from Nothing = [] - go from (Just count) = from:go (step + from) (count `maybeSubtract` step) - step = thn - from - enumFromThenTo from thn to | otherwise = go (from + step) where - go from | from >= to + step = let next = from - step in next:go next - | otherwise = [] - step = from - thn - enumFrom a = enumFromThenTo a (S a) infinity - enumFromThen a b = enumFromThenTo a b infinity - enumFromTo a c = enumFromThenTo a (S a) c - --- | The infinite natural number. -infinity :: Natural -infinity = S infinity +-- | Lazy natural numbers. +-- Addition and multiplication recurses over the first argument, i.e., +-- @1 + n@ is the way to write the constant time successor function. +-- +-- Note that (+) and (*) are not commutative for lazy natural numbers +-- when considering bottom. +module Data.Number.Natural(Natural, infinity) where + +import Data.Maybe + +data Natural = Z | S Natural + +instance Show Natural where + showsPrec p n = showsPrec p (toInteger n) + +instance Eq Natural where + x == y = x `compare` y == EQ + +instance Ord Natural where + Z `compare` Z = EQ + Z `compare` S _ = LT + S _ `compare` Z = GT + S x `compare` S y = x `compare` y + + -- (_|_) `compare` Z == _|_, but (_|_) >= Z = True + -- so for maximum laziness, we need a specialized version of (>=) and (<=) + _ >= Z = True + Z >= S _ = False + S a >= S b = a >= b + + (<=) = flip (>=) + + S x `max` S y = S (x `max` y) + x `max` y = x + y + + S x `min` S y = S (x `min` y) + _ `min` _ = Z + +maybeSubtract :: Natural -> Natural -> Maybe Natural +a `maybeSubtract` Z = Just a +S a `maybeSubtract` S b = a `maybeSubtract` b +_ `maybeSubtract` _ = Nothing + +instance Num Natural where + Z + y = y + S x + y = S (x + y) + + x - y = fromMaybe (error "Natural: (-)") (x `maybeSubtract` y) + + Z * y = Z + S x * y = y + x * y + + abs x = x + signum Z = Z + signum (S _) = S Z + + fromInteger x | x < 0 = error "Natural: fromInteger" + fromInteger 0 = Z + fromInteger x = S (fromInteger (x-1)) + +instance Integral Natural where + -- Not the most efficient version, but efficiency isn't the point of this module. :) + quotRem x y = + if x < y then + (0, x) + else + let (q, r) = quotRem (x-y) y + in (1+q, r) + div = quot + mod = rem + toInteger Z = 0 + toInteger (S x) = 1 + toInteger x + +instance Real Natural where + toRational = toRational . toInteger + +instance Enum Natural where + succ = S + pred Z = error "Natural: pred 0" + pred (S a) = a + toEnum = fromIntegral + fromEnum = fromIntegral + enumFromThenTo from thn to | from <= thn = go from (to `maybeSubtract` from) where + go from Nothing = [] + go from (Just count) = from:go (step + from) (count `maybeSubtract` step) + step = thn - from + enumFromThenTo from thn to | otherwise = go (from + step) where + go from | from >= to + step = let next = from - step in next:go next + | otherwise = [] + step = from - thn + enumFrom a = enumFromThenTo a (S a) infinity + enumFromThen a b = enumFromThenTo a b infinity + enumFromTo a c = enumFromThenTo a (S a) c + +-- | The infinite natural number. +infinity :: Natural +infinity = S infinity diff --git a/Data/Number/Symbolic.hs b/Data/Number/Symbolic.hs --- a/Data/Number/Symbolic.hs +++ b/Data/Number/Symbolic.hs @@ -1,179 +1,179 @@ --- | Symbolic number, i.e., these are not numbers at all, but just build --- a representation of the expressions. --- This implementation is incomplete in that it allows comnstruction, --- but not deconstruction of the expressions. It's mainly useful for --- debugging. -module Data.Number.Symbolic(Sym, var, con, subst, unSym) where - -import Data.Char(isAlpha) -import Data.Maybe(fromMaybe) - --- | Symbolic numbers over some base type for the literals. -data Sym a = Con a | App String ([a]->a) [Sym a] - -instance (Eq a) => Eq (Sym a) where - Con x == Con x' = x == x' - App f _ xs == App f' _ xs' = (f, xs) == (f', xs') - _ == _ = False - -instance (Ord a) => Ord (Sym a) where - Con x `compare` Con x' = x `compare` x' - Con _ `compare` App _ _ _ = LT - App _ _ _ `compare` Con _ = GT - App f _ xs `compare` App f' _ xs' = (f, xs) `compare` (f', xs') - --- | Create a variable. -var :: String -> Sym a -var s = App s undefined [] - --- | Create a constant (useful when it is not a literal). -con :: a -> Sym a -con = Con - --- | The expression @subst x v e@ substitutes the expression @v@ for each --- occurence of the variable @x@ in @e@. -subst :: (Num a, Eq a) => String -> Sym a -> Sym a -> Sym a -subst _ _ e@(Con _) = e -subst x v e@(App x' _ []) | x == x' = v - | otherwise = e -subst x v (App s f es) = - case map (subst x v) es of - [e] -> unOp (\ x -> f [x]) s e - [e1,e2] -> binOp (\ x y -> f [x,y]) e1 s e2 - es' -> App s f es' - --- Turn a symbolic number into a regular one if it is a constant, --- otherwise generate an error. -unSym :: (Show a) => Sym a -> a -unSym (Con c) = c -unSym e = error $ "unSym called: " ++ show e - -instance (Show a) => Show (Sym a) where - showsPrec p (Con c) = showsPrec p c - showsPrec _ (App s _ []) = showString s - showsPrec p (App op@(c:_) _ [x, y]) | not (isAlpha c) = - showParen (p>q) (showsPrec ql x . showString op . showsPrec qr y) - where (ql, q, qr) = fromMaybe (9,9,9) $ lookup op [ - ("**", (9,8,8)), - ("/", (7,7,8)), - ("*", (7,7,8)), - ("+", (6,6,7)), - ("-", (6,6,7))] - showsPrec p (App "negate" _ [x]) = - showParen (p>=6) (showString "-" . showsPrec 7 x) - showsPrec p (App f _ xs) = - showParen (p>10) (foldl (.) (showString f) (map (\ x -> showChar ' ' . showsPrec 11 x) xs)) - -instance (Num a, Eq a) => Num (Sym a) where - x + y = binOp (+) x "+" y - x - y = binOp (-) x "-" y - x * y = binOp (*) x "*" y - negate x = unOp negate "negate" x - abs x = unOp abs "abs" x - signum x = unOp signum "signum" x - fromInteger x = Con (fromInteger x) - -instance (Fractional a, Eq a) => Fractional (Sym a) where - x / y = binOp (/) x "/" y - fromRational x = Con (fromRational x) - --- Assume the numbers are a field and simplify a little -binOp :: (Num a, Eq a) => (a->a->a) -> Sym a -> String -> Sym a -> Sym a -binOp f (Con x) _ (Con y) = Con (f x y) -binOp _ x "+" 0 = x -binOp _ 0 "+" x = x -binOp _ x "+" (App "+" _ [y, z]) = (x + y) + z -binOp _ x "+" y | isCon y && not (isCon x) = y + x -binOp _ x "+" (App "negate" _ [y]) = x - y -binOp _ x "-" 0 = x -binOp _ x "-" x' | x == x' = 0 -binOp _ x "-" (Con y) | not (isCon x) = Con (-y) + x -binOp _ _ "*" 0 = 0 -binOp _ x "*" 1 = x -binOp _ x "*" (-1) = -x -binOp _ 0 "*" _ = 0 -binOp _ 1 "*" x = x -binOp _ (-1) "*" x = -x -binOp _ x "*" (App "*" _ [y, z]) = (x * y) * z -binOp _ x "*" y | isCon y && not (isCon x) = y * x -binOp _ x "*" (App "/" f [y, z]) = App "/" f [x*y, z] -{- -binOp _ x "*" (App "+" _ [y, z]) = x*y + x*z -binOp _ (App "+" _ [y, z]) "*" x = y*x + z*x --} -binOp _ x "/" 1 = x -binOp _ x "/" (-1) = -x -binOp _ x "/" x' | x == x' = 1 -binOp _ x "/" (App "/" f [y, z]) = App "/" f [x*z, y] -binOp f (App "**" _ [x, y]) "**" z = binOp f x "**" (y * z) -binOp _ _ "**" 0 = 1 -binOp _ 0 "**" _ = 0 -binOp f x op y = App op (\ [a,b] -> f a b) [x, y] - -unOp :: (Num a) => (a->a) -> String -> Sym a -> Sym a -unOp f _ (Con c) = Con (f c) -unOp _ "negate" (App "negate" _ [x]) = x -unOp _ "abs" e@(App "abs" _ _) = e -unOp _ "signum" e@(App "signum" _ _) = e -unOp f op x = App op (\ [a] -> f a) [x] - -isCon :: Sym a -> Bool -isCon (Con _) = True -isCon _ = False - - -instance (Integral a) => Integral (Sym a) where - quot x y = binOp quot x "quot" y - rem x y = binOp rem x "rem" y - quotRem x y = (quot x y, rem x y) - div x y = binOp div x "div" y - mod x y = binOp mod x "mod" y - toInteger (Con c) = toInteger c - -instance (Enum a) => Enum (Sym a) where - toEnum = Con . toEnum - fromEnum (Con a) = fromEnum a - -instance (Real a) => Real (Sym a) where - toRational (Con c) = toRational c - -instance (RealFrac a) => RealFrac (Sym a) where - properFraction (Con c) = (i, Con c') where (i, c') = properFraction c - -instance (Floating a, Eq a) => Floating (Sym a) where - pi = var "pi" - exp = unOp exp "exp" - sqrt = unOp sqrt "sqrt" - log = unOp log "log" - x ** y = binOp (**) x "**" y - logBase x y = binOp logBase x "logBase" y - sin = unOp sin "sin" - tan = unOp tan "tan" - cos = unOp cos "cos" - asin = unOp asin "asin" - atan = unOp atan "atan" - acos = unOp acos "acos" - sinh = unOp sinh "sinh" - tanh = unOp tanh "tanh" - cosh = unOp cosh "cosh" - asinh = unOp asinh "asinh" - atanh = unOp atanh "atanh" - acosh = unOp acosh "acosh" - -instance (RealFloat a, Show a) => RealFloat (Sym a) where - floatRadix = floatRadix . unSym - floatDigits = floatDigits . unSym - floatRange = floatRange . unSym - decodeFloat (Con c) = decodeFloat c - encodeFloat m e = Con (encodeFloat m e) - exponent (Con c) = exponent c - exponent _ = 0 - significand (Con c) = Con (significand c) - scaleFloat k (Con c) = Con (scaleFloat k c) - scaleFloat _ x = x - isNaN (Con c) = isNaN c - isInfinite (Con c) = isInfinite c - isDenormalized (Con c) = isDenormalized c - isNegativeZero (Con c) = isNegativeZero c - isIEEE = isIEEE . unSym - atan2 x y = binOp atan2 x "atan2" y +-- | Symbolic number, i.e., these are not numbers at all, but just build +-- a representation of the expressions. +-- This implementation is incomplete in that it allows comnstruction, +-- but not deconstruction of the expressions. It's mainly useful for +-- debugging. +module Data.Number.Symbolic(Sym, var, con, subst, unSym) where + +import Data.Char(isAlpha) +import Data.Maybe(fromMaybe) + +-- | Symbolic numbers over some base type for the literals. +data Sym a = Con a | App String ([a]->a) [Sym a] + +instance (Eq a) => Eq (Sym a) where + Con x == Con x' = x == x' + App f _ xs == App f' _ xs' = (f, xs) == (f', xs') + _ == _ = False + +instance (Ord a) => Ord (Sym a) where + Con x `compare` Con x' = x `compare` x' + Con _ `compare` App _ _ _ = LT + App _ _ _ `compare` Con _ = GT + App f _ xs `compare` App f' _ xs' = (f, xs) `compare` (f', xs') + +-- | Create a variable. +var :: String -> Sym a +var s = App s undefined [] + +-- | Create a constant (useful when it is not a literal). +con :: a -> Sym a +con = Con + +-- | The expression @subst x v e@ substitutes the expression @v@ for each +-- occurence of the variable @x@ in @e@. +subst :: (Num a, Eq a) => String -> Sym a -> Sym a -> Sym a +subst _ _ e@(Con _) = e +subst x v e@(App x' _ []) | x == x' = v + | otherwise = e +subst x v (App s f es) = + case map (subst x v) es of + [e] -> unOp (\ x -> f [x]) s e + [e1,e2] -> binOp (\ x y -> f [x,y]) e1 s e2 + es' -> App s f es' + +-- Turn a symbolic number into a regular one if it is a constant, +-- otherwise generate an error. +unSym :: (Show a) => Sym a -> a +unSym (Con c) = c +unSym e = error $ "unSym called: " ++ show e + +instance (Show a) => Show (Sym a) where + showsPrec p (Con c) = showsPrec p c + showsPrec _ (App s _ []) = showString s + showsPrec p (App op@(c:_) _ [x, y]) | not (isAlpha c) = + showParen (p>q) (showsPrec ql x . showString op . showsPrec qr y) + where (ql, q, qr) = fromMaybe (9,9,9) $ lookup op [ + ("**", (9,8,8)), + ("/", (7,7,8)), + ("*", (7,7,8)), + ("+", (6,6,7)), + ("-", (6,6,7))] + showsPrec p (App "negate" _ [x]) = + showParen (p>=6) (showString "-" . showsPrec 7 x) + showsPrec p (App f _ xs) = + showParen (p>10) (foldl (.) (showString f) (map (\ x -> showChar ' ' . showsPrec 11 x) xs)) + +instance (Num a, Eq a) => Num (Sym a) where + x + y = binOp (+) x "+" y + x - y = binOp (-) x "-" y + x * y = binOp (*) x "*" y + negate x = unOp negate "negate" x + abs x = unOp abs "abs" x + signum x = unOp signum "signum" x + fromInteger x = Con (fromInteger x) + +instance (Fractional a, Eq a) => Fractional (Sym a) where + x / y = binOp (/) x "/" y + fromRational x = Con (fromRational x) + +-- Assume the numbers are a field and simplify a little +binOp :: (Num a, Eq a) => (a->a->a) -> Sym a -> String -> Sym a -> Sym a +binOp f (Con x) _ (Con y) = Con (f x y) +binOp _ x "+" 0 = x +binOp _ 0 "+" x = x +binOp _ x "+" (App "+" _ [y, z]) = (x + y) + z +binOp _ x "+" y | isCon y && not (isCon x) = y + x +binOp _ x "+" (App "negate" _ [y]) = x - y +binOp _ x "-" 0 = x +binOp _ x "-" x' | x == x' = 0 +binOp _ x "-" (Con y) | not (isCon x) = Con (-y) + x +binOp _ _ "*" 0 = 0 +binOp _ x "*" 1 = x +binOp _ x "*" (-1) = -x +binOp _ 0 "*" _ = 0 +binOp _ 1 "*" x = x +binOp _ (-1) "*" x = -x +binOp _ x "*" (App "*" _ [y, z]) = (x * y) * z +binOp _ x "*" y | isCon y && not (isCon x) = y * x +binOp _ x "*" (App "/" f [y, z]) = App "/" f [x*y, z] +{- +binOp _ x "*" (App "+" _ [y, z]) = x*y + x*z +binOp _ (App "+" _ [y, z]) "*" x = y*x + z*x +-} +binOp _ x "/" 1 = x +binOp _ x "/" (-1) = -x +binOp _ x "/" x' | x == x' = 1 +binOp _ x "/" (App "/" f [y, z]) = App "/" f [x*z, y] +binOp f (App "**" _ [x, y]) "**" z = binOp f x "**" (y * z) +binOp _ _ "**" 0 = 1 +binOp _ 0 "**" _ = 0 +binOp f x op y = App op (\ [a,b] -> f a b) [x, y] + +unOp :: (Num a) => (a->a) -> String -> Sym a -> Sym a +unOp f _ (Con c) = Con (f c) +unOp _ "negate" (App "negate" _ [x]) = x +unOp _ "abs" e@(App "abs" _ _) = e +unOp _ "signum" e@(App "signum" _ _) = e +unOp f op x = App op (\ [a] -> f a) [x] + +isCon :: Sym a -> Bool +isCon (Con _) = True +isCon _ = False + + +instance (Integral a) => Integral (Sym a) where + quot x y = binOp quot x "quot" y + rem x y = binOp rem x "rem" y + quotRem x y = (quot x y, rem x y) + div x y = binOp div x "div" y + mod x y = binOp mod x "mod" y + toInteger (Con c) = toInteger c + +instance (Enum a) => Enum (Sym a) where + toEnum = Con . toEnum + fromEnum (Con a) = fromEnum a + +instance (Real a) => Real (Sym a) where + toRational (Con c) = toRational c + +instance (RealFrac a) => RealFrac (Sym a) where + properFraction (Con c) = (i, Con c') where (i, c') = properFraction c + +instance (Floating a, Eq a) => Floating (Sym a) where + pi = var "pi" + exp = unOp exp "exp" + sqrt = unOp sqrt "sqrt" + log = unOp log "log" + x ** y = binOp (**) x "**" y + logBase x y = binOp logBase x "logBase" y + sin = unOp sin "sin" + tan = unOp tan "tan" + cos = unOp cos "cos" + asin = unOp asin "asin" + atan = unOp atan "atan" + acos = unOp acos "acos" + sinh = unOp sinh "sinh" + tanh = unOp tanh "tanh" + cosh = unOp cosh "cosh" + asinh = unOp asinh "asinh" + atanh = unOp atanh "atanh" + acosh = unOp acosh "acosh" + +instance (RealFloat a, Show a) => RealFloat (Sym a) where + floatRadix = floatRadix . unSym + floatDigits = floatDigits . unSym + floatRange = floatRange . unSym + decodeFloat (Con c) = decodeFloat c + encodeFloat m e = Con (encodeFloat m e) + exponent (Con c) = exponent c + exponent _ = 0 + significand (Con c) = Con (significand c) + scaleFloat k (Con c) = Con (scaleFloat k c) + scaleFloat _ x = x + isNaN (Con c) = isNaN c + isInfinite (Con c) = isInfinite c + isDenormalized (Con c) = isDenormalized c + isNegativeZero (Con c) = isNegativeZero c + isIEEE = isIEEE . unSym + atan2 x y = binOp atan2 x "atan2" y diff --git a/Data/Number/Vectorspace.hs b/Data/Number/Vectorspace.hs --- a/Data/Number/Vectorspace.hs +++ b/Data/Number/Vectorspace.hs @@ -1,11 +1,11 @@ -{-# LANGUAGE - FunctionalDependencies, - MultiParamTypeClasses #-} -module Data.Number.Vectorspace(Vectorspace(..)) where - --- |Class of vector spaces /v/ with scalar /s/. -class Vectorspace s v | v -> s where - (*>) :: s -> v -> v - (<+>) :: v -> v -> v - vnegate :: v -> v - vzero :: v +{-# LANGUAGE + FunctionalDependencies, + MultiParamTypeClasses #-} +module Data.Number.Vectorspace(Vectorspace(..)) where + +-- |Class of vector spaces /v/ with scalar /s/. +class Vectorspace s v | v -> s where + (*>) :: s -> v -> v + (<+>) :: v -> v -> v + vnegate :: v -> v + vzero :: v diff --git a/LICENSE b/LICENSE --- a/LICENSE +++ b/LICENSE @@ -1,33 +1,33 @@ -Copyright (c) 2007-2012 -Lennart Augustsson, Russell O'Connor, Richard Smith, -Daniel Wagner, Dan Burton, Michael Orlitzky - -All rights reserved. - -Redistribution and use in source and binary forms, with or without -modification, are permitted provided that the following conditions are met: - - * Redistributions of source code must retain the above copyright - notice, this list of conditions and the following disclaimer. - - * Redistributions in binary form must reproduce the above - copyright notice, this list of conditions and the following - disclaimer in the documentation and/or other materials provided - with the distribution. - - * Neither the name of Dan Burton nor the names of other - contributors may be used to endorse or promote products derived - from this software without specific prior written permission. - -THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS -"AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT -LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR -A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT -OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, -SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT -LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, -DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY -THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT -(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE -OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. - +Copyright (c) 2007-2012 +Lennart Augustsson, Russell O'Connor, Richard Smith, +Daniel Wagner, Dan Burton, Michael Orlitzky + +All rights reserved. + +Redistribution and use in source and binary forms, with or without +modification, are permitted provided that the following conditions are met: + + * Redistributions of source code must retain the above copyright + notice, this list of conditions and the following disclaimer. + + * Redistributions in binary form must reproduce the above + copyright notice, this list of conditions and the following + disclaimer in the documentation and/or other materials provided + with the distribution. + + * Neither the name of Dan Burton nor the names of other + contributors may be used to endorse or promote products derived + from this software without specific prior written permission. + +THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS +"AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT +LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR +A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT +OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, +SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT +LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, +DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY +THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT +(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE +OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. + diff --git a/Setup.hs b/Setup.hs --- a/Setup.hs +++ b/Setup.hs @@ -1,3 +1,3 @@ -module Main where -import Distribution.Simple -main = defaultMain +module Main where +import Distribution.Simple +main = defaultMain diff --git a/Test/Data/Number/BigFloat.hs b/Test/Data/Number/BigFloat.hs --- a/Test/Data/Number/BigFloat.hs +++ b/Test/Data/Number/BigFloat.hs @@ -1,38 +1,38 @@ -module Test.Data.Number.BigFloat (bigfloat_properties) where - -import Data.Number.BigFloat (BigFloat, Prec50) - -import Test.Framework (Test, testGroup) -import Test.Framework.Providers.QuickCheck2 (testProperty) - - -prop_bigfloat_double_agree_equality :: Double -> Bool -prop_bigfloat_double_agree_equality dbl = - dbl == bf1 - where - -- Convert dbl to a BigFloat. - bf1' = realToFrac dbl :: BigFloat Prec50 - -- And convert it back. - bf1 = realToFrac bf1' :: Double - - -prop_bigfloat_double_agree_ordering :: Double -> Double -> Bool -prop_bigfloat_double_agree_ordering dbl1 dbl2 = - compare dbl1 dbl2 == compare bf1 bf2 - where - -- Convert dbl1,dbl2 to BigFloat. - bf1 = realToFrac dbl1 :: BigFloat Prec50 - bf2 = realToFrac dbl2 :: BigFloat Prec50 - - -bigfloat_properties :: Test.Framework.Test -bigfloat_properties = - testGroup "BigFloat Properties" [ - testProperty - "bigfloat/double agree (equality)" - prop_bigfloat_double_agree_equality, - - testProperty - "bigfloat/double agree (ordering)" - prop_bigfloat_double_agree_ordering - ] +module Test.Data.Number.BigFloat (bigfloat_properties) where + +import Data.Number.BigFloat (BigFloat, Prec50) + +import Test.Framework (Test, testGroup) +import Test.Framework.Providers.QuickCheck2 (testProperty) + + +prop_bigfloat_double_agree_equality :: Double -> Bool +prop_bigfloat_double_agree_equality dbl = + dbl == bf1 + where + -- Convert dbl to a BigFloat. + bf1' = realToFrac dbl :: BigFloat Prec50 + -- And convert it back. + bf1 = realToFrac bf1' :: Double + + +prop_bigfloat_double_agree_ordering :: Double -> Double -> Bool +prop_bigfloat_double_agree_ordering dbl1 dbl2 = + compare dbl1 dbl2 == compare bf1 bf2 + where + -- Convert dbl1,dbl2 to BigFloat. + bf1 = realToFrac dbl1 :: BigFloat Prec50 + bf2 = realToFrac dbl2 :: BigFloat Prec50 + + +bigfloat_properties :: Test.Framework.Test +bigfloat_properties = + testGroup "BigFloat Properties" [ + testProperty + "bigfloat/double agree (equality)" + prop_bigfloat_double_agree_equality, + + testProperty + "bigfloat/double agree (ordering)" + prop_bigfloat_double_agree_ordering + ] diff --git a/TestSuite.hs b/TestSuite.hs --- a/TestSuite.hs +++ b/TestSuite.hs @@ -1,15 +1,15 @@ -module Main -where - -import Test.Framework ( - Test, - defaultMain, - ) - -import Test.Data.Number.BigFloat (bigfloat_properties) - -main :: IO () -main = defaultMain tests - -tests :: [Test.Framework.Test] -tests = [ bigfloat_properties ] +module Main +where + +import Test.Framework ( + Test, + defaultMain, + ) + +import Test.Data.Number.BigFloat (bigfloat_properties) + +main :: IO () +main = defaultMain tests + +tests :: [Test.Framework.Test] +tests = [ bigfloat_properties ] diff --git a/numbers.cabal b/numbers.cabal --- a/numbers.cabal +++ b/numbers.cabal @@ -1,62 +1,65 @@ -Name: numbers -Version: 3000.1.0.1 -License: BSD3 -License-file: LICENSE -Author: Lennart Augustsson -Maintainer: Lennart Augustsson -Category: Data, Math -Synopsis: Various number types -Description: - Instances of the numerical classes for a variety of - different numbers: (computable) real numbers, - arbitrary precision fixed numbers, - arbitrary precision floating point numbers, - differentiable numbers, symbolic numbers, - natural numbers, interval arithmetic. -Build-type: Simple - -cabal-version: >= 1.8 - -homepage: https://github.com/DanBurton/numbers -source-repository head - type: git - location: git://github.com/DanBurton/numbers.git -source-repository this - type: git - location: git://github.com/DanBurton/numbers.git - tag: numbers-3000.0.0.0 - -Library - Build-Depends: - base >= 3 && < 5 - - Exposed-modules: - Data.Number.Symbolic Data.Number.Dif - Data.Number.CReal Data.Number.Fixed - Data.Number.Interval Data.Number.BigFloat - Data.Number.Natural - Other-modules: - Data.Number.Vectorspace - Data.Number.FixedFunctions - - Ghc-Options: - -Wall - -fno-warn-name-shadowing - -fno-warn-unused-binds - -fno-warn-unused-matches - -fno-warn-incomplete-patterns - -fno-warn-overlapping-patterns - -fno-warn-type-defaults - -test-suite testsuite - type: exitcode-stdio-1.0 - main-is: TestSuite.hs - build-depends: - base >= 3 && < 5, - -- Additional test dependencies. - QuickCheck == 2.*, - test-framework == 0.6.*, - test-framework-quickcheck2 == 0.2.* - - other-modules: - Test.Data.Number.BigFloat +Name: numbers +-- don't forget to bump the "this" source tag +Version: 3000.1.0.2 +License: BSD3 +License-file: LICENSE +Author: Lennart Augustsson +Maintainer: danburton.email@gmail.com +Category: Data, Math +Synopsis: Various number types +Description: + Instances of the numerical classes for a variety of + different numbers: (computable) real numbers, + arbitrary precision fixed numbers, + arbitrary precision floating point numbers, + differentiable numbers, symbolic numbers, + natural numbers, interval arithmetic. +Build-type: Simple + +cabal-version: >= 1.8 + +homepage: https://github.com/DanBurton/numbers#readme +bug-reports: https://github.com/DanBurton/numbers/issues + +source-repository head + type: git + location: git://github.com/DanBurton/numbers.git +source-repository this + type: git + location: git://github.com/DanBurton/numbers.git + tag: numbers-3000.1.0.2 + +Library + Build-Depends: + base >= 3 && < 5 + + Exposed-modules: + Data.Number.Symbolic Data.Number.Dif + Data.Number.CReal Data.Number.Fixed + Data.Number.Interval Data.Number.BigFloat + Data.Number.Natural + Other-modules: + Data.Number.Vectorspace + Data.Number.FixedFunctions + + Ghc-Options: + -Wall + -fno-warn-name-shadowing + -fno-warn-unused-binds + -fno-warn-unused-matches + -fno-warn-incomplete-patterns + -fno-warn-overlapping-patterns + -fno-warn-type-defaults + +test-suite testsuite + type: exitcode-stdio-1.0 + main-is: TestSuite.hs + build-depends: + base >= 3 && < 5, + -- Additional test dependencies. + QuickCheck == 2.*, + test-framework == 0.6.*, + test-framework-quickcheck2 == 0.2.* + + other-modules: + Test.Data.Number.BigFloat