Common.hs

{-# LANGUAGE BangPatterns #-}
{-# LANGUAGE DeriveFunctor #-}
{-# LANGUAGE NoImplicitPrelude #-}
module Numeric.Floating.IEEE.Internal.Rounding.Common where
import           Control.Exception (assert)
import           Data.Bits
import           Data.Functor.Product
import           Data.Int
import           GHC.Float (expt)
import           Math.NumberTheory.Logarithms (integerLog2')
import           MyPrelude
import           Numeric.Floating.IEEE.Internal.IntegerInternals

default ()

class Functor f => RoundingStrategy f where
  exact :: a -> f a
  inexact :: Ordering -- ^ LT -> toward-zero is the nearest, EQ -> midpoint, GT -> away-from-zero is the nearest
          -> Bool -- ^ negative (True -> negative, False -> positive)
          -> Int -- ^ parity (even -> toward-zero is even, odd -> toward-zero is odd)
          -> a -- ^ toward zero
          -> a -- ^ away from zero
          -> f a
  doRound :: Bool -- ^ exactness; if True, the Ordering must be LT
          -> Ordering -- ^ LT -> toward-zero is the nearest, EQ -> midpoint, GT -> away-from-zero is the nearest
          -> Bool -- ^ negative (True -> negative, False -> positive)
          -> Int -- ^ parity (even -> toward-zero is even, odd -> toward-zero is odd)
          -> a -- ^ toward zero
          -> a -- ^ away from zero
          -> f a
  exact x = doRound True LT False 0 x x
  inexact o neg parity zero away = doRound False o neg parity zero away

newtype RoundTiesToEven a = RoundTiesToEven { roundTiesToEven :: a }
  deriving (Functor)

instance RoundingStrategy RoundTiesToEven where
  exact = RoundTiesToEven
  inexact o _neg parity zero away = RoundTiesToEven $ case o of
                                                        LT -> zero
                                                        EQ | even parity -> zero
                                                           | otherwise -> away
                                                        GT -> away
  doRound _ex o _neg parity zero away = RoundTiesToEven $ case o of
    LT -> zero
    EQ | even parity -> zero
       | otherwise -> away
    GT -> away
  {-# INLINE exact #-}
  {-# INLINE inexact #-}
  {-# INLINE doRound #-}

newtype RoundTiesToAway a = RoundTiesToAway { roundTiesToAway :: a }
  deriving (Functor)

instance RoundingStrategy RoundTiesToAway where
  exact = RoundTiesToAway
  inexact o _neg _parity zero away = RoundTiesToAway $ case o of
                                                         LT -> zero
                                                         EQ -> away
                                                         GT -> away
  doRound _ex o _neg _parity zero away = RoundTiesToAway $ case o of
    LT -> zero
    EQ -> away
    GT -> away
  {-# INLINE exact #-}
  {-# INLINE inexact #-}
  {-# INLINE doRound #-}

newtype RoundTowardPositive a = RoundTowardPositive { roundTowardPositive :: a }
  deriving (Functor)

instance RoundingStrategy RoundTowardPositive where
  exact = RoundTowardPositive
  inexact _o neg _parity zero away | neg = RoundTowardPositive zero
                                   | otherwise = RoundTowardPositive away
  doRound ex _o neg _parity zero away | ex || neg = RoundTowardPositive zero
                                      | otherwise = RoundTowardPositive away
  {-# INLINE exact #-}
  {-# INLINE inexact #-}
  {-# INLINE doRound #-}

newtype RoundTowardNegative a = RoundTowardNegative { roundTowardNegative :: a }
  deriving (Functor)

instance RoundingStrategy RoundTowardNegative where
  exact = RoundTowardNegative
  inexact _o neg _parity zero away | neg = RoundTowardNegative away
                                   | otherwise = RoundTowardNegative zero
  doRound ex _o neg _parity zero away | not ex && neg = RoundTowardNegative away
                                      | otherwise = RoundTowardNegative zero
  {-# INLINE exact #-}
  {-# INLINE inexact #-}
  {-# INLINE doRound #-}

newtype RoundTowardZero a = RoundTowardZero { roundTowardZero :: a }
  deriving (Functor)

instance RoundingStrategy RoundTowardZero where
  exact = RoundTowardZero
  inexact _o _neg _parity zero _away = RoundTowardZero zero
  doRound _ex _o _neg _parity zero _away = RoundTowardZero zero
  {-# INLINE exact #-}
  {-# INLINE inexact #-}
  {-# INLINE doRound #-}

instance (RoundingStrategy f, RoundingStrategy g) => RoundingStrategy (Product f g) where
  exact x = Pair (exact x) (exact x)
  inexact o neg parity zero away = Pair (inexact o neg parity zero away) (inexact o neg parity zero away)
  doRound ex o neg parity zero away = Pair (doRound ex o neg parity zero away) (doRound ex o neg parity zero away)
  {-# INLINE exact #-}
  {-# INLINE inexact #-}
  {-# INLINE doRound #-}

{-
from GHC.Float:
expt :: Integer -> Int -> Integer
expt base n = base ^ n
-}

quotRemByExpt :: Integer -- ^ the dividend @x@
              -> Integer -- ^ base
              -> Int -- ^ the exponent @e@ (must be non-negative)
              -> (Integer, Integer) -- ^ @x \`'quotRem'\` (base ^ e)@
quotRemByExpt x 2 n    = assert (n >= 0) (x `unsafeShiftRInteger` n, x .&. (bit n - 1))
quotRemByExpt x base n = x `quotRem` expt base n
{-# INLINE quotRemByExpt #-}

multiplyByExpt :: Integer -- ^ the multiplicand @x@
               -> Integer -- ^ base
               -> Int -- ^ the exponent @e@ (must be non-negative)
               -> Integer -- ^ @x * base ^ e@
multiplyByExpt x 2 n    = assert (n >= 0) (x `unsafeShiftLInteger` n)
multiplyByExpt x base n = x * expt base n
{-# INLINE multiplyByExpt #-}

isDivisibleByExpt :: Integer -- ^ the dividend @x@
                  -> Integer -- ^ the base
                  -> Int -- ^ the exponent @e@ (must be non-negative)
                  -> Integer -- ^ the remainder @r@ (must be @x \`'rem'\` (base ^ e)@)
                  -> Bool -- ^ @r == 0@
isDivisibleByExpt x 2 e r = assert (r == x `rem` (2 ^ e)) $ x == 0 || Numeric.Floating.IEEE.Internal.IntegerInternals.countTrailingZerosInteger x >= e
isDivisibleByExpt x base e r = assert (r == x `rem` (base ^ e)) (r == 0)
{-# INLINE isDivisibleByExpt #-}

-- |
-- Assumption: @n >= 0@, @e >= 0@, and @r == n \`'rem'\` base^(e+1)@
--
-- Returns @compare r (base^e)@.
compareWithExpt :: Integer -- ^ base
                -> Integer -- ^ the number @n@ (must be non-negative)
                -> Integer -- ^ the remainder @r@ (must be @n \`'rem'\' base^(e+1)@)
                -> Int -- ^ the exponent @e@ (must be non-negative)
                -> Ordering
compareWithExpt 2 n r e = assert (r == n `rem` expt 2 (e+1)) $
  if n == 0 || integerLog2' n < e then
    -- If integerLog2 n < e (i.e. n < 2^e), it is trivial
    LT
  else
    -- In this branch, n > 0 && integerLog2' n >= e
    let result = Numeric.Floating.IEEE.Internal.IntegerInternals.roundingMode n e
        !_ = assert (result == compare r (expt 2 e)) ()
    in result
compareWithExpt base n r e = assert (r == n `rem` expt base (e+1)) $ compare r (expt base e)
{-# INLINE compareWithExpt #-}