{-# LANGUAGE ForeignFunctionInterface #-}
-----------------------------------------------------------------------------
-- |
-- Module     : Numeric.IEEE
-- Copyright  : Copyright (c) 2010, Patrick Perry <patperry@gmail.com>
-- License    : BSD3
-- Maintainer : Patrick Perry <patperry@gmail.com>
-- Stability  : experimental
--
-- Operations on IEEE floating point numbers.
--
module Numeric.IEEE (
    -- * IEEE type class
    IEEE(..),
    -- * NaN-aware minimum and maximum
    minNum,
    maxNum,
    minNaN,
    maxNaN,
    ) where

import Data.Word
import Foreign.C.Types( CFloat, CDouble )

-- | IEEE floating point types.
class (RealFloat a) => IEEE a where
    -- | Infinity value.
    infinity :: a

    -- | The smallest representable positive value.
    minDenormal :: a

    -- | The smallest representable positive normalized value.
    minNormal :: a

    -- | The largest representable finite value.
    maxFinite :: a

    -- | The smallest positive value @x@ such that @1 + x@ is representable.
    epsilon :: a

    -- | @copySign x y@ returns @x@ with its sign changed to @y@'s.
    copySign :: a -> a -> a

    -- | Return 'True' if two values are /exactly/ (bitwise) equal.
    identicalIEEE :: a -> a -> Bool

    -- | Return the next largest IEEE value (@Infinity@ and @NaN@ are
    -- unchanged).
    succIEEE :: a -> a

    -- | Return the next smallest IEEE value (@-Infinity@ and @NaN@ are
    -- unchanged).
    predIEEE :: a -> a

    -- | Given two values with the same sign, return the value halfway
    -- between them on the IEEE number line.  If the signs of the values
    -- differ or either is @NaN@, the value is undefined.
    bisectIEEE :: a -> a -> a

    -- | The number of significand bits which are equal in the two arguments
    -- (equivalent to @feqrel@ from the Tango Math library).  The result is
    -- between @0@ and @'floatDigits'@.
    sameSignificandBits :: a -> a -> Int

    -- | Default @NaN@ value.
    nan :: a

    -- | Quiet @NaN@ value with a positive integer payload.  Payload must be
    -- less than 'maxNaNPayload'.  Beware that while some platforms allow
    -- using @0@ as a payload, this behavior is not portable.
    nanWithPayload :: Word64 -> a

    -- | Maximum @NaN@ payload for type @a@.
    maxNaNPayload :: a -> Word64

    -- | The payload stored in a @NaN@ value.  Undefined if the argument
    -- is not @NaN@.
    nanPayload :: a -> Word64


-- | Return the maximum of two values; if one value is @NaN@, return the
-- other.  Prefer the first if both values are @NaN@.
maxNum :: (RealFloat a) => a -> a -> a
maxNum :: forall a. RealFloat a => a -> a -> a
maxNum a
x a
y | a
x a -> a -> Bool
forall a. Ord a => a -> a -> Bool
>= a
y Bool -> Bool -> Bool
|| a -> Bool
forall a. RealFloat a => a -> Bool
isNaN a
y = a
x
           | Bool
otherwise         = a
y
{-# INLINE maxNum #-}

-- | Return the minimum of two values; if one value is @NaN@, return the
-- other.  Prefer the first if both values are @NaN@.
minNum :: (RealFloat a) => a -> a -> a
minNum :: forall a. RealFloat a => a -> a -> a
minNum a
x a
y | a
x a -> a -> Bool
forall a. Ord a => a -> a -> Bool
<= a
y Bool -> Bool -> Bool
|| a -> Bool
forall a. RealFloat a => a -> Bool
isNaN a
y = a
x
           | Bool
otherwise         = a
y
{-# INLINE minNum #-}

-- | Return the maximum of two values; if one value is @NaN@, return it.
-- Prefer the first if both values are @NaN@.
maxNaN :: (RealFloat a) => a -> a -> a
maxNaN :: forall a. RealFloat a => a -> a -> a
maxNaN a
x a
y | a
x a -> a -> Bool
forall a. Ord a => a -> a -> Bool
>= a
y Bool -> Bool -> Bool
|| a -> Bool
forall a. RealFloat a => a -> Bool
isNaN a
x = a
x
           | Bool
otherwise         = a
y
{-# INLINE maxNaN #-}

-- | Return the minimum of two values; if one value is @NaN@, return it.
-- Prefer the first if both values are @NaN@.
minNaN :: (RealFloat a) =>  a -> a -> a
minNaN :: forall a. RealFloat a => a -> a -> a
minNaN a
x a
y | a
x a -> a -> Bool
forall a. Ord a => a -> a -> Bool
<= a
y Bool -> Bool -> Bool
|| a -> Bool
forall a. RealFloat a => a -> Bool
isNaN a
x = a
x
           | Bool
otherwise         = a
y
{-# INLINE minNaN #-}


instance IEEE Float where
    identicalIEEE :: Float -> Float -> Bool
identicalIEEE Float
x Float
y = Float -> Float -> Int
c_identicalf Float
x Float
y Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
/= Int
0
    {-# INLINE identicalIEEE #-}
    infinity :: Float
infinity = Float
1Float -> Float -> Float
forall a. Fractional a => a -> a -> a
/Float
0
    {-# INLINE infinity #-}
    nan :: Float
nan = (Float
0Float -> Float -> Float
forall a. Fractional a => a -> a -> a
/Float
0)
    {-# INLINE nan #-}
    nanWithPayload :: Word64 -> Float
nanWithPayload Word64
n = Word32 -> Float
c_mknanf (Word64 -> Word32
forall a b. (Integral a, Num b) => a -> b
fromIntegral Word64
n)
    {-# INLINE nanWithPayload #-}
    maxNaNPayload :: Float -> Word64
maxNaNPayload Float
_ = Word64
0x003FFFFF
    {-# INLINE maxNaNPayload #-}
    nanPayload :: Float -> Word64
nanPayload Float
x = Word32 -> Word64
forall a b. (Integral a, Num b) => a -> b
fromIntegral (Word32 -> Word64) -> Word32 -> Word64
forall a b. (a -> b) -> a -> b
$ Float -> Word32
c_getnanf Float
x
    {-# INLINE nanPayload #-}
    minDenormal :: Float
minDenormal = Float
1e-45
    {-# INLINE minDenormal #-}
    minNormal :: Float
minNormal = Float
1.17549435e-38
    {-# INLINE minNormal #-}
    maxFinite :: Float
maxFinite = Float
3.40282347e+38
    {-# INLINE maxFinite #-}
    epsilon :: Float
epsilon = Float
1.19209290e-07
    {-# INLINE epsilon #-}
    copySign :: Float -> Float -> Float
copySign = Float -> Float -> Float
c_copysignf
    {-# INLINE copySign #-}
    succIEEE :: Float -> Float
succIEEE = Float -> Float
c_ieeesuccf
    {-# INLINE succIEEE #-}
    predIEEE :: Float -> Float
predIEEE = Float -> Float
c_ieeepredf
    {-# INLINE predIEEE #-}
    bisectIEEE :: Float -> Float -> Float
bisectIEEE = Float -> Float -> Float
c_ieeemeanf
    {-# INLINE bisectIEEE #-}
    sameSignificandBits :: Float -> Float -> Int
sameSignificandBits = Float -> Float -> Int
c_feqrelf
    {-# INLINE sameSignificandBits #-}


instance IEEE CFloat where
    identicalIEEE :: CFloat -> CFloat -> Bool
identicalIEEE CFloat
x CFloat
y = Float -> Float -> Int
c_identicalf (CFloat -> Float
forall a b. (Real a, Fractional b) => a -> b
realToFrac CFloat
x) (CFloat -> Float
forall a b. (Real a, Fractional b) => a -> b
realToFrac CFloat
y) Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
/= Int
0
    {-# INLINE identicalIEEE #-}
    infinity :: CFloat
infinity = CFloat
1CFloat -> CFloat -> CFloat
forall a. Fractional a => a -> a -> a
/CFloat
0
    {-# INLINE infinity #-}
    nan :: CFloat
nan = (CFloat
0CFloat -> CFloat -> CFloat
forall a. Fractional a => a -> a -> a
/CFloat
0)
    {-# INLINE nan #-}
    nanWithPayload :: Word64 -> CFloat
nanWithPayload Word64
n = Float -> CFloat
forall a b. (Real a, Fractional b) => a -> b
realToFrac (Float -> CFloat) -> Float -> CFloat
forall a b. (a -> b) -> a -> b
$ Word32 -> Float
c_mknanf (Word64 -> Word32
forall a b. (Integral a, Num b) => a -> b
fromIntegral Word64
n)
    {-# INLINE nanWithPayload #-}
    maxNaNPayload :: CFloat -> Word64
maxNaNPayload CFloat
_ = Word64
0x003FFFFF
    {-# INLINE maxNaNPayload #-}
    nanPayload :: CFloat -> Word64
nanPayload CFloat
x = Word32 -> Word64
forall a b. (Integral a, Num b) => a -> b
fromIntegral (Word32 -> Word64) -> Word32 -> Word64
forall a b. (a -> b) -> a -> b
$ Float -> Word32
c_getnanf (CFloat -> Float
forall a b. (Real a, Fractional b) => a -> b
realToFrac CFloat
x)
    {-# INLINE nanPayload #-}
    minDenormal :: CFloat
minDenormal = CFloat
1e-45
    {-# INLINE minDenormal #-}
    minNormal :: CFloat
minNormal = CFloat
1.17549435e-38
    {-# INLINE minNormal #-}
    maxFinite :: CFloat
maxFinite = CFloat
3.40282347e+38
    {-# INLINE maxFinite #-}
    epsilon :: CFloat
epsilon = CFloat
1.19209290e-07
    {-# INLINE epsilon #-}
    copySign :: CFloat -> CFloat -> CFloat
copySign CFloat
x CFloat
y = Float -> CFloat
forall a b. (Real a, Fractional b) => a -> b
realToFrac (Float -> CFloat) -> Float -> CFloat
forall a b. (a -> b) -> a -> b
$ Float -> Float -> Float
c_copysignf (CFloat -> Float
forall a b. (Real a, Fractional b) => a -> b
realToFrac CFloat
x) (CFloat -> Float
forall a b. (Real a, Fractional b) => a -> b
realToFrac CFloat
y)
    {-# INLINE copySign #-}
    succIEEE :: CFloat -> CFloat
succIEEE CFloat
x = Float -> CFloat
forall a b. (Real a, Fractional b) => a -> b
realToFrac (Float -> CFloat) -> Float -> CFloat
forall a b. (a -> b) -> a -> b
$ Float -> Float
c_ieeesuccf (CFloat -> Float
forall a b. (Real a, Fractional b) => a -> b
realToFrac CFloat
x)
    {-# INLINE succIEEE #-}
    predIEEE :: CFloat -> CFloat
predIEEE CFloat
x = Float -> CFloat
forall a b. (Real a, Fractional b) => a -> b
realToFrac (Float -> CFloat) -> Float -> CFloat
forall a b. (a -> b) -> a -> b
$ Float -> Float
c_ieeepredf (CFloat -> Float
forall a b. (Real a, Fractional b) => a -> b
realToFrac CFloat
x)
    {-# INLINE predIEEE #-}
    bisectIEEE :: CFloat -> CFloat -> CFloat
bisectIEEE CFloat
x CFloat
y = Float -> CFloat
forall a b. (Real a, Fractional b) => a -> b
realToFrac (Float -> CFloat) -> Float -> CFloat
forall a b. (a -> b) -> a -> b
$ Float -> Float -> Float
c_ieeemeanf (CFloat -> Float
forall a b. (Real a, Fractional b) => a -> b
realToFrac CFloat
x) (CFloat -> Float
forall a b. (Real a, Fractional b) => a -> b
realToFrac CFloat
y)
    {-# INLINE bisectIEEE #-}
    sameSignificandBits :: CFloat -> CFloat -> Int
sameSignificandBits CFloat
x CFloat
y = Float -> Float -> Int
c_feqrelf (CFloat -> Float
forall a b. (Real a, Fractional b) => a -> b
realToFrac CFloat
x) (CFloat -> Float
forall a b. (Real a, Fractional b) => a -> b
realToFrac CFloat
y)
    {-# INLINE sameSignificandBits #-}


instance IEEE Double where
    identicalIEEE :: Double -> Double -> Bool
identicalIEEE Double
x Double
y = Double -> Double -> Int
c_identical Double
x Double
y Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
/= Int
0
    {-# INLINE identicalIEEE #-}
    infinity :: Double
infinity = Double
1Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/Double
0
    {-# INLINE infinity #-}
    nan :: Double
nan = (Double
0Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/Double
0)
    {-# INLINE nan #-}
    nanWithPayload :: Word64 -> Double
nanWithPayload Word64
n = Word64 -> Double
c_mknan Word64
n
    {-# INLINE nanWithPayload #-}
    maxNaNPayload :: Double -> Word64
maxNaNPayload Double
_ = Word64
0x0007FFFFFFFFFFFF
    {-# INLINE maxNaNPayload #-}
    nanPayload :: Double -> Word64
nanPayload Double
x = Double -> Word64
c_getnan Double
x
    {-# INLINE nanPayload #-}
    minDenormal :: Double
minDenormal = Double
5e-324
    {-# INLINE minDenormal #-}
    minNormal :: Double
minNormal = Double
2.2250738585072014e-308
    {-# INLINE minNormal #-}
    maxFinite :: Double
maxFinite = Double
1.7976931348623157e+308
    {-# INLINE maxFinite #-}
    epsilon :: Double
epsilon = Double
2.2204460492503131e-16
    {-# INLINE epsilon #-}
    copySign :: Double -> Double -> Double
copySign = Double -> Double -> Double
c_copysign
    {-# INLINE copySign #-}
    succIEEE :: Double -> Double
succIEEE = Double -> Double
c_ieeesucc
    {-# INLINE succIEEE #-}
    predIEEE :: Double -> Double
predIEEE = Double -> Double
c_ieeepred
    {-# INLINE predIEEE #-}
    bisectIEEE :: Double -> Double -> Double
bisectIEEE = Double -> Double -> Double
c_ieeemean
    {-# INLINE bisectIEEE #-}
    sameSignificandBits :: Double -> Double -> Int
sameSignificandBits = Double -> Double -> Int
c_feqrel
    {-# INLINE sameSignificandBits #-}


instance IEEE CDouble where
    identicalIEEE :: CDouble -> CDouble -> Bool
identicalIEEE CDouble
x CDouble
y = Double -> Double -> Int
c_identical (CDouble -> Double
forall a b. (Real a, Fractional b) => a -> b
realToFrac CDouble
x) (CDouble -> Double
forall a b. (Real a, Fractional b) => a -> b
realToFrac CDouble
y) Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
/= Int
0
    {-# INLINE identicalIEEE #-}
    infinity :: CDouble
infinity = CDouble
1CDouble -> CDouble -> CDouble
forall a. Fractional a => a -> a -> a
/CDouble
0
    {-# INLINE infinity #-}
    nan :: CDouble
nan = (CDouble
0CDouble -> CDouble -> CDouble
forall a. Fractional a => a -> a -> a
/CDouble
0)
    {-# INLINE nan #-}
    nanWithPayload :: Word64 -> CDouble
nanWithPayload Word64
n = Double -> CDouble
forall a b. (Real a, Fractional b) => a -> b
realToFrac (Double -> CDouble) -> Double -> CDouble
forall a b. (a -> b) -> a -> b
$ Word64 -> Double
c_mknan Word64
n
    {-# INLINE nanWithPayload #-}
    maxNaNPayload :: CDouble -> Word64
maxNaNPayload CDouble
_ = Word64
0x0007FFFFFFFFFFFF
    {-# INLINE maxNaNPayload #-}
    nanPayload :: CDouble -> Word64
nanPayload CDouble
x = Double -> Word64
c_getnan (CDouble -> Double
forall a b. (Real a, Fractional b) => a -> b
realToFrac CDouble
x)
    {-# INLINE nanPayload #-}
    minDenormal :: CDouble
minDenormal = CDouble
5e-324
    {-# INLINE minDenormal #-}
    minNormal :: CDouble
minNormal = CDouble
2.2250738585072014e-308
    {-# INLINE minNormal #-}
    maxFinite :: CDouble
maxFinite = CDouble
1.7976931348623157e+308
    {-# INLINE maxFinite #-}
    epsilon :: CDouble
epsilon = CDouble
2.2204460492503131e-16
    {-# INLINE epsilon #-}
    succIEEE :: CDouble -> CDouble
succIEEE CDouble
x = Double -> CDouble
forall a b. (Real a, Fractional b) => a -> b
realToFrac (Double -> CDouble) -> Double -> CDouble
forall a b. (a -> b) -> a -> b
$ Double -> Double
c_ieeesucc (CDouble -> Double
forall a b. (Real a, Fractional b) => a -> b
realToFrac CDouble
x)
    {-# INLINE succIEEE #-}
    copySign :: CDouble -> CDouble -> CDouble
copySign CDouble
x CDouble
y = Double -> CDouble
forall a b. (Real a, Fractional b) => a -> b
realToFrac (Double -> CDouble) -> Double -> CDouble
forall a b. (a -> b) -> a -> b
$ Double -> Double -> Double
c_copysign (CDouble -> Double
forall a b. (Real a, Fractional b) => a -> b
realToFrac CDouble
x) (CDouble -> Double
forall a b. (Real a, Fractional b) => a -> b
realToFrac CDouble
y)
    {-# INLINE copySign #-}
    predIEEE :: CDouble -> CDouble
predIEEE CDouble
x = Double -> CDouble
forall a b. (Real a, Fractional b) => a -> b
realToFrac (Double -> CDouble) -> Double -> CDouble
forall a b. (a -> b) -> a -> b
$ Double -> Double
c_ieeepred (CDouble -> Double
forall a b. (Real a, Fractional b) => a -> b
realToFrac CDouble
x)
    {-# INLINE predIEEE #-}
    bisectIEEE :: CDouble -> CDouble -> CDouble
bisectIEEE CDouble
x CDouble
y = Double -> CDouble
forall a b. (Real a, Fractional b) => a -> b
realToFrac (Double -> CDouble) -> Double -> CDouble
forall a b. (a -> b) -> a -> b
$ Double -> Double -> Double
c_ieeemean (CDouble -> Double
forall a b. (Real a, Fractional b) => a -> b
realToFrac CDouble
x) (CDouble -> Double
forall a b. (Real a, Fractional b) => a -> b
realToFrac CDouble
y)
    {-# INLINE bisectIEEE #-}
    sameSignificandBits :: CDouble -> CDouble -> Int
sameSignificandBits CDouble
x CDouble
y = Double -> Double -> Int
c_feqrel (CDouble -> Double
forall a b. (Real a, Fractional b) => a -> b
realToFrac CDouble
x) (CDouble -> Double
forall a b. (Real a, Fractional b) => a -> b
realToFrac CDouble
y)
    {-# INLINE sameSignificandBits #-}


foreign import ccall unsafe "identical"
    c_identical :: Double -> Double -> Int

foreign import ccall unsafe "identicalf"
    c_identicalf :: Float -> Float -> Int

foreign import ccall unsafe "feqrel"
    c_feqrel :: Double -> Double -> Int
foreign import ccall unsafe "feqrelf"
    c_feqrelf :: Float -> Float -> Int

foreign import ccall unsafe "ieeesucc"
    c_ieeesucc :: Double -> Double

foreign import ccall unsafe "ieeesuccf"
    c_ieeesuccf :: Float -> Float

foreign import ccall unsafe "ieeepred"
    c_ieeepred :: Double -> Double

foreign import ccall unsafe "ieeepredf"
    c_ieeepredf :: Float -> Float

foreign import ccall unsafe "ieeemean"
    c_ieeemean :: Double -> Double -> Double

foreign import ccall unsafe "ieeemeanf"
    c_ieeemeanf :: Float -> Float -> Float

foreign import ccall unsafe "copysign"
    c_copysign :: Double -> Double -> Double

foreign import ccall unsafe "copysignf"
    c_copysignf :: Float -> Float -> Float

foreign import ccall unsafe "mknan"
    c_mknan :: Word64 -> Double

foreign import ccall unsafe "getnan"
    c_getnan :: Double -> Word64

foreign import ccall unsafe "mknanf"
    c_mknanf :: Word32 -> Float

foreign import ccall unsafe "getnanf"
    c_getnanf :: Float -> Word32