module Numeric.QD.QuadDouble
( QuadDouble(QuadDouble)
, toDouble
, fromDouble
, toDoubleDouble
, fromDoubleDouble
, sqr
) where
import Foreign (Ptr, alloca, castPtr, Storable(..), unsafePerformIO, with)
import Foreign.C (CDouble, CInt)
import Data.Bits (shiftL, shiftR)
import Data.Ratio ((%))
import Data.Typeable (Typeable(..))
import Numeric (showFloat, readSigned, readFloat)
import Numeric.QD.DoubleDouble (DoubleDouble(DoubleDouble))
import Numeric.QD.QuadDouble.Raw
( c_qd_add
, c_qd_sub
, c_qd_mul
, c_qd_div
, c_qd_pi
, c_qd_exp
, c_qd_sqrt
, c_qd_log
, c_qd_sin
, c_qd_cos
, c_qd_tan
, c_qd_asin
, c_qd_acos
, c_qd_atan
, c_qd_sinh
, c_qd_cosh
, c_qd_tanh
, c_qd_asinh
, c_qd_acosh
, c_qd_atanh
, c_qd_comp
, c_qd_neg
, c_qd_abs
, c_qd_aint
, c_qd_nint
, c_qd_ceil
, c_qd_floor
, c_qd_atan2
, c_qd_sqr
)
data QuadDouble = QuadDouble !CDouble !CDouble !CDouble !CDouble deriving Typeable
toDouble :: QuadDouble -> Double
toDouble !(QuadDouble a _ _ _) = realToFrac a
fromDouble :: Double -> QuadDouble
fromDouble !a = QuadDouble (realToFrac a) 0 0 0
toDoubleDouble :: QuadDouble -> DoubleDouble
toDoubleDouble !(QuadDouble a b _ _) = DoubleDouble a b
fromDoubleDouble :: DoubleDouble -> QuadDouble
fromDoubleDouble !(DoubleDouble a b) = QuadDouble a b 0 0
instance Eq QuadDouble where
(!a) == (!b) = a `compare` b == EQ
(!a) /= (!b) = a `compare` b /= EQ
instance Ord QuadDouble where
(!a) `compare` (!b) = unsafePerformIO $ with a $ \p -> with b $ \q -> alloca $ \r -> do
c_qd_comp (castPtr p) (castPtr q) (castPtr r)
!i <- peek r
return $ i `compare` (0 :: CInt)
instance Num QuadDouble where
(+) = lift_qd_qd_qd c_qd_add
(*) = lift_qd_qd_qd c_qd_mul
() = lift_qd_qd_qd c_qd_sub
negate = lift_qd_qd c_qd_neg
abs = lift_qd_qd c_qd_abs
signum !a = case a `compare` 0 of { LT -> 1 ; EQ -> 0 ; GT -> 1 }
fromInteger !i = fromRational (i % 1)
sqr :: QuadDouble -> QuadDouble
sqr = lift_qd_qd c_qd_sqr
instance Fractional QuadDouble where
(/) = lift_qd_qd_qd c_qd_div
recip !b = 1 / b
fromRational !k = let a = fromRational k
ka = k toRational a
b = fromRational ka
kb = ka toRational b
c = fromRational kb
kc = kb toRational c
d = fromRational kc
in QuadDouble a b c d
instance Real QuadDouble where
toRational (QuadDouble a b c d) = toRational a + toRational b + toRational c + toRational d
instance RealFrac QuadDouble where
properFraction k = let (n, r) = properFraction (toRational k)
in (n, fromRational r)
truncate = truncate . toRational . lift_qd_qd c_qd_aint
round = round . toRational . lift_qd_qd c_qd_nint
ceiling = ceiling . toRational . lift_qd_qd c_qd_ceil
floor = floor . toRational . lift_qd_qd c_qd_floor
instance Floating QuadDouble where
pi = unsafePerformIO $ alloca $ \r -> c_qd_pi (castPtr r) >> peek r
exp = lift_qd_qd c_qd_exp
sqrt = lift_qd_qd c_qd_sqrt
log = lift_qd_qd c_qd_log
sin = lift_qd_qd c_qd_sin
cos = lift_qd_qd c_qd_cos
tan = lift_qd_qd c_qd_tan
asin = lift_qd_qd c_qd_asin
acos = lift_qd_qd c_qd_acos
atan = lift_qd_qd c_qd_atan
sinh = lift_qd_qd c_qd_sinh
cosh = lift_qd_qd c_qd_cosh
tanh = lift_qd_qd c_qd_tanh
asinh = lift_qd_qd c_qd_asinh
acosh = lift_qd_qd c_qd_acosh
atanh = lift_qd_qd c_qd_atanh
instance RealFloat QuadDouble where
floatRadix _ = 2
floatDigits _ = 4 * floatDigits (error "Numeric.QD.QuadDouble.floatDigits" :: CDouble)
floatRange _ = floatRange (error "Numeric.QD.QuadDouble.floatRange" :: CDouble)
decodeFloat !x = case toRational x of
0 -> (0, 0)
r -> let k = floor $ fromIntegral ff logBase (fromIntegral $ floatRadix x) (abs x)
i = round $ (fromIntegral $ floatRadix x) ^ k * r
fixup m e =
if abs m < mMin
then fixup (m `shiftL` 1) (e 1)
else if abs m >= mMax
then fixup (m `shiftR` 1) (e + 1)
else (m, e)
mMin = 1 `shiftL` (ff 1)
mMax = 1 `shiftL` ff
ff = floatDigits x
g = k
in fixup i g
encodeFloat m e = scaleFloat e (fromInteger m)
scaleFloat !n !(QuadDouble a b c d) = QuadDouble (scaleFloat n a) (scaleFloat n b) (scaleFloat n c) (scaleFloat n d)
isNaN !(QuadDouble a b c d) = isNaN a || isNaN b || isNaN c || isNaN d
isInfinite !(QuadDouble a b c d) = isInfinite a || isInfinite b || isInfinite c || isInfinite d
isDenormalized !(QuadDouble a b c d) = isDenormalized a || isDenormalized b || isDenormalized c || isDenormalized d
isNegativeZero !(QuadDouble a b c d) = isNegativeZero a || (a == 0 && (isNegativeZero b || (b == 0 && (isNegativeZero c || (c == 0 && isNegativeZero d)))))
isIEEE _ = False
atan2 = lift_qd_qd_qd c_qd_atan2
instance Show QuadDouble where
show = flip showFloat ""
instance Read QuadDouble where
readsPrec _ = readSigned readFloat
instance Storable QuadDouble where
sizeOf _ = 4 * sizeOf (error "Numeric.QD.QuadDouble.sizeOf" :: CDouble)
alignment _ = alignment (error "Numeric.QD.QuadDouble.alignment" :: CDouble)
peek !p = do
let !q = castPtr p
a <- peekElemOff q 0
b <- peekElemOff q 1
c <- peekElemOff q 2
d <- peekElemOff q 3
return $ QuadDouble a b c d
poke !p !(QuadDouble a b c d) = do
let !q = castPtr p
pokeElemOff q 0 a
pokeElemOff q 1 b
pokeElemOff q 2 c
pokeElemOff q 3 d
lift_qd_qd :: (Ptr CDouble -> Ptr CDouble -> IO ()) -> QuadDouble -> QuadDouble
lift_qd_qd f !a = unsafePerformIO $ with a $ \p -> alloca $ \r -> f (castPtr p) (castPtr r) >> peek r
lift_qd_qd_qd :: (Ptr CDouble -> Ptr CDouble -> Ptr CDouble -> IO ()) -> QuadDouble -> QuadDouble -> QuadDouble
lift_qd_qd_qd f !a !b = unsafePerformIO $ with a $ \p -> with b $ \q -> alloca $ \r -> f (castPtr p) (castPtr q) (castPtr r) >> peek r