QuadDouble.hs

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

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

High-level interface to libqd for quad-double numbers.
-}
{-# LANGUAGE DeriveDataTypeable #-}
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)
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
  )

-- | @'QuadDouble' a b c d@ represents the unevaluated sum @a + b + c + d@.
data QuadDouble = QuadDouble {-# UNPACK #-} !CDouble {-# UNPACK #-} !CDouble {-# UNPACK #-} !CDouble {-# UNPACK #-} !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 -- FIXME somewhat arbitrary

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

-- | Convert to 'Double'.
toDouble :: QuadDouble -> Double
toDouble (QuadDouble a _ _ _) = realToFrac a

-- | Convert from 'Double'.
fromDouble :: Double -> QuadDouble
fromDouble a = QuadDouble (realToFrac a) 0 0 0

-- | Convert to 'DoubleDouble'.
toDoubleDouble :: QuadDouble -> DoubleDouble
toDoubleDouble (QuadDouble a b _ _) = DoubleDouble a b

-- | Convert from 'DoubleDouble'.
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)

-- | Square a 'QuadDouble' number.
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 Enum QuadDouble -- FIXME

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