SIMD.hs

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

{-# OPTIONS_GHC -Wno-orphans #-}

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

-- base
import Prelude hiding ( Num(..), Fractional(..), Floating(..) )
import qualified Prelude
import GHC.Exts
  ( (<=##), isTrue#
  , Double(D#)
  )

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

-- rounded-hw
import Numeric.Rounded.Hardware
  ( Rounded(..) )
import qualified Numeric.Rounded.Hardware.Interval.NonEmpty as Interval
  ( Interval(..) )

-- brush-strokes
import Math.Ring

-- brush-strokes:simd-interval
import Math.Interval.Internal.SIMD.Internal

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

deriving via ViaPrelude 𝕀
  instance AbelianGroup 𝕀
deriving via ViaPrelude 𝕀
  instance Field 𝕀

instance Ring 𝕀 where
  (*) = (Prelude.*)
  iv ^ 1 = iv
  PackI m_lo hi ^ n
    | odd n || isTrue# (m_lo <=## 0.0##)
    , let !(D# m_lo') = D# m_lo Prelude.^ n
          !(D# hi'  ) = D# hi   Prelude.^ n
    = PackI m_lo' hi'
    | isTrue# (hi <=## 0.0##)
    , let !(D# m_lo') = Prelude.negate $ (D# hi) Prelude.^ n
          !(D# hi')   = (D# m_lo) Prelude.^ n
    = PackI m_lo' hi'
    | otherwise
    , let !(D# hi1) = (D# m_lo) Prelude.^ n
          !(D# hi2) = (D# hi  ) Prelude.^ n
          hi' = maxDouble# hi1 hi2
    = PackI 0.0## hi'

instance Floating 𝕀 where
  sqrt = withHW Prelude.sqrt

instance Transcendental 𝕀 where
  pi = 𝕀 3.141592653589793 3.1415926535897936
  cos = withHW Prelude.cos
  sin = withHW Prelude.sin
  atan = withHW Prelude.atan

{-
TODO: consider alternatives for sin/cos, such as:

  - https://github.com/JishinMaster/simd_utils/blob/160c50f07e08d2077ae4368f0aed2f10f7173c67/simd_utils_sse_double.h#L530
  - Intel SVML
  - metalibm
  - sleef

-}

{-# INLINE withHW #-}
-- | Internal function: use @rounded-hw@ to define a function on intervals.
withHW :: (Interval.Interval Double -> Interval.Interval Double) -> 𝕀 -> 𝕀
withHW f = \ ( 𝕀 lo hi ) ->
  case f ( Interval.I ( Rounded lo ) ( Rounded hi ) ) of
    Interval.I ( Rounded x ) ( Rounded y ) -> 𝕀 x y