module Numeric.QD.QuadDouble
( QuadDouble(QuadDouble)
, toDouble
, fromDouble
, toDoubleDouble
, fromDoubleDouble
, sqr
) where
import Foreign (Ptr, alloca, castPtr, Storable(..), with)
import System.IO.Unsafe (unsafePerformIO)
import Foreign.C (CDouble, CInt)
import Data.Bits (shiftL)
import Data.Ratio ((%))
import Data.Typeable (Typeable(..))
import Numeric (readSigned, readFloat)
import Text.FShow.Raw (BinDecode(..), DecimalFormat(..), FormatStyle(Generic), decimalFormat)
import Numeric.QD.DoubleDouble (DoubleDouble(DoubleDouble))
import Numeric.QD.QuadDouble.Raw.Safe
( 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_sqr
)
data QuadDouble = QuadDouble !CDouble !CDouble !CDouble !CDouble deriving Typeable
instance BinDecode QuadDouble where
decode (QuadDouble a b c d) =
let repack (0, 0 ) (bm, be) = (bm, be)
repack (am, ae) (0, 0 ) = (am, ae)
repack (am, ae) (bm, be) = (am `shiftL` (ae be) + bm, be)
in foldl1 repack $ map decodeFloat [a, b, c, d]
showDigits _ = 70
instance DecimalFormat QuadDouble where
nanTest (QuadDouble a _ _ _) = isNaN a
infTest (QuadDouble a _ _ _) = isInfinite a
negTest x = x < 0
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
round = round . toRational
ceiling = ceiling . toRational
floor = floor . toRational
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 Show QuadDouble where
show = decimalFormat (Generic Nothing) Nothing
instance Read QuadDouble where
readsPrec _ = readSigned readFloat
instance Storable QuadDouble where
sizeOf _ = 4 * sizeOf (0 :: CDouble)
alignment _ = alignment (0 :: 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