-----------------------------------------------------------------------------
-- |
-- Copyright   :  (C) 2015 Anselm Jonas Scholl, (C) 2023 Julia Longtin
-- License     :  BSD3
-- Maintainer  :  Julia Longtin <Julia.longtin@gmail.com>
-- Stability   :  experimental
-- Portability :  GHC-specific
--
-- Provides increment-by-ulp, decrement-by-ulp, and get-ulp functions for
-- Doubles, and Floats.
----------------------------------------------------------------------------
{-# LANGUAGE ScopedTypeVariables #-}

module Data.Bits.Floating.Ulp (
     doubleNextUlp
    ,doublePrevUlp
    ,doubleUlp
    ,floatNextUlp
    ,floatPrevUlp
    ,floatUlp
    ) where

import Prelude (Double, Float, Integral, Num, RealFloat, Show, (<), (>), ($), (+), (-), (>=), (&&), (==), (||), abs, fromIntegral, isInfinite, isNaN, isNegativeZero, otherwise)
import Data.Int (Int64, Int32)
import Data.Bits (Bits, (.&.), shiftL, shiftR)
import Data.Bits.Floating.Prim (double2WordBitwise, float2WordBitwise, word2DoubleBitwise, word2FloatBitwise)

-- Tell HLint to ignore suggestions to eta-reduce expressions in this module.
{-# ANN module "HLint: ignore Eta reduce" #-}

---------------------
-- * Double constants
---------------------

doubleMinValue :: Double
doubleMinValue :: Double
doubleMinValue = Double
4.94065645841246544176568792868E-324

doubleMaxExponent :: Int64
doubleMaxExponent :: Int64
doubleMaxExponent = Int64
1023

doubleMinExponent :: Int64
doubleMinExponent :: Int64
doubleMinExponent = -Int64
1022

doubleSignificandWidth :: Int64
doubleSignificandWidth :: Int64
doubleSignificandWidth = Int64
53

doubleExpBias :: Int64
doubleExpBias :: Int64
doubleExpBias = Int64
1023

doubleExpBitMask :: Int64
doubleExpBitMask :: Int64
doubleExpBitMask = Int64
0x7FF0000000000000

--------------------
-- * Float constants
--------------------

floatMinValue :: Float
floatMinValue :: Float
floatMinValue = Float
1.40129846432481707092372958329E-45

floatMaxExponent :: Int32
floatMaxExponent :: Int32
floatMaxExponent = Int32
127

floatMinExponent :: Int32
floatMinExponent :: Int32
floatMinExponent = -Int32
126

floatSignificandWidth :: Int32
floatSignificandWidth :: Int32
floatSignificandWidth = Int32
24

floatExpBias :: Int32
floatExpBias :: Int32
floatExpBias = Int32
127

floatExpBitMask :: Int32
floatExpBitMask :: Int32
floatExpBitMask = Int32
0x7F800000

---------------------------
-- * Generic implementation
---------------------------

-- | Advance a RealFloat by one unit in the last place (ULP).
--
-- If the argument is NaN, NaN is returned.
-- If the argument is +INF, +INF is returned.
-- If the argument is (-)0.0, the minimum value greater than 0.0 is returned.
-- if the argument is -INF, -INF is returned.
{-# INLINE genericUp #-}
genericUp :: (RealFloat f, Num w) => (f -> w) -> (w -> f) -> f -> f
genericUp :: forall f w. (RealFloat f, Num w) => (f -> w) -> (w -> f) -> f -> f
genericUp f -> w
mkW w -> f
mkF f
f
    | forall a. RealFloat a => a -> Bool
isNaN f
f Bool -> Bool -> Bool
|| (forall a. RealFloat a => a -> Bool
isInfinite f
f Bool -> Bool -> Bool
&& f
f forall a. Ord a => a -> a -> Bool
> f
0.0) = f
f
    | forall a. RealFloat a => a -> Bool
isNegativeZero f
f = forall f w. (RealFloat f, Num w) => (f -> w) -> (w -> f) -> f -> f
genericUp f -> w
mkW w -> f
mkF f
0.0
    | Bool
otherwise = w -> f
mkF forall a b. (a -> b) -> a -> b
$ f -> w
mkW f
f forall a. Num a => a -> a -> a
+ (if f
f forall a. Ord a => a -> a -> Bool
>= f
0.0 then w
1 else -w
1)

{-# INLINE genericDown #-}
genericDown :: (RealFloat f, Num w) => (f -> w) -> (w -> f) -> f -> f -> f
genericDown :: forall f w.
(RealFloat f, Num w) =>
(f -> w) -> (w -> f) -> f -> f -> f
genericDown f -> w
mkW w -> f
mkF f
minValue f
f
    | forall a. RealFloat a => a -> Bool
isNaN f
f Bool -> Bool -> Bool
|| (forall a. RealFloat a => a -> Bool
isInfinite f
f Bool -> Bool -> Bool
&& f
f forall a. Ord a => a -> a -> Bool
< f
0.0) = f
f
    | f
f forall a. Eq a => a -> a -> Bool
== f
0.0  = -f
minValue
    | Bool
otherwise = w -> f
mkF forall a b. (a -> b) -> a -> b
$ f -> w
mkW f
f forall a. Num a => a -> a -> a
+ (if f
f forall a. Ord a => a -> a -> Bool
> f
0.0 then -w
1 else w
1)

{-# INLINE genericUlp #-}
genericUlp :: forall f i w . (RealFloat f, Bits w, Integral w, Bits i, Integral i, Show i) =>
    (f -> w) -> (w -> f) -> i -> i -> i -> i -> i -> f -> f -> f
genericUlp :: forall f i w.
(RealFloat f, Bits w, Integral w, Bits i, Integral i, Show i) =>
(f -> w) -> (w -> f) -> i -> i -> i -> i -> i -> f -> f -> f
genericUlp f -> w
mkW w -> f
mkF i
expBitMask i
significandWidth i
expBias i
maxExponent i
minExponent f
minValue f
f
    | i
expnt forall a. Eq a => a -> a -> Bool
== (i
maxExponent forall a. Num a => a -> a -> a
+ i
1) = forall a. Num a => a -> a
abs f
f -- NaN or INF
    | i
expnt forall a. Eq a => a -> a -> Bool
== (i
minExponent forall a. Num a => a -> a -> a
- i
1) = f
minValue
    | i
expnt2 forall a. Ord a => a -> a -> Bool
>= i
minExponent      = i -> f
powerOfTwo i
expnt2
--    | expnt > maxExponent        = error $ "Exponent too large: " ++ show expnt ++ " > " ++ show maxExponent
--    | expnt < minExponent        = error $ "Exponent too small: " ++ show expnt ++ " < " ++ show minExponent
    | Bool
otherwise                  = w -> f
mkF forall a b. (a -> b) -> a -> b
$ w
1 forall a. Bits a => a -> Int -> a
`shiftL` forall a b. (Integral a, Num b) => a -> b
fromIntegral (i
expnt2 forall a. Num a => a -> a -> a
- (i
minExponent forall a. Num a => a -> a -> a
- (i
significandWidth forall a. Num a => a -> a -> a
- i
1)))
    where
        expnt :: i
        expnt :: i
expnt = ((forall a b. (Integral a, Num b) => a -> b
fromIntegral (f -> w
mkW f
f) forall a. Bits a => a -> a -> a
.&. i
expBitMask) forall a. Bits a => a -> Int -> a
`shiftR` forall a b. (Integral a, Num b) => a -> b
fromIntegral (i
significandWidth forall a. Num a => a -> a -> a
- i
1)) forall a. Num a => a -> a -> a
- i
expBias
        expnt2 :: i
        expnt2 :: i
expnt2 = i
expnt forall a. Num a => a -> a -> a
- (i
significandWidth forall a. Num a => a -> a -> a
- i
1)
        powerOfTwo :: i -> f
        powerOfTwo :: i -> f
powerOfTwo i
n = w -> f
mkF forall a b. (a -> b) -> a -> b
$ forall a b. (Integral a, Num b) => a -> b
fromIntegral forall a b. (a -> b) -> a -> b
$ ((i
n forall a. Num a => a -> a -> a
+ i
expBias) forall a. Bits a => a -> Int -> a
`shiftL` forall a b. (Integral a, Num b) => a -> b
fromIntegral (i
significandWidth forall a. Num a => a -> a -> a
- i
1)) forall a. Bits a => a -> a -> a
.&. i
expBitMask

-----------------------------
-- * Specific implementations
-----------------------------

-- | Advance a 'Double' by one ULP.
{-# INLINABLE doubleNextUlp #-}
doubleNextUlp :: Double -> Double
doubleNextUlp :: Double -> Double
doubleNextUlp Double
d = forall f w. (RealFloat f, Num w) => (f -> w) -> (w -> f) -> f -> f
genericUp Double -> Word64
double2WordBitwise Word64 -> Double
word2DoubleBitwise Double
d

-- | Subtract one ULP from a 'Double'.
{-# INLINABLE doublePrevUlp #-}
doublePrevUlp :: Double -> Double
doublePrevUlp :: Double -> Double
doublePrevUlp Double
d = forall f w.
(RealFloat f, Num w) =>
(f -> w) -> (w -> f) -> f -> f -> f
genericDown Double -> Word64
double2WordBitwise Word64 -> Double
word2DoubleBitwise Double
doubleMinValue Double
d

-- | Return the distance to the next floating point number.
{-# INLINABLE doubleUlp #-}
doubleUlp :: Double -> Double
doubleUlp :: Double -> Double
doubleUlp Double
d = forall f i w.
(RealFloat f, Bits w, Integral w, Bits i, Integral i, Show i) =>
(f -> w) -> (w -> f) -> i -> i -> i -> i -> i -> f -> f -> f
genericUlp Double -> Word64
double2WordBitwise Word64 -> Double
word2DoubleBitwise Int64
doubleExpBitMask Int64
doubleSignificandWidth Int64
doubleExpBias Int64
doubleMaxExponent Int64
doubleMinExponent Double
doubleMinValue Double
d

-- | Advance a 'Float' by one ULP.
{-# INLINABLE floatNextUlp #-}
floatNextUlp :: Float -> Float
floatNextUlp :: Float -> Float
floatNextUlp Float
f = forall f w. (RealFloat f, Num w) => (f -> w) -> (w -> f) -> f -> f
genericUp Float -> Word32
float2WordBitwise Word32 -> Float
word2FloatBitwise Float
f

-- | Subtract one ULP from a 'Float'.
{-# INLINABLE floatPrevUlp #-}
floatPrevUlp :: Float -> Float
floatPrevUlp :: Float -> Float
floatPrevUlp Float
f = forall f w.
(RealFloat f, Num w) =>
(f -> w) -> (w -> f) -> f -> f -> f
genericDown Float -> Word32
float2WordBitwise Word32 -> Float
word2FloatBitwise Float
floatMinValue Float
f

-- | Return the distance to the next floating point number.
{-# INLINABLE floatUlp #-}
floatUlp :: Float -> Float
floatUlp :: Float -> Float
floatUlp Float
f = forall f i w.
(RealFloat f, Bits w, Integral w, Bits i, Integral i, Show i) =>
(f -> w) -> (w -> f) -> i -> i -> i -> i -> i -> f -> f -> f
genericUlp Float -> Word32
float2WordBitwise Word32 -> Float
word2FloatBitwise Int32
floatExpBitMask Int32
floatSignificandWidth Int32
floatExpBias Int32
floatMaxExponent Int32
floatMinExponent Float
floatMinValue Float
f