packages feed

brush-strokes-0.1.0.0: src/lib/Math/Interval/Internal/FMA.hs

{-# LANGUAGE AllowAmbiguousTypes #-}
{-# LANGUAGE PolyKinds #-}
{-# LANGUAGE ScopedTypeVariables #-}

module Math.Interval.Internal.FMA
  ( 𝕀(..)
  , scaleInterval
  , addI, subI, prodI, divI
  , withError1
  ) where

-- base
import Prelude hiding ( Num(..), Fractional(..), Floating(..) )
import qualified Prelude
import GHC.Exts
  ( Double(D#), fmsubDouble#, fnmaddDouble#, TYPE )

-- brush-strokes
import Math.Float.Utils
  ( prevFP, succFP )
import Math.Ring

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

-- rounded-hw
import Numeric.Rounded.Hardware
  ( Rounded(..) )
import Numeric.Rounded.Hardware.Class
  ( intervalFromRational, intervalFromInteger )
import qualified Numeric.Rounded.Hardware.Interval.NonEmpty as Interval
  ( Interval(..) )

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

{-# INLINE withError1 #-}
withError1 :: Double -> Double -> ( Double, Double )
withError1 f e =
  case compare e 0 of
    LT -> ( prevFP f, f )
    EQ -> ( f, f )
    GT -> ( f, succFP f )

addI :: forall rep (r :: TYPE rep). ( Double -> Double -> r ) -> Double -> Double -> r
addI withErr a b =
  let !s  = a + b
      !a' = s - b
      !b' = s - a'
      !da = a - a'
      !db = b - b'
      !e = da + db
  in s `withErr` e
{-# SPECIALISE addI withError1 #-}
{-# INLINEABLE addI #-}

subI :: forall rep (r :: TYPE rep). ( Double -> Double -> r ) -> Double -> Double -> r
subI withErr a b =
  let !s  = a - b
      !a' = s - b
      !b' = s - a'
      !da = a - a'
      !db = b - b'
      !e = da + db
  in s `withErr` e
{-# SPECIALISE subI withError1 #-}
{-# INLINEABLE subI #-}

{-# INLINE fmsubDouble #-}
fmsubDouble :: Double -> Double -> Double -> Double
fmsubDouble ( D# x ) ( D# y ) ( D# z ) = D# ( fmsubDouble# x y z )
{-# INLINE fnmaddDouble #-}
fnmaddDouble :: Double -> Double -> Double -> Double
fnmaddDouble ( D# x ) ( D# y ) ( D# z ) = D# ( fnmaddDouble# x y z )

prodI :: forall rep (r :: TYPE rep). ( Double -> Double -> r ) -> Double -> Double -> r
prodI withErr a b =
  let !p = a * b
      !e = fmsubDouble a b p
  in p `withErr` e
{-# SPECIALISE prodI withError1 #-}
{-# INLINEABLE prodI #-}

divI :: forall rep (r :: TYPE rep). ( Double -> Double -> r ) -> Double -> Double -> r
divI withErr a b =
  let !r = a / b
      !e = fnmaddDouble r b a
  in r `withErr` e
{-# SPECIALISE divI withError1 #-}
{-# INLINEABLE divI #-}

-- | Power of a **non-negative** number to a natural power.
posPowI :: Double -- ^ Assumed to be non-negative!
        -> Word -> ( Double, Double )
posPowI _ 0 = ( 1, 1 )
posPowI f 1 = ( f, f )
posPowI f 2 = prodI withError1 f f
posPowI f n
  | even n
  , let m = n `quot` 2
        ( fΒ²_lo, fΒ²_hi ) = prodI withError1 f f
  = ( fst $ posPowI fΒ²_lo m, snd $ posPowI fΒ²_hi m )
  | otherwise
  , let m = n `quot` 2
        ( fΒ²_lo, fΒ²_hi ) = prodI withError1 f f
  = ( fst $ posPowAcc fΒ²_lo m f, snd $ posPowAcc fΒ²_hi m f )

posPowAcc :: Double -> Word -> Double -> ( Double, Double )
posPowAcc f 1 x = prodI withError1 f x
posPowAcc f n x
  | even n
  , let m = n `quot` 2
        ( fΒ²_lo, fΒ²_hi ) = prodI withError1 f f
  = ( fst $ posPowAcc fΒ²_lo m x, snd $ posPowAcc fΒ²_hi m x )
  | otherwise
  , let m = n `quot` 2
        ( fΒ²_lo, fΒ²_hi ) = prodI withError1 f f
        ( y_lo, y_hi ) = prodI withError1 f x
  = ( fst $ posPowAcc fΒ²_lo m y_lo, snd $ posPowAcc fΒ²_hi m y_hi )

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

-- | A non-empty interval of real numbers (possibly unbounded).
data 𝕀 = 𝕀 { inf, sup :: !Double }

instance NFData 𝕀 where
  rnf ( 𝕀 lo hi ) = rnf lo `seq` rnf hi

instance Prelude.Num 𝕀 where
  𝕀 x_lo x_hi + 𝕀 y_lo y_hi
    | let !z_lo = fst $ addI withError1 x_lo y_lo
          !z_hi = snd $ addI withError1 x_hi y_hi
    = 𝕀 z_lo z_hi
  𝕀 x_lo x_hi - 𝕀 y_lo y_hi
    | let !z_lo = fst $ subI withError1 x_lo y_hi
          !z_hi = snd $ subI withError1 x_hi y_lo
    = 𝕀 z_lo z_hi
  negate (𝕀 lo hi) = 𝕀 (Prelude.negate hi) (Prelude.negate lo)
  (*) = (*)
  fromInteger i =
    case intervalFromInteger i of
       ( Rounded lo, Rounded hi ) -> 𝕀 lo hi
  abs (𝕀 lo hi)
    | 0 <= lo
    = 𝕀 lo hi
    | hi <= 0
    = 𝕀 (Prelude.negate hi) (Prelude.negate lo)
    | otherwise
    = 𝕀 0 (max (Prelude.negate lo) hi)
  signum _ = error "No implementation of signum for intervals"

instance Ring 𝕀 where
  𝕀 lo1 hi1 * 𝕀 lo2 hi2
    | let !( x_min, x_max ) = prodI withError1 lo1 lo2
          !( y_min, y_max ) = prodI withError1 lo1 hi2
          !( z_min, z_max ) = prodI withError1 hi1 lo2
          !( w_min, w_max ) = prodI withError1 hi1 hi2
    = 𝕀 ( min ( min x_min y_min ) ( min z_min w_min ) )
        ( max ( max x_max y_max ) ( max z_max w_max ) )
  _ ^ 0 = 𝕀 1 1
  iv ^ 1 = iv
  𝕀 lo hi ^ n
    | odd n || 0 <= lo
    , let !lo' = fst $ posPowI lo n
          !hi' = snd $ posPowI hi n
    = 𝕀 lo' hi'
    | hi <= 0
    , let !lo' = fst $ posPowI (negate hi) n
          !hi' = snd $ posPowI (negate lo) n
    = 𝕀 lo' hi'
    | otherwise
    , let !hi1 = snd $ posPowI (negate lo) n
    , let !hi2 = snd $ posPowI hi n
    = 𝕀 0 ( max hi1 hi2 )

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

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

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

{-# 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

scaleInterval :: Double -> 𝕀 -> 𝕀
scaleInterval s ( 𝕀 lo hi ) =
  case compare s 0 of
    LT -> 𝕀 ( fst $ prodI withError1 s hi ) ( snd $ prodI withError1 s lo )
    EQ -> 𝕀 0 0
    GT -> 𝕀 ( fst $ prodI withError1 s lo ) ( snd $ prodI withError1 s hi )