numbers 3000.1.0.0 → 3000.1.0.1
raw patch · 13 files changed
+1405/−1350 lines, 13 filessetup-changedPVP ok
version bump matches the API change (PVP)
API changes (from Hackage documentation)
Files
- Data/Number/BigFloat.hs +110/−110
- Data/Number/Dif.hs +183/−183
- Data/Number/Fixed.hs +158/−158
- Data/Number/FixedFunctions.hs +471/−471
- Data/Number/Interval.hs +45/−45
- Data/Number/Natural.hs +97/−97
- Data/Number/Symbolic.hs +179/−179
- Data/Number/Vectorspace.hs +11/−11
- LICENSE +33/−33
- Setup.hs +3/−3
- Test/Data/Number/BigFloat.hs +38/−0
- TestSuite.hs +15/−0
- numbers.cabal +62/−60
Data/Number/BigFloat.hs view
@@ -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
Data/Number/Dif.hs view
@@ -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 <http://users.info.unicaen.fr/~karczma/arpap/diffalg.pdf>.------ 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 <http://users.info.unicaen.fr/~karczma/arpap/diffalg.pdf>. +-- +-- 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
Data/Number/Fixed.hs view
@@ -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)
Data/Number/FixedFunctions.hs view
@@ -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------ <jans@numeric-quest.com>------ 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 +-- +-- <jans@numeric-quest.com> +-- +-- http://www.numeric-quest.com +-- +-- License: +-- +-- GNU General Public License, GPL +-- +-----------------------------------------------------------------------------
Data/Number/Interval.hs view
@@ -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
Data/Number/Natural.hs view
@@ -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
Data/Number/Symbolic.hs view
@@ -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
Data/Number/Vectorspace.hs view
@@ -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
LICENSE view
@@ -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. +
Setup.hs view
@@ -1,3 +1,3 @@-module Main where-import Distribution.Simple-main = defaultMain+module Main where +import Distribution.Simple +main = defaultMain
+ Test/Data/Number/BigFloat.hs view
@@ -0,0 +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 + ]
+ TestSuite.hs view
@@ -0,0 +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 ]
numbers.cabal view
@@ -1,60 +1,62 @@-Name: numbers-Version: 3000.1.0.0-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- hs-source-dirs: . test- main-is: TestSuite.hs- build-depends:- base >= 3 && < 5,- -- Additional test dependencies.- QuickCheck == 2.*,- test-framework == 0.6.*,- test-framework-quickcheck2 == 0.2.*+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