Permalink
Browse files

initial commit

  • Loading branch information...
0 parents commit c40a53af96d4cf818b063e02ec8f6a4089bbeb2a @DanBurton DanBurton committed Feb 21, 2012
@@ -0,0 +1,105 @@
+-- | 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, Ord)
+
+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) => Fractional (BigFloat e) where
+ recip (BF m e) = bf (base / m) (-(e + 1))
+ -- Take care not to lose precision for small numbers
+ fromRational x = if abs x < 1 then recip $ bf (fromRational (recip x)) 0
+ else bf (fromRational 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
@@ -0,0 +1,239 @@
+{-# 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)
+
+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)
Oops, something went wrong.

0 comments on commit c40a53a

Please sign in to comment.