{- |
Module      :  Numeric.QD.DoubleDouble
Copyright   :  (c) Claude Heiland-Allen 2011
License     :  BSD3

Maintainer  :  claudiusmaximus@goto10.org
Stability   :  unstable
Portability :  portable

High-level interface to libqd for double-double numbers.
-}
{-# LANGUAGE DeriveDataTypeable #-}
module Numeric.QD.DoubleDouble
  ( DoubleDouble(DoubleDouble)
  , toDouble
  , fromDouble
  , sqr
  ) where

import Foreign (Ptr, alloca, castPtr, Storable(..), unsafePerformIO, with)
import Foreign.C (CDouble, CInt)
import Data.Ratio ((%))
import Data.Bits (shiftL)
import Data.Typeable (Typeable(..))
import Numeric (readSigned, readFloat)
import Text.FShow.Raw (BinDecode(..), DecimalFormat(..), FormatStyle(Generic), decimalFormat)

import Numeric.QD.DoubleDouble.Raw.Safe
  ( 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_neg
  , c_dd_abs
  , c_dd_sqr
  )

-- | @'DoubleDouble' a b@ represents the unevaluated sum @a + b@.
data DoubleDouble = DoubleDouble {-# UNPACK #-} !CDouble {-# UNPACK #-} !CDouble deriving Typeable

instance BinDecode DoubleDouble where
  decode (DoubleDouble a b) =
    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]
  showDigits _ = 35 -- FIXME somewhat arbitrary

-- FIXME check these
instance DecimalFormat DoubleDouble where
  nanTest (DoubleDouble a _b) = isNaN a
  infTest (DoubleDouble a _b) = isInfinite a
  negTest x = x < 0

-- | Extract the first component and convert to 'Double'.
toDouble :: DoubleDouble -> Double
toDouble (DoubleDouble a _) = realToFrac a

-- | Convert from 'Double' by pairing with 0.
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 = decimalFormat (Generic Nothing) Nothing

instance Read DoubleDouble where
  readsPrec _ = readSigned readFloat

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 = fromRational (i % 1)

-- | Square a 'DoubleDouble' number.
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 a = fromRational k
                        k' = k - toRational a
                        b = fromRational k'
                    in  DoubleDouble a b

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

instance RealFrac DoubleDouble 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 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 Enum DoubleDouble -- FIXME

instance Storable DoubleDouble where
  sizeOf _ = 2 * sizeOf (0 :: CDouble)
  alignment _ = alignment (0 :: 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
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
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