{-# LANGUAGE BangPatterns, DeriveDataTypeable #-}
module Numeric.QD.DoubleDouble
  ( DoubleDouble(DoubleDouble)
  , toDouble
  , fromDouble
  , sqr
  ) 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.Raw
  ( c_dd_add
  , c_dd_sub
  , c_dd_mul
  , c_dd_div
  , c_dd_pi
  , c_dd_exp
  , c_dd_sqrt
  , c_dd_log
  , c_dd_sin
  , c_dd_cos
  , c_dd_tan
  , c_dd_asin
  , c_dd_acos
  , c_dd_atan
  , c_dd_sinh
  , c_dd_cosh
  , c_dd_tanh
  , c_dd_asinh
  , c_dd_acosh
  , c_dd_atanh
  , c_dd_comp
  , c_dd_swrite
  , c_dd_neg
  , c_dd_abs
  , c_dd_copy_d
  , c_dd_ceil
  , c_dd_floor
  , c_dd_atan2
  , c_dd_sqr
  )

data DoubleDouble = DoubleDouble {-# UNPACK #-} !CDouble {-# UNPACK #-} !CDouble deriving Typeable

toDouble :: DoubleDouble -> Double
toDouble !(DoubleDouble a _) = realToFrac a

fromDouble :: Double -> DoubleDouble
fromDouble !a = DoubleDouble (realToFrac a) 0

instance Eq DoubleDouble where
  (!a) == (!b) = a `compare` b == EQ
  (!a) /= (!b) = a `compare` b /= EQ

instance Ord DoubleDouble where
  (!a) `compare` (!b) = unsafePerformIO $ with a $ \p -> with b $ \q -> alloca $ \r -> do
                          c_dd_comp (castPtr p) (castPtr q) (castPtr r)
                          !i <- peek r
                          return $ i `compare` (0 :: CInt)

instance Show DoubleDouble where
  show !a = unsafePerformIO $ with a $ \p -> allocaArray (fromIntegral l) $ \r -> do
              c_dd_swrite (castPtr p) (l`div`2) (castPtr r) l
              peekCString r
            where l = 64
  
instance Num DoubleDouble where
  (+) = lift_dd_dd_dd c_dd_add
  (*) = lift_dd_dd_dd c_dd_mul
  (-) = lift_dd_dd_dd c_dd_sub
  negate = lift_dd_dd c_dd_neg
  abs = lift_dd_dd c_dd_abs
  signum !a = case a `compare` 0 of { LT -> -1 ; EQ -> 0 ; GT -> 1 }
  fromInteger !i = unsafePerformIO $ alloca $ \r -> c_dd_copy_d (fromInteger i) (castPtr r) >> peek r

sqr :: DoubleDouble -> DoubleDouble
sqr = lift_dd_dd c_dd_sqr

instance Fractional DoubleDouble where
  (/) = lift_dd_dd_dd c_dd_div
  recip !b = 1 / b
  fromRational !k = let { !n = numerator k ; !d = denominator k } in fromInteger n / fromInteger d

instance Real DoubleDouble where
  toRational !(DoubleDouble a b) = toRational a + toRational b

instance RealFrac DoubleDouble where
  properFraction = error "Numeric.QD.DoubleDouble.properFraction: not yet implemented" -- FIXME
  truncate = error "Numeric.QD.DoubleDouble.truncate: not yet implemented" -- FIXME
  round = error "Numeric.QD.DoubleDouble.round: not yet implemented" -- FIXME
  ceiling = ceiling . toDouble . lift_dd_dd c_dd_ceil
  floor = floor . toDouble . lift_dd_dd c_dd_floor

instance Floating DoubleDouble where
  pi = unsafePerformIO $ alloca $ \r -> c_dd_pi (castPtr r) >> peek r
  exp = lift_dd_dd c_dd_exp
  sqrt = lift_dd_dd c_dd_sqrt
  log = lift_dd_dd c_dd_log
  sin = lift_dd_dd c_dd_sin
  cos = lift_dd_dd c_dd_cos
  tan = lift_dd_dd c_dd_tan
  asin = lift_dd_dd c_dd_asin
  acos = lift_dd_dd c_dd_acos
  atan = lift_dd_dd c_dd_atan
  sinh = lift_dd_dd c_dd_sinh
  cosh = lift_dd_dd c_dd_cosh
  tanh = lift_dd_dd c_dd_tanh
  asinh = lift_dd_dd c_dd_asinh
  acosh = lift_dd_dd c_dd_acosh
  atanh = lift_dd_dd c_dd_atanh

instance RealFloat DoubleDouble where
  floatRadix _ = 2
  floatDigits _ = 2 * floatDigits (error "Numeric.QD.DoubleDouble.floatDigits" :: CDouble)
  floatRange _ = floatRange (error "Numeric.QD.DoubleDouble.floatRange" :: CDouble)
  decodeFloat = error "Numeric.QD.DoubleDouble.decodeFloat: not yet implemented" -- FIXME
  encodeFloat = error "Numeric.QD.DoubleDouble.encodeFloat: not yet implemented" -- FIXME
  exponent _ = error "Numeric.QD.DoubleDouble.exponent: not yet implemented" -- FIXME
  significand _ = error "Numeric.QD.DoubleDouble.significand: not yet implemented" -- FIXME
  scaleFloat !n !x = x * 2 ** fromIntegral n
  isNaN !(DoubleDouble a b) = isNaN a || isNaN b
  isInfinite !(DoubleDouble a b) = isInfinite a || isInfinite b
  isDenormalized !(DoubleDouble a b) = isDenormalized a || isDenormalized b
  isNegativeZero !(DoubleDouble a b) = isNegativeZero a || (a == 0 && isNegativeZero b)
  isIEEE _ = False -- FIXME what does this imply?
  atan2 = lift_dd_dd_dd c_dd_atan2

-- instance Enum DoubleDouble -- FIXME
-- instance Read DoubleDouble -- FIXME

instance Storable DoubleDouble where
  sizeOf _ = 2 * sizeOf (error "Numeric.QD.DoubleDouble.sizeOf" :: CDouble)
  alignment _ = alignment (error "Numeric.QD.DoubleDouble.alignment" :: CDouble)
  peek !p = do
    let !q = castPtr p
    a <- peekElemOff q 0
    b <- peekElemOff q 1
    return $ DoubleDouble a b
  poke !p !(DoubleDouble a b) = do
    let !q = castPtr p
    pokeElemOff q 0 a
    pokeElemOff q 1 b

lift_dd_dd :: (Ptr CDouble -> Ptr CDouble -> IO ()) -> DoubleDouble -> DoubleDouble
{-# INLINE lift_dd_dd #-}
lift_dd_dd f !a = unsafePerformIO $ with a $ \p -> alloca $ \r -> f (castPtr p) (castPtr r) >> peek r

lift_dd_dd_dd :: (Ptr CDouble -> Ptr CDouble -> Ptr CDouble -> IO ()) -> DoubleDouble -> DoubleDouble -> DoubleDouble
{-# INLINE lift_dd_dd_dd #-}
lift_dd_dd_dd f !a !b = unsafePerformIO $ with a $ \p -> with b $ \q -> alloca $ \r -> f (castPtr p) (castPtr q) (castPtr r) >> peek r