{-# LANGUAGE DeriveDataTypeable #-}
module Numeric.QD.QuadDouble
  ( QuadDouble(QuadDouble)
  , toDouble
  , fromDouble
  , toDoubleDouble
  , fromDoubleDouble
  ) where

import Foreign (Ptr, alloca, allocaArray, castPtr, Storable(..), unsafePerformIO, with)
import Foreign.C (CDouble, CInt, peekCString)
import Data.Ratio (denominator, numerator)
import Data.Typeable (Typeable(..))

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_swrite
  , c_qd_neg
  , c_qd_abs
  , c_qd_copy_d
  , c_qd_ceil
  , c_qd_floor
  , c_qd_atan2
  )

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 Show QuadDouble where
  show a =  unsafePerformIO $ with a $ \p -> allocaArray (fromIntegral l) $ \r -> do
              c_qd_swrite (castPtr p) (l`div`2) (castPtr r) l
              peekCString r
            where l = 128
  
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 = unsafePerformIO $ alloca $ \r -> c_qd_copy_d (fromInteger i) (castPtr r) >> peek r

instance Fractional QuadDouble where
  (/) = lift_qd_qd_qd c_qd_div
  recip b = 1 / b
  fromRational k = let { n = numerator k ; d = denominator k } in fromInteger n / fromInteger d

instance Real QuadDouble where
  toRational (QuadDouble a b c d) = toRational a + toRational b + toRational c + toRational d

instance RealFrac QuadDouble where
  properFraction = error "Numeric.QD.QuadDouble.properFraction: not yet implemented" -- FIXME
  truncate = error "Numeric.QD.QuadDouble.truncate: not yet implemented" -- FIXME
  round = error "Numeric.QD.QuadDouble.round: not yet implemented" -- FIXME
  ceiling = ceiling . toDouble . lift_qd_qd c_qd_ceil
  floor = floor . toDouble . 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 = error "Numeric.QD.QuadDouble.decodeFloat: not yet implemented" -- FIXME
  encodeFloat = error "Numeric.QD.QuadDouble.encodeFloat: not yet implemented" -- FIXME
  exponent _ = error "Numeric.QD.QuadDouble.exponent: not yet implemented" -- FIXME
  significand _ = error "Numeric.QD.QuadDouble.significand: not yet implemented" -- FIXME
  scaleFloat n x = x * 2 ** fromIntegral n
  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 -- FIXME what does this imply?
  atan2 = lift_qd_qd_qd c_qd_atan2

-- instance Enum QuadDouble -- FIXME
-- instance Read QuadDouble -- FIXME

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