{-# LANGUAGE DeriveDataTypeable #-} {-| Module : Data.Number.ER.Real.Base Description : arbitrary precision floating point numbers Copyright : (c) Michal Konecny License : BSD3 Maintainer : mik@konecny.aow.cz Stability : experimental Portability : portable This module defines an arbitrary precision floating point type and its operations. It should be viewed more abstractly as an instance of 'B.ERRealBase' when used as interval endpoints. -} module Data.Number.ER.Real.Base.Float ( ERFloat ) where import qualified Data.Number.ER.ExtendedInteger as EI import Data.Number.ER.PlusMinus import Data.Number.ER.Misc import Data.Number.ER.BasicTypes import qualified Data.Number.ER.Real.Base as B import Data.Ratio import Data.Typeable import Data.Generics.Basics import Data.Binary -- import BinaryDerive --debugMsg = unsafePrint debugMsg msg = id {-| A floating point number with a given but arbitrary precision represented by its 'Granularity'. * base: 2. * granularity specifies the bit-size of both the significand and the exponent * special values: NaN, signed Infinity and signed Zero * no denormalised numbers * operations unify the granularity of their operands to the maximum 'Granularity' * Rounding is always towards +Infinity. For field operations, the rounded result is as close as possible to the exact result. -} data ERFloat = ERFloatNaN -- any number / bottom { apfltGran :: Granularity -- >0 } | ERFloatInfty { apfltGran :: Granularity, -- >0 apfltSign :: PlusMinus } | ERFloatZero { apfltGran :: Granularity, -- >0 apfltSign :: PlusMinus } | ERFloat { -- represents: -- sign * (1 + (mant/2^gran)) * (2 ^ exp) apfltGran :: Granularity, -- >0 granularity apfltSign :: PlusMinus, apfltMant :: Integer, -- 0 .. (2^gran - 1) apfltExp :: Integer -- -2^gran..2^gran } deriving (Typeable, Data) zero = ERFloatZero 10 Plus {- the following has been generated by BinaryDerive -} instance Binary ERFloat where put (ERFloatNaN a) = putWord8 0 >> put a put (ERFloatInfty a b) = putWord8 1 >> put a >> put b put (ERFloatZero a b) = putWord8 2 >> put a >> put b put (ERFloat a b c d) = putWord8 3 >> put a >> put b >> put c >> put d get = do tag_ <- getWord8 case tag_ of 0 -> get >>= \a -> return (ERFloatNaN a) 1 -> get >>= \a -> get >>= \b -> return (ERFloatInfty a b) 2 -> get >>= \a -> get >>= \b -> return (ERFloatZero a b) 3 -> get >>= \a -> get >>= \b -> get >>= \c -> get >>= \d -> return (ERFloat a b c d) _ -> fail "no parse" {- the above has been generated by BinaryDerive -} {-| normalisation * ensures that the components are within their regions * possibly turning the number into a zero or infinity -} normaliseERFloat :: ERFloat -> ERFloat normaliseERFloat flt@(ERFloat gr s m e) | m < 0 = normaliseERFloat $ ERFloat gr s (2*m + grmax) (e - 1) | m >= grmax = normaliseERFloat $ ERFloat gr s ((m - grmax + (rndCorr s)) `div` 2) (e + 1) | e > grmax = case s of Plus -> ERFloatInfty gr Plus Minus -> minERFloat gr -- round upwards! | e < -grmax = case s of Plus -> ulpERFloat gr -- round upwards! Minus -> ERFloatZero gr Minus | otherwise = flt where grmax = 2^gr normaliseERFloat flt = flt ulpERFloat gr = ERFloat gr Plus 0 (-2^gr) minERFloat gr = ERFloat gr Minus (grmax - 1) grmax where grmax = 2^gr maxERFloat gr = ERFloat gr Plus (grmax - 1) grmax where grmax = 2^gr rndCorr Plus = 1 rndCorr Minus = 0 increaseERFloatExp e flt@(ERFloat gr s m eOld) = ERFloat gr s mNew e where mNew = -grmax + ((m + grmax + (rndCorr s) * (ediff - 1)) `div` ediff) -- .^^^^^^^^^^^^^^^^^^^^^^^^^ round upwards grmax = 2^gr ediff = 2^(e - eOld) -- assuming e >= eOld increaseERFloatExp _ flt = flt decreaseERFloatExp e flt@(ERFloat gr s m eOld) = ERFloat gr s mNew e where mNew = -grmax + ediff * (m + grmax) grmax = 2^gr ediff = 2^(eOld - e) -- assuming e <= eOld decreaseERFloatExp _ flt = flt apFloatExponent :: ERFloat -> EI.ExtendedInteger apFloatExponent (ERFloatInfty _ _) = EI.PlusInfinity apFloatExponent (ERFloatZero _ _) = EI.MinusInfinity apFloatExponent (ERFloatNaN _) = EI.PlusInfinity -- includes infinity apFloatExponent flt = EI.Finite $ apfltExp flt setERFloatGranularity :: Granularity -> ERFloat -> ERFloat setERFloatGranularity gr flt@(ERFloat oldGr s m e) | gr > 0 = normaliseERFloat $ ERFloat gr s newM e | otherwise = flt where newM = (m * (2^gr) + ((rndCorr s)*(2^oldGr - 1))) -- round upwards! `div` (2^oldGr) setERFloatGranularity gr f = f { apfltGran = gr } setERFloatMinGranularity :: Granularity -> ERFloat -> ERFloat setERFloatMinGranularity gr flt | gr > oldGr = setERFloatGranularity gr flt | otherwise = flt where oldGr = apfltGran flt apfltGranularity = apfltGran {-^ see the documentation of 'ERRealBase.getBaseMaxRounding' -} apfltGetMaxRounding :: ERFloat -> ERFloat apfltGetMaxRounding f@(ERFloatNaN _) = f apfltGetMaxRounding f@(ERFloatInfty _ _) = f apfltGetMaxRounding (ERFloatZero gr _) = ERFloat gr Plus 0 (-(2^gr)) apfltGetMaxRounding (ERFloat gr s m e) = ERFloat gr Plus 0 (max (e - (toInteger gr)) (-(2^gr))) instance Show ERFloat where show = showERFloat 6 True False showERFloat numDigits showGran showComponents flt = showERF flt where maybeGran gr | showGran = "{g=" ++ show gr ++ "}" | otherwise = "" showPM Plus = "" showPM Minus = "-" showERF (ERFloatNaN gr) = "NaN" ++ (maybeGran gr) showERF (ERFloatZero gr pm) = showPM pm ++ "0.0" ++ (maybeGran gr) showERF (ERFloatInfty gr pm) = showPM pm ++ "oo" ++ (maybeGran gr) showERF f@(ERFloat gr s m e) = decimal ++ (maybeGran gr) ++ maybeComps where maybeComps | showComponents = "{val="++ show (s,m,e) ++ "}" | otherwise = "" decimal = showPM s ++ show digit1 ++ "." ++ (concat $ map show $ take numDigits digits) ++ (if dexp == 0 then "" else "e" ++ show dexp) (dexp, digit1 : digits) | noLeadingZerosDexp == -1 = (0, 0 : noLeadingZeroDigits) | otherwise = (noLeadingZerosDexp, noLeadingZeroDigits) noLeadingZerosDexp = dexpBound - zerosCount (zerosCount, noLeadingZeroDigits) = stripCountZeros 0 preDigits stripCountZeros prevZeros ds@(d : dsRest) | d == 0 = stripCountZeros (prevZeros + 1) dsRest | otherwise = (prevZeros, ds) dexpBound -- upper bound of dexp: f/10^dexpBound < 1 | e >= 0 = intLogUp 10 (2^e) | e < 0 = 2 - (intLogUp 10 (2^(-e))) preDigits = getDigits $ (abs $ setERFloatGranularity gran f) / (ten ^^ dexpBound) ten = setERFloatGranularity gran 10 gran = 10 + (max (4 * numDigits) gr) getDigits ff = digit : digits where digit :: Integer digit = truncate ff digits = getDigits ((ff - (fromInteger digit)) * ten) {- Beware: cannot use List.elem with ERFloat because of the intensional nature of Eq (eg ERFloatNaN /= ERFloatNaN). -} instance Eq ERFloat where (ERFloatNaN _) == _ = False -- error "cannot compare NaN" _ == (ERFloatNaN _) = False -- error "cannot compare NaN" (ERFloatZero _ _) == (ERFloatZero _ _) = True (ERFloatInfty _ pm1) == (ERFloatInfty _ pm2) = (pm1 == pm2) f1@(ERFloat gr1 s1 m1 e1) == f2@(ERFloat gr2 s2 m2 e2) | gr1 < gr2 = (setERFloatGranularity gr2 f1) == f2 | gr1 > gr2 = f1 == (setERFloatGranularity gr1 f2) | otherwise = s1 == s2 && m1 == m2 && e1 == e2 _ == _ = False isERFloatNaN (ERFloatNaN _) = True isERFloatNaN _ = False instance Ord ERFloat where {- compare NaN -} compare a b@(ERFloatNaN _) = unsafePrint ("ERFloat: comparing NaN: " ++ show a ++ " vs. " ++ show b) EQ -- error $ "ERFloat: comparing NaN: " ++ show a ++ " vs. " ++ show b compare a@(ERFloatNaN _) b = unsafePrint ("ERFloat: comparing NaN: " ++ show a ++ " vs. " ++ show b) EQ -- error $ "ERFloat: comparing NaN: " ++ show a ++ " vs. " ++ show b {- compare infty -} compare (ERFloatInfty gr1 pm1) (ERFloatInfty gr2 pm2) = compare pm1 pm2 compare _ (ERFloatInfty _ Plus) = LT compare _ (ERFloatInfty _ Minus) = GT compare (ERFloatInfty _ Plus) _ = GT compare (ERFloatInfty _ Minus) _ = LT {- compare zero -} compare (ERFloatZero gr1 pm1) (ERFloatZero gr2 pm2) = EQ compare (ERFloatZero _ _) (ERFloat _ Plus _ _) = LT compare (ERFloatZero _ _) (ERFloat _ Minus _ _) = GT compare (ERFloat _ Minus _ _) (ERFloatZero _ _) = LT compare (ERFloat _ Plus _ _) (ERFloatZero _ _) = GT {- compare regular -} compare (ERFloat _ Minus _ _) (ERFloat _ Plus _ _) = LT compare (ERFloat _ Plus _ _) (ERFloat _ Minus _ _) = GT compare (ERFloat gr1 Plus m1 e1) (ERFloat gr2 _ m2 e2) | e1 < e2 = LT | e1 > e2 = GT | gr1 == gr2 = compare m1 m2 | otherwise = compare ((2^gr2)*m1) ((2^gr1)*m2) compare f1@(ERFloat _ Minus _ _) f2@(ERFloat _ _ _ _) = compare (-f2) (-f1) instance Num ERFloat where fromInteger n | n == 0 = ERFloatZero (B.defaultGranularity zero) Plus | n < 0 = normaliseERFloat $ ERFloat gr Minus m e | otherwise = normaliseERFloat $ ERFloat gr Plus m e where gr = fromInteger e e = max (toInteger (B.defaultGranularity zero)) $ (intLogUp 2 $ abs n) - 1 m = (abs n) - 2^gr abs f@(ERFloatNaN _) = f abs f = f { apfltSign = Plus } signum f@(ERFloatNaN _) = f signum (ERFloatZero gr Plus) = setERFloatMinGranularity gr 1 signum (ERFloatZero gr Minus) = setERFloatMinGranularity gr (-1) signum (ERFloatInfty gr Plus) = setERFloatMinGranularity gr 1 signum (ERFloatInfty gr Minus) = setERFloatMinGranularity gr (-1) signum flt = case apfltSign flt of { Plus -> 1; Minus -> -1 } negate (ERFloat gr s m e) = ERFloat gr (signNeg s) m e negate (ERFloatZero gr s) = ERFloatZero gr (signNeg s) negate (ERFloatInfty gr s) = ERFloatInfty gr (signNeg s) negate f@(ERFloatNaN _) = f {- addition -} f1 + f2 -- ensure equal granularity: | gr1 > gr2 = f1 + (setERFloatGranularity gr1 f2) | gr1 < gr2 = (setERFloatGranularity gr2 f1) + f2 where gr1 = apfltGran f1 gr2 = apfltGran f2 f@(ERFloatNaN _) + _ = f _ + f@(ERFloatNaN _) = f (ERFloatZero _ _) + f = f f + (ERFloatZero _ _) = f (ERFloatInfty gr Plus) + (ERFloatInfty _ Minus) = debugMsg ("ERFloat: infty - infty -> NaN\n") $ ERFloatNaN gr (ERFloatInfty gr Minus) + (ERFloatInfty _ Plus) = debugMsg ("ERFloat: -infty + infty -> NaN\n") $ ERFloatNaN gr f@(ERFloatInfty gr s) + _ = f _ + f@(ERFloatInfty gr s) = f f1@(ERFloat gr s1 m1 e1) + f2@(ERFloat _ s2 m2 e2) -- equalise the exponents: | e1 < e2 = f1 + (decreaseERFloatExp e1 f2) | e1 > e2 = (decreaseERFloatExp e2 f1) + f2 -- ensure positive comes before negative: | s1 == Minus && s2 == Plus = f2 + f1 -- opposite signs: | s1 == Plus && s2 == Minus && m1 == m2 = ERFloatZero gr Plus | s1 == Plus && s2 == Minus && m1 > m2 = normaliseERFloat $ ERFloat gr s1 (m1 - m2 - 2^gr) e1 | s1 == Plus && s2 == Minus && m1 < m2 = normaliseERFloat $ ERFloat gr s2 (m2 - m1 - 2^gr) e1 -- equal signs: | otherwise = normaliseERFloat $ ERFloat gr s1 (m1 + m2 + 2^gr) e1 {- multiplication -} -- ensure equal granularity: f1 * f2 | gr1 > gr2 = f1 * (setERFloatGranularity gr1 f2) | gr1 < gr2 = (setERFloatGranularity gr2 f1) * f2 where gr1 = apfltGran f1 gr2 = apfltGran f2 -- NaN: f@(ERFloatNaN _) * _ = f _ * f@(ERFloatNaN _) = f -- Infty (ERFloatInfty gr _) * (ERFloatZero _ _) = debugMsg ("ERFloat: infty * 0 -> NaN\n") $ ERFloatNaN gr (ERFloatZero gr _) * (ERFloatInfty _ _) = debugMsg ("ERFloat: 0 * infty -> NaN\n") $ ERFloatNaN gr f * (ERFloatInfty gr s) = ERFloatInfty gr $ signMult s (apfltSign f) (ERFloatInfty gr s) * f = ERFloatInfty gr $ signMult s (apfltSign f) -- Zero (ERFloatZero gr s) * f = ERFloatZero gr $ signMult s (apfltSign f) f * (ERFloatZero gr s) = ERFloatZero gr $ signMult s (apfltSign f) -- regular f1@(ERFloat gr s1 m1 e1) * f2@(ERFloat _ s2 m2 e2) = normaliseERFloat $ ERFloat gr s mNew (e1 + e2) where s = signMult s1 s2 mNew = m1 + m2 + ((m1 * m2 + (rndCorr s) * (2^gr - 1)) `div` 2^gr) instance Fractional ERFloat where fromRational rat = -- error "ERFloat: fromRational cannot be implemented reliably: use apfloatFromRational instead" (fromInteger $ numerator rat) / (fromInteger $ denominator rat) f1 / f2 | gr1 > gr2 = f1 / (setERFloatGranularity gr1 f2) | gr1 < gr2 = (setERFloatGranularity gr2 f1) / f2 where gr1 = apfltGran f1 gr2 = apfltGran f2 f@(ERFloatNaN _) / _ = f f1 / f2 = case apfltSign f1 of Plus -> f1 * (recip f2) Minus -> (- f1) * (recip (- f2)) -- rounding upwards! recip f@(ERFloatNaN _) = f recip (ERFloatZero gr s) = ERFloatInfty gr s recip (ERFloatInfty gr s) = ERFloatZero gr s recip (ERFloat gr s m e) = normaliseERFloat $ ERFloat gr s mNew (-e) where mNew = (- grmax * m + (rndCorr s) * (grmax + m -1)) -- rounding upwards! `div` (grmax + m) grmax = 2^gr apfloatFromRational :: Granularity -> Rational -> ERFloat apfloatFromRational gran rat = (setERFloatMinGranularity gran (fromInteger $ numerator rat)) / (fromInteger $ denominator rat) instance Real ERFloat where toRational (ERFloat gr s m e) = case s of Plus -> r Minus -> -r where r = (eOn2R) * (1 + mR/(grOn2R)) mR = toRational m eOn2R = toRational $ 2 ^^ e grOn2R = toRational $ 2 ^ gr toRational (ERFloatZero _ _) = 0 toRational f = error $ "cannot covert " ++ show f ++ " to a rational" instance RealFrac ERFloat where properFraction (ERFloatNaN _) = error "no integral part in ERFloatNaN" properFraction (ERFloatZero _ _) = (0, 0) properFraction (ERFloatInfty _ _) = error "no integral part in ERFloatInfty" properFraction f@(ERFloat gr s m e) | e < 0 = (0, f) | s == Plus = (n, frac) | s == Minus = (-n, frac) where n = fromInteger $ 2^e + (m*(2^e) `div` 2^gr) frac = f - (fromInteger $ toInteger n) instance B.ERRealBase ERFloat where typeName _ = "ERFloat (pure Haskell implementation)" defaultGranularity _ = 10 getApproxBinaryLog = apFloatExponent getGranularity = apfltGran setMinGranularity = setERFloatMinGranularity setGranularity = setERFloatGranularity getMaxRounding = apfltGetMaxRounding isERNaN (ERFloatNaN _) = True isERNaN _ = False erNaN = ERFloatNaN (B.defaultGranularity zero) isPlusInfinity (ERFloatInfty _ Plus) = True isPlusInfinity _ = False plusInfinity = ERFloatInfty (B.defaultGranularity zero) Plus fromDouble d | isNaN d = ERFloatNaN (B.defaultGranularity zero) | otherwise = (fromRational . toRational) d toDouble (ERFloatInfty _ s) = signToNum s * (1/0) toDouble (ERFloatNaN _) = 0/0 toDouble flt = (fromInteger $ numerator rat) / (fromInteger $ denominator rat) where rat = toRational flt fromFloat f | isNaN f = ERFloatNaN (B.defaultGranularity zero) | otherwise = (fromRational . toRational) f toFloat (ERFloatInfty _ s) = signToNum s * (1/0) toFloat (ERFloatNaN _) = 0/0 toFloat flt = (fromInteger $ numerator rat) / (fromInteger $ denominator rat) where rat = toRational flt showDiGrCmp = showERFloat