Internal.hs

{-# LANGUAGE CPP #-}
{-# LANGUAGE ForeignFunctionInterface #-}
{-# LANGUAGE GHCForeignImportPrim #-}
{-# LANGUAGE UnliftedFFITypes #-}

module Math.Interval.Internal.SIMD.Internal
  ( 𝕀(𝕀, MkI, PackI, inf, sup)
  , scaleInterval
  )
  where

-- base
import Prelude hiding ( Num(..), Fractional(..), Floating(..) )
import qualified Prelude
import GHC.Exts
  ( Double(D#), Double#
  , negateDouble#
  , DoubleX2#
  , packDoubleX2#, unpackDoubleX2#
  , broadcastDoubleX2#
  , plusDoubleX2#, timesDoubleX2#

  , indexDoubleX2Array#
  , readDoubleX2Array#
  , writeDoubleX2Array#
  , indexDoubleX2OffAddr#
  , readDoubleX2OffAddr#
  , writeDoubleX2OffAddr#
  )

-- primitive
import Data.Primitive

-- rounded-hw
import Numeric.Rounded.Hardware
  ( Rounded(..) )
import Numeric.Rounded.Hardware.Internal
  ( intervalFromInteger, intervalFromRational )

-- ghc-prim
import GHC.Prim
  ( maxDouble#, shuffleDoubleX2# )

-- deepseq
import Control.DeepSeq
  ( NFData(..) )

--------------------------------------------------------------------------------

-- | A non-empty interval of real numbers (possibly unbounded).
data 𝕀 = MkI DoubleX2#
  -- Interval [lo..hi] represented as (-lo, hi)

instance NFData 𝕀 where
  rnf ( MkI !_ ) = ()

foreign import prim "in2_mul8"
  mulI :: DoubleX2# -> DoubleX2# -> DoubleX2#

{-# COMPLETE PackI #-}
{-# INLINE PackI #-}
pattern PackI :: Double# -> Double# -> 𝕀
pattern PackI m_lo hi <- ( ( \ ( MkI x ) -> unpackDoubleX2# x ) -> (# m_lo, hi #) )
  where
    PackI m_lo hi = MkI ( packDoubleX2# (# m_lo, hi #) )


{-# COMPLETE 𝕀 #-}
{-# INLINE 𝕀 #-}
pattern 𝕀 :: Double -> Double -> 𝕀
pattern 𝕀 { inf, sup } <- ( ( \ ( PackI m_lo hi ) -> (# D# ( negateDouble# m_lo ), D# hi #) )
                            -> (# inf, sup #) )
  where
    𝕀 (D# lo) (D# hi) =
      let m_lo = negateDouble# lo
      in  PackI m_lo hi

instance Prelude.Num 𝕀 where
  MkI x + MkI y = MkI ( x `plusDoubleX2#` y )
  negate ( MkI x ) = MkI ( shuffleDoubleX2# x x (# 1#, 0# #) )
  ( MkI x ) * ( MkI y ) = MkI ( mulI x y )
    -- This uses the C code from Lockless Inc.
    -- TODO: compare with
    --   "Fast and Correct SIMD Algorithms for Interval Arithmetic"
    --   Frédéric Goualard, 2008
  fromInteger i =
    case intervalFromInteger i of
       ( Rounded lo, Rounded hi ) -> 𝕀 lo hi
  abs x@(PackI m_lo hi)
    | D# m_lo <= 0
    = PackI m_lo hi
    | D# hi <= 0
    = Prelude.negate x
    | otherwise
    = 𝕀 0 ( D# ( maxDouble# m_lo hi ) )
  signum _ = error "No implementation of signum for intervals"

instance Prelude.Fractional 𝕀 where
  fromRational r =
    case intervalFromRational r of
       ( Rounded lo, Rounded hi ) -> 𝕀 lo hi
  recip ( 𝕀 lo hi )
#ifdef ASSERTS
    | lo == 0
    = 𝕀 ( Prelude.recip hi ) ( 1 Prelude./ 0 )
    | hi == 0
    = 𝕀 ( -1 Prelude./ 0 ) ( Prelude.recip lo )
    | lo > 0 || hi < 0
#endif
    = 𝕀 ( Prelude.recip hi ) ( Prelude.recip lo )
#ifdef ASSERTS
    | otherwise
    = error "BAD interval recip; should use extendedRecip instead"
#endif
  p / q = p Prelude.* Prelude.recip q

scaleInterval :: Double -> 𝕀 -> 𝕀
scaleInterval s@(D# s#) i =
  case compare s 0 of
    LT -> case Prelude.negate i of { MkI x -> MkI ( broadcastDoubleX2# ( negateDouble# s# ) `timesDoubleX2#` x ) }
    EQ -> 𝕀 0 0
    GT -> case i of { MkI x -> MkI ( broadcastDoubleX2# s# `timesDoubleX2#` x ) }

--------------------------------------------------------------------------------

instance Prim 𝕀 where

  sizeOf# _ = 16#

  alignment# _ = 16#

  indexByteArray# arr# i# =
    let v = indexDoubleX2Array# arr# i#
    in MkI v

  readByteArray# arr# i# s0 =
    case readDoubleX2Array# arr# i# s0 of
      (# s1, v #) -> (# s1, MkI v #)

  writeByteArray# arr# i# (MkI v) s0 =
    writeDoubleX2Array# arr# i# v s0

  indexOffAddr# addr# i# =
    let v = indexDoubleX2OffAddr# addr# i#
    in MkI v

  readOffAddr# addr# i# s0 =
    case readDoubleX2OffAddr# addr# i# s0 of
      (# s1, v #) -> (# s1, MkI v #)

  writeOffAddr# addr# i# (MkI v) s0 =
    writeDoubleX2OffAddr# addr# i# v s0