packages feed

fp-ieee-0.1.0.6: src/Numeric/Floating/IEEE/Internal/Float128.hs

{-# LANGUAGE CPP #-}
{-# LANGUAGE MagicHash #-}
{-# LANGUAGE NoImplicitPrelude #-}
{-# LANGUAGE NumericUnderscores #-}
{-# OPTIONS_GHC -Wno-orphans -Wno-unused-imports #-}
module Numeric.Floating.IEEE.Internal.Float128 where
import           Data.Bits
import           Data.Word
import           GHC.Exts (Int#)
import           MyPrelude
import           Numeric.Float128 (Float128 (F128))
import qualified Numeric.Float128
import           Numeric.Floating.IEEE.Internal.Base
import           Numeric.Floating.IEEE.Internal.Classify
import           Numeric.Floating.IEEE.Internal.Conversion
import           Numeric.Floating.IEEE.Internal.FMA
import           Numeric.Floating.IEEE.Internal.NaN (RealFloatNaN)
import qualified Numeric.Floating.IEEE.Internal.NaN as NaN
import           Numeric.Floating.IEEE.Internal.NextFloat
import           Numeric.Floating.IEEE.Internal.Rounding
import           Numeric.Floating.IEEE.Internal.RoundToIntegral

default ()

{-
Float128:
- exponent = 15 bits
- precision = 113 bits
-}

float128ToWord64Hi, float128ToWord64Lo :: Float128 -> Word64
float128ToWord64Hi (F128 hi _lo) = hi
float128ToWord64Lo (F128 _hi lo) = lo
{-# INLINE float128ToWord64Hi #-}
{-# INLINE float128ToWord64Lo #-}

float128ToWord64Pair :: Float128 -> (Word64, Word64)
float128ToWord64Pair (F128 hi lo) = (hi, lo)
{-# INLINE float128ToWord64Pair #-}

float128FromWord64Pair :: Word64 -- ^ higher 64 bits
                       -> Word64 -- ^ lower 64 bits
                       -> Float128
float128FromWord64Pair hi lo = F128 hi lo
{-# INLINE float128FromWord64Pair #-}

succWord64Pair :: Word64 -> Word64 -> (Word64, Word64)
succWord64Pair hi lo | lo + 1 == 0 = (hi + 1, 0)
                     | otherwise = (hi, lo + 1)

predWord64Pair :: Word64 -> Word64 -> (Word64, Word64)
predWord64Pair hi lo | lo == 0 = (hi - 1, fromInteger (-1))
                     | otherwise = (hi, lo - 1)

nextUpF128 :: Float128 -> Float128
nextUpF128 x =
  case float128ToWord64Pair x of
    (hi, lo) | hi .&. 0x7fff_0000_0000_0000 == 0x7fff_0000_0000_000
             , (hi, lo) /= (0xffff_0000_0000_0000, 0) -> x + x -- NaN or positive infinity -> itself
    (0x8000_0000_0000_0000, 0x0000_0000_0000_0000) -> minPositive -- -0 -> min positive
    (hi, lo) | testBit hi 63 -> -- negative
                 case predWord64Pair hi lo of
                   (hi', lo') -> float128FromWord64Pair hi' lo'
             | otherwise -> -- positive
                 case succWord64Pair hi lo of
                   (hi', lo') -> float128FromWord64Pair hi' lo'

nextDownF128 :: Float128 -> Float128
nextDownF128 x =
  case float128ToWord64Pair x of
    (hi, lo) | hi .&. 0x7fff_0000_0000_0000 == 0x7fff_0000_0000_000
             , (hi, lo) /= (0x7fff_0000_0000_0000, 0) -> x + x -- NaN or negative infinity -> itself
    (0x0000_0000_0000_0000, 0x0000_0000_0000_0000) -> - minPositive -- +0 -> max negative
    (hi, lo) | testBit hi 63 -> -- negative
                 case succWord64Pair hi lo of
                   (hi', lo') -> float128FromWord64Pair hi' lo'
             | otherwise -> -- positive
                 case predWord64Pair hi lo of
                   (hi', lo') -> float128FromWord64Pair hi' lo'

nextTowardZeroF128 :: Float128 -> Float128
nextTowardZeroF128 x =
  case float128ToWord64Pair x of
    (hi, lo) | hi .&. 0x7fff_0000_0000_0000 == 0x7fff_0000_0000_000
             , (lo, hi .&. 0x0000_ffff_ffff_ffff) /= (0, 0) -> x + x -- NaN -> itself
    (0x8000_0000_0000_0000, 0x0000_0000_0000_0000) -> x -- -0 -> itself
    (0x0000_0000_0000_0000, 0x0000_0000_0000_0000) -> x -- +0 -> itself
    (hi, lo) -> -- positive / negative
      case predWord64Pair hi lo of
        (hi', lo') -> float128FromWord64Pair hi' lo'

isNormalF128 :: Float128 -> Bool
isNormalF128 x = case float128ToWord64Pair x of
                   (hi, _) -> let hi' = hi .&. 0x7fff_0000_0000_0000
                              in hi' /= 0 && hi' /= 0x7fff_0000_0000_0000

isFiniteF128 :: Float128 -> Bool
isFiniteF128 x = case float128ToWord64Pair x of
                   (hi, _) -> let hi' = hi .&. 0x7fff_0000_0000_0000
                              in hi' /= 0 && hi' /= 0x7fff_0000_0000_0000

classifyF128DiscardingSignalingNaNs :: Float128 -> Class
classifyF128DiscardingSignalingNaNs x =
  let hi = float128ToWord64Hi x
      s = testBit hi 63
      e = (hi `unsafeShiftR` 48) .&. 0x7fff -- exponent (15 bits)
      m_hi = hi .&. 0x0000_ffff_ffff_ffff
      m_lo = float128ToWord64Lo x
  in case (s, e, m_hi, m_lo) of
       (True,  0,      0, 0) -> NegativeZero
       (False, 0,      0, 0) -> PositiveZero
       (True,  0,      _, _) -> NegativeSubnormal
       (False, 0,      _, _) -> PositiveSubnormal
       (True,  0x7fff, 0, 0) -> NegativeInfinity
       (False, 0x7fff, 0, 0) -> PositiveInfinity
       (_,     0x7fff, _, _) -> QuietNaN -- treat all NaNs as quiet
       (True,  _,      _, _) -> NegativeNormal
       (False, _,      _, _) -> PositiveNormal

instance RealFloatNaN Float128 where
  copySign x y = let (x_hi, x_lo) = float128ToWord64Pair x
                     y_hi = float128ToWord64Hi y
                 in float128FromWord64Pair ((x_hi .&. 0x7fff_ffff_ffff_ffff) .|. (y_hi .&. 0x8000_0000_0000_0000)) x_lo
  isSignMinus x = let hi = float128ToWord64Hi x
                  in testBit hi 63
  isSignaling x = let hi = float128ToWord64Hi x
                  in isNaN x && not (testBit hi 47)

  getPayload x
    | not (isNaN x) = -1
    | otherwise = let hi = fromIntegral (float128ToWord64Hi x .&. 0x0000_7fff_ffff_ffff)
                      lo = fromIntegral (float128ToWord64Lo x)
                  in hi * 0x1_0000_0000_0000_0000 + lo

  setPayload x
    | 0 <= x && x <= 0x0000_7fff_ffff_ffff_ffff_ffff_ffff_ffff
    = let payloadI = round x
          hi = fromInteger (payloadI `shiftR` 64) .|. 0x7fff_8000_0000_0000
          lo = fromInteger (payloadI .&. 0xffff_ffff_ffff_ffff)
      in float128FromWord64Pair hi lo
    | otherwise = 0

  setPayloadSignaling x
    | 0 < x && x <= 0x0000_7fff_ffff_ffff_ffff_ffff_ffff_ffff
    = let payloadI = round x
          hi = fromInteger (payloadI `shiftR` 64) .|. 0x7fff_0000_0000_0000
          lo = fromInteger (payloadI .&. 0xffff_ffff_ffff_ffff)
      in float128FromWord64Pair hi lo
    | otherwise = 0

  classify x =
    let hi = float128ToWord64Hi x
        s = testBit hi 63
        e = (hi `unsafeShiftR` 48) .&. 0x7fff -- exponent (15 bits)
        m_hi = hi .&. 0x0000_ffff_ffff_ffff
        m_lo = float128ToWord64Lo x
    in case (s, e, m_hi, m_lo) of
         (True,  0,      0, 0) -> NegativeZero
         (False, 0,      0, 0) -> PositiveZero
         (True,  0,      _, _) -> NegativeSubnormal
         (False, 0,      _, _) -> PositiveSubnormal
         (True,  0x7fff, 0, 0) -> NegativeInfinity
         (False, 0x7fff, 0, 0) -> PositiveInfinity
         (_,     0x7fff, _, _) -> if testBit m_hi 47 then
                                    QuietNaN
                                  else
                                    SignalingNaN
         (True,  _,      _, _) -> NegativeNormal
         (False, _,      _, _) -> PositiveNormal

  compareByTotalOrder x y =
    let (x_hi, x_lo) = float128ToWord64Pair x
        (y_hi, y_lo) = float128ToWord64Pair y
    in compare (testBit y_hi 63) (testBit x_hi 63) -- sign bit
       <> if testBit x_hi 63 then
            compare y_hi x_hi <> compare y_lo x_lo -- negative
          else
            compare x_hi y_hi <> compare x_lo y_lo -- positive

{-# RULES
"nextUp/Float128" nextUp = nextUpF128
"nextDown/Float128" nextDown = nextDownF128
"nextTowardZero/Float128" nextTowardZero = nextTowardZeroF128
"isNormal/F128" isNormal = isNormalF128
"isFinite/F128" isFinite = isFiniteF128
"classify/F128" classify = classifyF128DiscardingSignalingNaNs
"isMantissaEven/F128"
  isMantissaEven = \x -> case x :: Float128 of F128 _hi lo -> even lo
"roundAway'/Float128" roundAway' = Numeric.Float128.round'
"ceiling'/Float128" ceiling' = Numeric.Float128.ceiling'
"floor'/Float128" floor' = Numeric.Float128.floor'
"truncate'/Float128" truncate' = Numeric.Float128.truncate'
  #-}

-- TODO: Write directly?
{-# SPECIALIZE minPositive :: Float128 #-}
{-# SPECIALIZE minPositiveNormal :: Float128 #-}
{-# SPECIALIZE maxFinite :: Float128 #-}

-- We shouldn't need specializations of positiveWordToBinaryFloatR# as long as WORD_SIZE_IN_BITS <= 113
{-# SPECIALIZE fromPositiveIntegerR :: RoundingStrategy f => Bool -> Integer -> f Float128 #-}
{-# SPECIALIZE fromPositiveIntegerR :: Bool -> Integer -> RoundTiesToEven Float128 #-}
{-# SPECIALIZE fromPositiveIntegerR :: Bool -> Integer -> RoundTiesToAway Float128 #-}
{-# SPECIALIZE fromPositiveIntegerR :: Bool -> Integer -> RoundTowardPositive Float128 #-}
{-# SPECIALIZE fromPositiveIntegerR :: Bool -> Integer -> RoundTowardNegative Float128 #-}
{-# SPECIALIZE fromPositiveIntegerR :: Bool -> Integer -> RoundTowardZero Float128 #-}

{-# SPECIALIZE fromPositiveRatioR :: RoundingStrategy f => Bool -> Integer -> Integer -> f Float128 #-}
{-# SPECIALIZE fromPositiveRatioR :: Bool -> Integer -> Integer -> RoundTiesToEven Float128 #-}
{-# SPECIALIZE fromPositiveRatioR :: Bool -> Integer -> Integer -> RoundTiesToAway Float128 #-}
{-# SPECIALIZE fromPositiveRatioR :: Bool -> Integer -> Integer -> RoundTowardPositive Float128 #-}
{-# SPECIALIZE fromPositiveRatioR :: Bool -> Integer -> Integer -> RoundTowardNegative Float128 #-}
{-# SPECIALIZE fromPositiveRatioR :: Bool -> Integer -> Integer -> RoundTowardZero Float128 #-}

{-# SPECIALIZE encodePositiveFloatR# :: RoundingStrategy f => Bool -> Integer -> Int# -> f Float128 #-}
{-# SPECIALIZE encodePositiveFloatR# :: Bool -> Integer -> Int# -> RoundTiesToEven Float128 #-}
{-# SPECIALIZE encodePositiveFloatR# :: Bool -> Integer -> Int# -> RoundTiesToAway Float128 #-}
{-# SPECIALIZE encodePositiveFloatR# :: Bool -> Integer -> Int# -> RoundTowardPositive Float128 #-}
{-# SPECIALIZE encodePositiveFloatR# :: Bool -> Integer -> Int# -> RoundTowardNegative Float128 #-}
{-# SPECIALIZE encodePositiveFloatR# :: Bool -> Integer -> Int# -> RoundTowardZero Float128 #-}

{-# SPECIALIZE scaleFloatR# :: RoundingStrategy f => Int# -> Float128 -> f Float128 #-}
{-# SPECIALIZE scaleFloatR# :: Int# -> Float128 -> RoundTiesToEven Float128 #-}
{-# SPECIALIZE scaleFloatR# :: Int# -> Float128 -> RoundTiesToAway Float128 #-}
{-# SPECIALIZE scaleFloatR# :: Int# -> Float128 -> RoundTowardPositive Float128 #-}
{-# SPECIALIZE scaleFloatR# :: Int# -> Float128 -> RoundTowardNegative Float128 #-}
{-# SPECIALIZE scaleFloatR# :: Int# -> Float128 -> RoundTowardZero Float128 #-}