numbers (empty) → 2007.4.29
raw patch · 9 files changed
+1255/−0 lines, 9 filesdep +basebuild-type:Customsetup-changed
Dependencies added: base
Files
- Data/Number/CReal.hs +240/−0
- Data/Number/Dif.hs +179/−0
- Data/Number/Fixed.hs +129/−0
- Data/Number/FixedFunctions.hs +469/−0
- Data/Number/Interval.hs +44/−0
- Data/Number/Symbolic.hs +167/−0
- Data/Number/Vectorspace.hs +9/−0
- Setup.hs +3/−0
- numbers.cabal +15/−0
+ Data/Number/CReal.hs view
@@ -0,0 +1,240 @@+{-# OPTIONS -fglasgow-exts #-} +-- ERA: Exact Real Arithmetic (version 1.0) +-- +-- A tolerably efficient and possibly correct implementation of the computable +-- reals using Haskell 1.2. +-- +-- David Lester, Department of Computer Science, Manchester University, M13 9PL. +-- (2000-2001) + +-- #hide +module Data.Number.CReal(CReal, showCReal) where +import Data.Ratio +import Numeric(readFloat, readSigned) + +-- |The 'CReal' type implements (constructive) real numbers. +-- +-- Note that the comparison operations on 'CReal' may diverge +-- since it is (by necessity) impossible to implementent them +-- correctly and always terminating. +-- +-- This implementation is really David Lester's ERA package. +data CReal = CR (Int -> Integer) + +instance Eq CReal where + x == y = s' (digitsToBits digits) == 0 where (CR s') = x-y + +instance Ord CReal where + x <= y = s' (digitsToBits digits) <= 0 where (CR s') = x-y + x < y = s' (digitsToBits digits) < 0 where (CR s') = x-y + x >= y = s' (digitsToBits digits) >= 0 where (CR s') = x-y + x > y = s' (digitsToBits digits) > 0 where (CR s') = x-y + max (CR x') (CR y') = CR (\p -> max (x' p) (y' p)) + min (CR x') (CR y') = CR (\p -> min (x' p) (y' p)) + +instance Num CReal where + (CR x') + (CR y') = CR (\p -> round_uk ((x' (p+2) + y' (p+2)) % 4)) + (CR x') * (CR y') = CR (\p -> round_uk ((x' (p+sy)*y' (p+sx)) % 2^(p+sx+sy))) + where x0 = abs (x' 0)+2; y0 = abs (y' 0)+2 + sx = sizeinbase x0 2+3; sy = sizeinbase y0 2+3 + negate (CR x') = CR (\p -> negate (x' p)) + abs x = max x (negate x) + signum (CR x') = fromInteger (signum (x' (digitsToBits digits))) + fromInteger n = CR (\p -> n*2^p) + +instance Fractional CReal where + recip (CR x') = CR (\p -> let s = head [n | n <- [0..], 3 <= abs (x' n)] + in round_uk (2^(2*p+2*s+2) % (x' (p+2*s+2)))) + fromRational x = fromInteger (numerator x) / fromInteger (denominator x) + +-- two useful scaling functions: + +div2n :: CReal -> Int -> CReal +div2n (CR x') n = CR (\p -> if p >= n then x' (p-n) else round_uk (x' p % 2^n)) + +mul2n :: CReal -> Int -> CReal +mul2n (CR x') n = CR (\p -> x' (p+n)) + +-- transcendental functions (mostly range reductions): + +instance Floating CReal where + pi = 16 * atan (fromRational (1 % 5)) + - 4 * atan (fromRational (1 % 239)) + sqrt x = CR (\p -> floorsqrt (x' (2*p))) where (CR x') = x + + log x = if t < 0 then error "log of negative number\n" else + if t < 4 then - log (recip x) else + if t < 8 then log_dr x else + {- 7 < t -} log_dr (div2n x n) + fromIntegral n * log2 + where (CR x') = x; t = x' 2; n = sizeinbase t 2 - 3 + exp x = if n < 0 then div2n (exp_dr s) (fromInteger (-n)) else + if n > 0 then mul2n (exp_dr s) (fromInteger n) else exp_dr s + where (CR u') = x/log2; n = u' 0; s = x-fromInteger n*log2 + sin x = if n == 0 then sin_dr y else + if n == 1 then sqrt1By2 * (cos_dr y + sin_dr y) else + if n == 2 then cos_dr y else + if n == 3 then sqrt1By2 * (cos_dr y - sin_dr y) else + if n == 4 then - sin_dr y else + if n == 5 then - sqrt1By2 * (cos_dr y + sin_dr y) else + if n == 6 then - cos_dr y else + {- n == 7 -} - sqrt1By2 * (cos_dr y - sin_dr y) + where (CR z') = x/piBy4; s = round_uk (z' 2 % 4); n = s `mod` 8 + y = x - piBy4 * fromInteger s + cos x = if n == 0 then cos_dr y else + if n == 1 then sqrt1By2 * (cos_dr y - sin_dr y) else + if n == 2 then sin_dr y else + if n == 3 then sqrt1By2 * (cos_dr y + sin_dr y) else + if n == 4 then - cos_dr y else + if n == 5 then - sqrt1By2 * (cos_dr y - sin_dr y) else + if n == 6 then - sin_dr y else + {- n == 7 -} - sqrt1By2 * (cos_dr y + sin_dr y) + where (CR z') = x/piBy4; s = round_uk (z' 2 % 4); n = s `mod` 8 + y = x - piBy4 * fromInteger s + atan x = if t < -5 then atan_dr (negate (recip x)) - piBy2 else + if t == -4 then -piBy4 - atan_dr (xp1/xm1) else + if t < 4 then atan_dr x else + if t == 4 then piBy4 + atan_dr (xm1/xp1) else + {- t > 4 -} piBy2 - atan_dr (recip x) + where (CR x') = x; t = x' 2 + xp1 = x+1; xm1 = x-1 + asin x = if x0 > 0 then pi / 2 - atan (s/x) else + if x0 == 0 then atan (x/s) else + {- x0 < 0 -} atan (s/x) - pi / 2 + where (CR x') = x; x0 = x' 0; s = sqrt (1 - x*x) + acos x = pi / 2 - asin x + sinh x = (y - recip y) / 2 where y = exp x + cosh x = (y + recip y) / 2 where y = exp x + tanh x = (y - y') / (y + y') where y = exp x; y' = recip y + asinh x = log (x + sqrt (x*x + 1)) + acosh x = log (x + sqrt (x*x - 1)) + atanh x = log ((1 + x) / (1 - x)) / 2 + + +acc_seq :: (Rational -> Integer -> Rational) -> [Rational] +acc_seq f = scanl f (1 % 1) [1..] + +exp_dr :: CReal -> CReal +exp_dr = power_series (acc_seq (\a n -> a*(1 % n))) id + +log_dr :: CReal -> CReal +log_dr x = y * log_drx y where y = (x - 1) / x + +log_drx :: CReal -> CReal +log_drx = power_series [1 % n | n <- [1..]] (+1) + +sin_dr :: CReal -> CReal +sin_dr x = x*power_series (acc_seq (\a n -> -a*(1 % (2*n*(2*n+1))))) id (x*x) + +cos_dr :: CReal -> CReal +cos_dr x = power_series (acc_seq (\a n -> -a*(1 % (2*n*(2*n-1))))) id (x*x) + +atan_dr :: CReal -> CReal +atan_dr x = (x/y) * atan_drx ((x*x)/y) where y = x*x+1 + +atan_drx :: CReal -> CReal +atan_drx = power_series (acc_seq (\a n -> a*((2*n) % (2*n+1)))) (+1) + +-- power_series takes as arguments: +-- a (rational) list of the coefficients of the power series +-- a function from the desired accuracy to the number of terms needed +-- the argument x + +power_series :: [Rational] -> (Int -> Int) -> CReal -> CReal +power_series ps terms (CR x') + = CR (\p -> let t = terms p; l2t = 2*sizeinbase (toInteger t+1) 2+6; p' = p + l2t + xr = x' p'; xn = 2^p'; g yn = round_uk ((yn*xr) % (2^p')) + in round_uk (accumulate (iterate g xn) (take t ps) % (2^l2t))) + where accumulate _ [] = 0 + accumulate [] _ = error "CReal.power_series.accumulate" + accumulate (x:xs) (c:cs) = let t = round_uk (c*(x % 1)) in + if t == 0 then 0 else t + accumulate xs cs + +-- Some useful constants: + +piBy2 :: CReal +piBy2 = div2n pi 1 + +piBy4 :: CReal +piBy4 = div2n pi 2 + +log2 :: CReal +log2 = div2n (log_drx (recip 2)) 1 + +sqrt1By2 :: CReal +sqrt1By2 = sqrt (recip 2) + +instance Enum CReal where + toEnum i = fromIntegral i + fromEnum _ = error "Cannot fromEnum CReal" + enumFrom = iterate (+ 1) + enumFromTo n e = takeWhile (<= e) $ iterate (+ 1)n + enumFromThen n m = iterate (+(m-n)) n + enumFromThenTo n m e = if m >= n then takeWhile (<= e) $ iterate (+(m-n)) n + else takeWhile (>= e) $ iterate (+(m-n)) n + +instance Real CReal where + -- toRational x@(CR x') = x' n % 2^n where n = digitsToBits digits + toRational _ = error "CReal.toRational" +instance RealFrac CReal where + properFraction x@(CR x') = (fromInteger n, x - fromInteger n) where n = x' 0 + +instance RealFloat CReal where + floatRadix _ = error "CCeal.floatRadix" + floatDigits _ = error "CReal.floatDigits" + floatRange _ = error "CReal.floatRange" + decodeFloat _ = error "CReal.decodeFloat" + encodeFloat _ _ = error "CReal.encodeFloat" + exponent _ = 0 + scaleFloat 0 x = x + significand x = x + isNaN _ = False + isInfinite _ = False + isDenormalized _ = False + isNegativeZero _ = False + isIEEE _ = False + +-- printing and reading the reals: + +-- |The 'showCReal' function connverts a 'CReal' to a 'String'. +showCReal :: Int -- ^ The number of decimals + -> CReal -- ^ The real number + -> String -- ^ The resulting string +showCReal d (CR x') + = (if s then "-" else "") ++ zs ++ (if d /= 0 then '.':fs' else "") + where b = digitsToBits d + n = x' b + ds = show (round_uk ((n*10^d) % 2^b)) + (s,ds') = let sgn = head ds == '-' in (sgn, if sgn then tail ds else ds) + ds'' = take (max (d+1-length ds') 0) (repeat '0') ++ ds' + (zs,fs) = splitAt (length ds'' -d) ds'' + fs' = case reverse $ dropWhile (== '0') $ reverse fs of + "" -> "0" + xs -> xs + +digitsToBits :: Int -> Int +digitsToBits d = ceiling (fromIntegral d * (logBase 2.0 10.0 :: Double)) + 4 + +digits :: Int +digits = 40 + +instance Read CReal where + readsPrec _p = readSigned readFloat + +instance Show CReal where + showsPrec p x = let xs = showCReal digits x in + if head xs == '-' then showParen (p > 6) (showString xs) + else showString xs + +-- GMP functions not provided by Haskell + +sizeinbase :: Integer -> Int -> Int +sizeinbase i b = f (abs i) + where f n = if n <= 1 then 1 else 1 + f (n `div` toInteger b) + +floorsqrt :: Integer -> Integer +floorsqrt x = until satisfy improve x + where improve y = floor ((y*y+x) % (2*y)) + satisfy y = y*y <= x && x <= (y+1)*(y+1) + +round_uk :: Rational -> Integer +round_uk x = floor (x+1 % 2)
+ Data/Number/Dif.hs view
@@ -0,0 +1,179 @@+-- | 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) => 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) => 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) => (Dif a -> Dif b) -> (a -> b) +deriv f = val . df . f . dVar + +-- |Convert a 'Dif' function to an ordinary function. +unDif :: (Num 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) => 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) => 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) => [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) => 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
+ Data/Number/Fixed.hs view
@@ -0,0 +1,129 @@+{-# OPTIONS_GHC -fglasgow-exts -fscoped-type-variables #-}+module Data.Number.Fixed(Fixed, Epsilon, Eps1, EpsDiv10, Prec10, Prec50, PrecPlus20,+ convertFixed, dynamicEps) where+import Numeric+import Data.Char+import Data.Ratio+import qualified Data.Number.FixedFunctions as F++class Epsilon e where+ eps :: e -> Rational++data Eps1+instance Epsilon Eps1 where+ eps _ = 1++data EpsDiv10 p+instance (Epsilon e) => Epsilon (EpsDiv10 e) where+ eps e = eps (un e) / 10+ where un :: EpsDiv10 e -> e+ un = undefined++data Prec10+instance Epsilon Prec10 where+ eps _ = 1e-10++data Prec50+instance Epsilon Prec50 where+ eps _ = 1e-50++data Prec500+instance Epsilon Prec500 where+ eps _ = 1e-500++data PrecPlus20 e+instance (Epsilon e) => Epsilon (PrecPlus20 e) where+ eps e = 1e-20 * eps (un e)+ where un :: PrecPlus20 e -> e+ un = undefined++-----------++newtype Fixed e = F Rational deriving (Eq, Ord, Enum, Real, RealFrac)++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) eps++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++-----------+++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
@@ -0,0 +1,469 @@+-- 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+ ++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
@@ -0,0 +1,44 @@+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/Symbolic.hs view
@@ -0,0 +1,167 @@+module Data.Number.Symbolic(Sym, var, con, subst, unSym) where++import Data.Char(isAlpha)+import Data.Maybe(fromMaybe)+import Debug.Trace++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')++var :: String -> Sym a+var s = App s undefined []++con :: a -> Sym a+con = Con++subst :: (Num 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'++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) => 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) => 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) => (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+ 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) => 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) => 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
@@ -0,0 +1,9 @@+{-# OPTIONS_GHC -fglasgow-exts #-}+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
+ Setup.hs view
@@ -0,0 +1,3 @@+module Main where+import Distribution.Simple+main = defaultMain
+ numbers.cabal view
@@ -0,0 +1,15 @@+Name: numbers+Version: 2007.4.29+License: BSD3+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+ precion fixed numbers, differentiable numbers, symbolic numbers.+Build-Depends: base+Exposed-modules: Data.Number.Symbolic Data.Number.Dif+ Data.Number.CReal Data.Number.Fixed+ Data.Number.Interval+Other-modules: Data.Number.Vectorspace Data.Number.FixedFunctions