module Numeric.VariablePrecision.Float
( VFloat()
, Normed(norm1, norm2, norm2Squared, normInfinity)
, effectivePrecisionWith
, effectivePrecision
, (-@?)
, DFloat(..)
, toDFloat
, fromDFloat
, withDFloat
) where
import Data.Data (Data())
import Data.Typeable (Typeable())
import Data.Bits (bit, shiftL, shiftR)
import Data.Ratio ((%), numerator, denominator)
import GHC.Float (showSignedFloat)
import Numeric (readSigned, readFloat)
import Text.FShow.RealFloat (DispFloat(), FShow(fshowsPrec), fshowFloat)
import Numeric.VariablePrecision.Algorithms
import Numeric.VariablePrecision.Precision
import Numeric.VariablePrecision.Precision.Reify
import Numeric.VariablePrecision.Integer.Logarithm
data VFloat p
= F !Integer !Int
| FZero
| FPosInf
| FNegInf
| FNaN
deriving (Data, Typeable)
encodeVFloat :: NaturalNumber p => VFloat p -> Integer -> Int -> VFloat p
encodeVFloat witness = self
where
b = fromIntegral $ precision (undefined `asTypeOf` witness)
b' = b 1
self 0 !_ = FZero
self m e = checkVFloat "encodeFloat'" $ encodeFloat'' (m > 0) m' (e sh) l
where
absm = abs m
m' = absm `shift` sh
e2 = integerLog2 absm
sh = b e2
l = integerLog2 m'
encodeFloat'' !s' !m' !e' !l
| m' <= 0 = failed
| b' == l = F (if s' then m' else negate m') e'
| b' < l = encodeFloat'' s' (m' `shiftR` 1) (e' + 1) (l 1)
| b' > l = encodeFloat'' s' (m' `shiftL` 1) (e' 1) (l + 1)
| otherwise = failed
where
failed = error $ "Numeric.VariablePrecision.Float.encodeVFloat: internal error (please report this bug): "
++ show (b, b', l, s', m', e')
instance NaturalNumber p => DispFloat (VFloat p) where
instance NaturalNumber p => FShow (VFloat p) where
fshowsPrec p = showSignedFloat fshowFloat p
instance NaturalNumber p => Show (VFloat p) where
showsPrec = fshowsPrec
instance NaturalNumber p => Read (VFloat p) where
readsPrec _ = readSigned readFloat
instance HasPrecision VFloat
minimumExponent, maximumExponent :: Int
minimumExponent = negate (bit 20)
maximumExponent = bit 20
asTypeIn :: (a -> b) -> a
asTypeIn _ = undefined
asTypeOut :: (a -> b) -> b
asTypeOut _ = undefined
asTypeOut2 :: (a -> b -> c) -> c
asTypeOut2 _ = undefined
instance VariablePrecision VFloat where
adjustPrecision = self
where
p = asTypeIn self
q = asTypeOut self
np = floatDigits p
nq = floatDigits q
n = nq np
self FZero = FZero
self FPosInf = FPosInf
self FNegInf = FNegInf
self FNaN = FNaN
self (F m e)
| n > 0 = encodeVFloat q (m `shiftL` n) (e n)
| n == 0 = encodeVFloat q m e
| n < 0 = encodeVFloat q (m `shiftR` negate n) (e + negate n)
| otherwise = unreachable
instance Eq (VFloat p) where
FZero == FZero = True
F a b == F x y = a == x && b == y
FPosInf == FPosInf = True
FNegInf == FNegInf = True
_ == _ = False
a /= x = not (a == x)
instance Ord (VFloat p) where
FZero < FZero = False
FZero < F x _ = 0 < x
F a _ < FZero = a < 0
F a b < F x y
| a > 0 && x > 0 && b < y = True
| a > 0 && x > 0 && b == y = a < x
| a > 0 && x > 0 && b > y = False
| a > 0 && x < 0 = False
| a < 0 && x > 0 = True
| a < 0 && x < 0 && b < y = False
| a < 0 && x < 0 && b == y = a < x
| a < 0 && x < 0 && b > y = True
| otherwise = unreachable
FNaN < _ = False
_ < FNaN = False
FPosInf < _ = False
_ < FPosInf = True
_ < FNegInf = False
FNegInf < _ = True
a > x = x < a
a <= x = a < x || a == x
a >= x = a > x || a == x
min a@FNaN !_ = a
min !_ x@FNaN = x
min a x
| a <= x = a
| otherwise = x
max a@FNaN !_ = a
max !_ x@FNaN = x
max a x
| a >= x = a
| otherwise = x
instance NaturalNumber p => Num (VFloat p) where
f@(F a b) + F x y
| b > y = encodeVFloat f (a + (x `shiftR` (b y))) b
| b == y = encodeVFloat f (a + x) b
| b < y = encodeVFloat f ((a `shiftR` (y b)) + x) y
| otherwise = unreachable
a@FNaN + _ = a
_ + x@FNaN = x
FZero + x = x
a + FZero = a
FPosInf + FNegInf = FNaN
FNegInf + FPosInf = FNaN
FPosInf + _ = FPosInf
_ + FPosInf = FPosInf
FNegInf + _ = FNegInf
_ + FNegInf = FNegInf
f@(F a b) F x y
| b > y = encodeVFloat f (a (x `shiftR` (b y))) b
| b == y = encodeVFloat f (a x) b
| b < y = encodeVFloat f ((a `shiftR` (y b)) x) y
| otherwise = unreachable
a@FNaN _ = a
_ x@FNaN = x
FZero x = negate x
a FZero = a
FPosInf FPosInf = FNaN
FNegInf FNegInf = FNaN
FPosInf _ = FPosInf
_ FPosInf = FNegInf
FNegInf _ = FNegInf
_ FNegInf = FPosInf
negate (F a b) = checkVFloat "negate" $ F (negate a) b
negate FZero = FZero
negate FPosInf = FNegInf
negate FNegInf = FPosInf
negate a@FNaN = a
abs !a
| a < 0 = negate a
| otherwise = a
signum !a
| a < 0 = 1
| a > 0 = 1
| otherwise = a
f@(F a b) * F x y = encodeVFloat f ((a * x) `shiftR` (k 1)) (b + y + k 1) where k = fromIntegral $ precision f
a@FNaN * _ = a
_ * x@FNaN = x
FZero * FPosInf = FNaN
FZero * FNegInf = FNaN
FZero * _ = FZero
FPosInf * FZero = FNaN
FNegInf * FZero = FNaN
_ * FZero = FZero
a * x
| sameSign a x = FPosInf
| otherwise = FNegInf
fromInteger !i = encodeFloat i 0
instance NaturalNumber p => Real (VFloat p) where
toRational FZero = 0
toRational (F m e)
| e > 0 = fromInteger (m `shiftL` e)
| e == 0 = fromInteger m
| e < 0 = m % bit (negate e)
| otherwise = unreachable
toRational FPosInf = 1 % 0
toRational FNegInf = (1) % 0
toRational FNaN = 0 % 0
instance NaturalNumber p => Fractional (VFloat p) where
f@(F _ _) / g@(F _ _) = f * recip g
a@FNaN / _ = a
_ / x@FNaN = x
FPosInf / FPosInf = FNaN
FPosInf / FNegInf = FNaN
FNegInf / FPosInf = FNaN
FNegInf / FNegInf = FNaN
_ / FPosInf = FZero
_ / FNegInf = FZero
a / FZero
| a > 0 = FPosInf
| a < 0 = FNegInf
| otherwise = FNaN
FZero / _ = FZero
a / x
| a `sameSign` x = FPosInf
| otherwise = FNegInf
recip a@FNaN = a
recip FZero = FPosInf
recip FPosInf = FZero
recip FNegInf = FZero
recip f@(F m e) = encodeVFloat f (bit k `quot` m) (negate (k + e)) where k = 2 * fromIntegral (precision f)
fromRational r = fromInteger (numerator r) / fromInteger (denominator r)
instance NaturalNumber p => RealFrac (VFloat p) where
properFraction = self
where
p = fromIntegral $ precision (asTypeIn self)
self FZero = (0, FZero)
self me@(F m e)
| e >= 0 = (fromInteger m, FZero)
| e < negate p = (0, me)
| otherwise = (fromInteger n', f')
where
n = m `shiftR` (negate e)
d = checkVFloat "properFraction" $ F (n `shiftL` (negate e)) e
f = me d
(n', f')
| (m >= 0) == (f >= 0) = (n, f)
| otherwise = (n + 1, f 1)
self f = (error $ "Numeric.VariablePrecision.Float.properFraction: not finite: " ++ show f, f)
instance NaturalNumber p => RealFloat (VFloat p) where
floatRadix _ = 2
floatDigits = self
where
prec = fromIntegral $ precision (asTypeIn self)
self = const prec
floatRange = const (minimumExponent, maximumExponent)
isNaN FNaN = True
isNaN _ = False
isInfinite FPosInf = True
isInfinite FNegInf = True
isInfinite _ = False
isDenormalized _ = False
isNegativeZero _ = False
isIEEE _ = False
decodeFloat FZero = (0, 0)
decodeFloat (F m e) = (m, e)
decodeFloat f = error $ "Numeric.VariablePrecision.Float.decodeFloat: not finite: " ++ show f
encodeFloat = self
where
self = encodeVFloat (undefined `asTypeOf` asTypeOut2 self)
exponent = self
where
prec = fromIntegral $ precision (asTypeIn self)
self FZero = 0
self (F _ e) = e + prec
self f = error $ "Numeric.VariablePrecision.Float.exponent: not finite: " ++ show f
significand = self
where
prec = fromIntegral $ precision (asTypeIn self)
e = negate prec
self (F m _) = checkVFloat "significand" $ F m e
self f = f
scaleFloat n (F m e) = checkVFloat "scaleFloat" $ F m (e + n)
scaleFloat _ f = f
shift :: Integer -> Int -> Integer
shift !n !k
| k > 0 = n `shiftL` k
| k == 0 = n
| k < 0 = n `shiftR` (negate k)
| otherwise = unreachable
instance NaturalNumber p => Floating (VFloat p) where
pi = genericPi 2
sqrt = genericSqrt 2
exp = genericExp 2
log = self
where
log2 = genericLog2 2
self = genericLog' 2 log2
sin = viaDouble sin
cos = viaDouble cos
tan = viaDouble tan
sinh = viaDouble sinh
cosh = viaDouble cosh
tanh = viaDouble tanh
asin = viaDouble asin
acos = viaDouble acos
atan = viaDouble atan
asinh = viaDouble asinh
acosh = viaDouble acosh
atanh = viaDouble atanh
checkVFloat :: NaturalNumber p => String -> VFloat p -> VFloat p
checkVFloat = self
where
prec = fromIntegral $ precision (asTypeOut2 self)
prec' = prec 1
elo = minimumExponent
ehi = maximumExponent
self s x@(F m e)
| not mok = error $ "Numeric.VariablePrecision.Float.checkVFloat." ++ s ++ ": internal error (please report this bug): " ++ show ((m, am, lm, prec, prec', mok), (elo, e, ehi, eok))
| eok = x
| e < elo = FZero
| m > 0 = FPosInf
| m < 0 = FNegInf
| otherwise = unreachable
where
eok = elo <= e && e <= ehi
mok = lm == prec'
lm = integerLog2 am
am = abs m
self _ x = x
class HasPrecision t => Normed t where
norm1 :: NaturalNumber p => t p -> VFloat p
norm2 :: NaturalNumber p => t p -> VFloat p
norm2Squared :: NaturalNumber p => t p -> VFloat p
normInfinity :: NaturalNumber p => t p -> VFloat p
instance Normed VFloat where
norm1 = abs
norm2 = abs
norm2Squared x = x * x
normInfinity = abs
effectivePrecisionWith :: (Num t, RealFloat r) => (t -> r) -> t -> t -> Int
effectivePrecisionWith n i j
| t a && t b && t c = p (d `max` (e d))
| otherwise = 0
where
t k = k > 0 && not (isInfinite k)
d = (x `max` y) z
e = abs (x y) `min` p
p = floatDigits a
x = exponent a
y = exponent b
z = exponent c
a = n i
b = n j
c = n (i j)
effectivePrecision :: (NaturalNumber p, HasPrecision t, Normed t, Num (t p)) => t p -> t p -> Int
effectivePrecision = effectivePrecisionWith normInfinity
infix 6 `effectivePrecision`
(-@?) :: (NaturalNumber p, HasPrecision t, Normed t, Num (t p)) => t p -> t p -> Int
(-@?) = effectivePrecision
infix 6 -@?
unreachable :: a
unreachable = error "Numeric.VariablePrecision.Float: internal error (please report this bug): unreachable code was reached"
data DFloat
= DFloat { dPrecision :: !Word, dMantissa :: !Integer, dExponent :: !Int }
| DZero { dPrecision :: !Word }
| DPositiveInfinity { dPrecision :: !Word }
| DNegativeInfinity { dPrecision :: !Word }
| DNotANumber { dPrecision :: !Word }
deriving (Eq, Ord, Read, Show, Data, Typeable)
toDFloat :: NaturalNumber p => VFloat p -> DFloat
toDFloat f@(F m e) = DFloat (precision f) m e
toDFloat f@FZero = DZero (precision f)
toDFloat f@FPosInf = DPositiveInfinity (precision f)
toDFloat f@FNegInf = DNegativeInfinity (precision f)
toDFloat f@FNaN = DNotANumber (precision f)
fromDFloat :: NaturalNumber p => DFloat -> Maybe (VFloat p)
fromDFloat d
| dPrecision d == precision result = Just result
| otherwise = Nothing
where
result = case d of
DFloat _ m e -> encodeVFloat undefined m e
DZero _ -> FZero
DPositiveInfinity _ -> FPosInf
DNegativeInfinity _ -> FNegInf
DNotANumber _ -> FNaN
withDFloat :: DFloat -> (forall p . NaturalNumber p => VFloat p -> r) -> r
withDFloat (DFloat p m e) f = reifyPrecision p $ \prec -> f (encodeVFloat undefined m e `atPrecision` prec)
withDFloat d f = unsafeWithDFloat d f
unsafeWithDFloat :: DFloat -> (forall p . NaturalNumber p => VFloat p -> r) -> r
unsafeWithDFloat (DFloat p m e) f = reifyPrecision p $ \prec -> f (F m e `atPrecision` prec)
unsafeWithDFloat (DZero p) f = reifyPrecision p $ \prec -> f (FZero `atPrecision` prec)
unsafeWithDFloat (DPositiveInfinity p) f = reifyPrecision p $ \prec -> f (FPosInf `atPrecision` prec)
unsafeWithDFloat (DNegativeInfinity p) f = reifyPrecision p $ \prec -> f (FNegInf `atPrecision` prec)
unsafeWithDFloat (DNotANumber p) f = reifyPrecision p $ \prec -> f (FNaN `atPrecision` prec)