{-# LANGUAGE CPP #-}

module MAlonzo.RTE.Float where

import Prelude (Bool, Double, Eq (..), Floating (..), Fractional (..), Functor (..), Int, Integer,
                Integral (..), Maybe (..), Num (..), Ord (..), Ordering (..), Real (..),
                RealFloat (..), RealFrac (..), even, fromIntegral, fst, otherwise, snd, uncurry,
                undefined, ($), (&&), (.), (^))

import Data.Bifunctor (bimap, second)
import Data.Function (on)
import Data.Maybe (fromMaybe)
import Data.Ratio (denominator, numerator, (%))
import Data.Word (Word64)

#if __GLASGOW_HASKELL__ >= 804
import GHC.Float (castDoubleToWord64, castWord64ToDouble)
#else
import Foreign qualified as F
import Foreign.Storable qualified as F
import System.IO.Unsafe (unsafePerformIO)
#endif

#if __GLASGOW_HASKELL__ < 804
castDoubleToWord64 :: Double -> Word64
castDoubleToWord64 float = unsafePerformIO $ F.alloca $ \buf -> do
  F.poke (F.castPtr buf) float
  F.peek buf

castWord64ToDouble :: Word64 -> Double
castWord64ToDouble word = unsafePerformIO $ F.alloca $ \buf -> do
  F.poke (F.castPtr buf) word
  F.peek buf
#endif

{-# INLINE doubleEq #-}
doubleEq :: Double -> Double -> Bool
doubleEq :: Double -> Double -> Bool
doubleEq = Double -> Double -> Bool
forall a. Eq a => a -> a -> Bool
(==)

{-# INLINE doubleLe #-}
doubleLe :: Double -> Double -> Bool
doubleLe :: Double -> Double -> Bool
doubleLe = Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
(<=)

{-# INLINE doubleLt #-}
doubleLt :: Double -> Double -> Bool
doubleLt :: Double -> Double -> Bool
doubleLt = Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
(<)

truncateDouble :: Double -> Double
truncateDouble :: Double -> Double
truncateDouble = Word64 -> Double
castWord64ToDouble (Word64 -> Double) -> (Double -> Word64) -> Double -> Double
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Double -> Word64
castDoubleToWord64

{-# INLINE intToDouble #-}
intToDouble :: Integral a => a -> Double
intToDouble :: forall a. Integral a => a -> Double
intToDouble = Double -> Double
truncateDouble (Double -> Double) -> (a -> Double) -> a -> Double
forall b c a. (b -> c) -> (a -> b) -> a -> c
. a -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral

{-# INLINE doublePlus #-}
doublePlus :: Double -> Double -> Double
doublePlus :: Double -> Double -> Double
doublePlus Double
x Double
y = Double -> Double
truncateDouble (Double
x Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
y)

{-# INLINE doubleMinus #-}
doubleMinus :: Double -> Double -> Double
doubleMinus :: Double -> Double -> Double
doubleMinus Double
x Double
y = Double -> Double
truncateDouble (Double
x Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
y)

{-# INLINE doubleTimes #-}
doubleTimes :: Double -> Double -> Double
doubleTimes :: Double -> Double -> Double
doubleTimes Double
x Double
y = Double -> Double
truncateDouble (Double
x Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
y)

{-# INLINE doubleNegate #-}
doubleNegate :: Double -> Double
doubleNegate :: Double -> Double
doubleNegate = Double -> Double
forall a. Num a => a -> a
negate -- NOTE: doesn't cause underflow/overflow

{-# INLINE doubleDiv #-}
doubleDiv :: Double -> Double -> Double
doubleDiv :: Double -> Double -> Double
doubleDiv = Double -> Double -> Double
forall a. Fractional a => a -> a -> a
(/) -- NOTE: doesn't cause underflow/overflow

{-# INLINE doublePow #-}
doublePow :: Double -> Double -> Double
doublePow :: Double -> Double -> Double
doublePow Double
x Double
y = Double -> Double
truncateDouble (Double
x Double -> Double -> Double
forall a. Floating a => a -> a -> a
** Double
y)

{-# INLINE doubleSqrt #-}
doubleSqrt :: Double -> Double
doubleSqrt :: Double -> Double
doubleSqrt = Double -> Double
forall a. Floating a => a -> a
sqrt -- NOTE: doesn't cause underflow/overflow

{-# INLINE doubleExp #-}
doubleExp :: Double -> Double
doubleExp :: Double -> Double
doubleExp Double
x = Double -> Double
truncateDouble (Double -> Double
forall a. Floating a => a -> a
exp Double
x)

{-# INLINE doubleLog #-}
doubleLog :: Double -> Double
doubleLog :: Double -> Double
doubleLog = Double -> Double
forall a. Floating a => a -> a
log -- NOTE: doesn't cause underflow/overflow

{-# INLINE doubleSin #-}
doubleSin :: Double -> Double
doubleSin :: Double -> Double
doubleSin = Double -> Double
forall a. Floating a => a -> a
sin -- NOTE: doesn't cause underflow/overflow

{-# INLINE doubleCos #-}
doubleCos :: Double -> Double
doubleCos :: Double -> Double
doubleCos = Double -> Double
forall a. Floating a => a -> a
cos -- NOTE: doesn't cause underflow/overflow

{-# INLINE doubleTan #-}
doubleTan :: Double -> Double
doubleTan :: Double -> Double
doubleTan = Double -> Double
forall a. Floating a => a -> a
tan -- NOTE: doesn't cause underflow/overflow

{-# INLINE doubleASin #-}
doubleASin :: Double -> Double
doubleASin :: Double -> Double
doubleASin = Double -> Double
forall a. Floating a => a -> a
asin -- NOTE: doesn't cause underflow/overflow

{-# INLINE doubleACos #-}
doubleACos :: Double -> Double
doubleACos :: Double -> Double
doubleACos = Double -> Double
forall a. Floating a => a -> a
acos -- NOTE: doesn't cause underflow/overflow

{-# INLINE doubleATan #-}
doubleATan :: Double -> Double
doubleATan :: Double -> Double
doubleATan = Double -> Double
forall a. Floating a => a -> a
atan -- NOTE: doesn't cause underflow/overflow

{-# INLINE doubleATan2 #-}
doubleATan2 :: Double -> Double -> Double
doubleATan2 :: Double -> Double -> Double
doubleATan2 = Double -> Double -> Double
forall a. RealFloat a => a -> a -> a
atan2 -- NOTE: doesn't cause underflow/overflow

{-# INLINE doubleSinh #-}
doubleSinh :: Double -> Double
doubleSinh :: Double -> Double
doubleSinh = Double -> Double
forall a. Floating a => a -> a
sinh -- NOTE: doesn't cause underflow/overflow

{-# INLINE doubleCosh #-}
doubleCosh :: Double -> Double
doubleCosh :: Double -> Double
doubleCosh = Double -> Double
forall a. Floating a => a -> a
cosh -- NOTE: doesn't cause underflow/overflow

{-# INLINE doubleTanh #-}
doubleTanh :: Double -> Double
doubleTanh :: Double -> Double
doubleTanh = Double -> Double
forall a. Floating a => a -> a
tanh -- NOTE: doesn't cause underflow/overflow

{-# INLINE doubleASinh #-}
doubleASinh :: Double -> Double
doubleASinh :: Double -> Double
doubleASinh = Double -> Double
forall a. Floating a => a -> a
asinh -- NOTE: doesn't cause underflow/overflow

{-# INLINE doubleACosh #-}
doubleACosh :: Double -> Double
doubleACosh :: Double -> Double
doubleACosh = Double -> Double
forall a. Floating a => a -> a
acosh -- NOTE: doesn't cause underflow/overflow

{-# INLINE doubleATanh #-}
doubleATanh :: Double -> Double
doubleATanh :: Double -> Double
doubleATanh = Double -> Double
forall a. Floating a => a -> a
atanh -- NOTE: doesn't cause underflow/overflow

{-# INLINE negativeZero #-}
negativeZero :: Double
negativeZero :: Double
negativeZero = -Double
0.0

positiveInfinity :: Double
positiveInfinity :: Double
positiveInfinity = Double
1.0 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
0.0

negativeInfinity :: Double
negativeInfinity :: Double
negativeInfinity = -Double
positiveInfinity

nan :: Double
nan :: Double
nan = Double
0.0 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
0.0

isPosInf :: Double -> Bool
isPosInf :: Double -> Bool
isPosInf Double
x = Double
x Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
> Double
0.0 Bool -> Bool -> Bool
&& Double -> Bool
forall a. RealFloat a => a -> Bool
isInfinite Double
x

isNegInf :: Double -> Bool
isNegInf :: Double -> Bool
isNegInf Double
x = Double
x Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
< Double
0.0 Bool -> Bool -> Bool
&& Double -> Bool
forall a. RealFloat a => a -> Bool
isInfinite Double
x

isPosZero :: Double -> Bool
isPosZero :: Double -> Bool
isPosZero Double
x = Double -> Double -> Bool
doubleDenotEq Double
x Double
0.0

isNegZero :: Double -> Bool
isNegZero :: Double -> Bool
isNegZero Double
x = Double -> Double -> Bool
doubleDenotEq Double
x (-Double
0.0)

doubleRound :: Double -> Maybe Integer
doubleRound :: Double -> Maybe Integer
doubleRound = (Double -> Integer) -> Maybe Double -> Maybe Integer
forall a b. (a -> b) -> Maybe a -> Maybe b
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
fmap Double -> Integer
forall b. Integral b => Double -> b
forall a b. (RealFrac a, Integral b) => a -> b
round (Maybe Double -> Maybe Integer)
-> (Double -> Maybe Double) -> Double -> Maybe Integer
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Double -> Maybe Double
asFinite

doubleFloor :: Double -> Maybe Integer
doubleFloor :: Double -> Maybe Integer
doubleFloor = (Double -> Integer) -> Maybe Double -> Maybe Integer
forall a b. (a -> b) -> Maybe a -> Maybe b
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
fmap Double -> Integer
forall b. Integral b => Double -> b
forall a b. (RealFrac a, Integral b) => a -> b
floor (Maybe Double -> Maybe Integer)
-> (Double -> Maybe Double) -> Double -> Maybe Integer
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Double -> Maybe Double
asFinite

doubleCeiling :: Double -> Maybe Integer
doubleCeiling :: Double -> Maybe Integer
doubleCeiling = (Double -> Integer) -> Maybe Double -> Maybe Integer
forall a b. (a -> b) -> Maybe a -> Maybe b
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
fmap Double -> Integer
forall b. Integral b => Double -> b
forall a b. (RealFrac a, Integral b) => a -> b
ceiling (Maybe Double -> Maybe Integer)
-> (Double -> Maybe Double) -> Double -> Maybe Integer
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Double -> Maybe Double
asFinite

normaliseNaN :: Double -> Double
normaliseNaN :: Double -> Double
normaliseNaN Double
x
  | Double -> Bool
forall a. RealFloat a => a -> Bool
isNaN Double
x   = Double
nan
  | Bool
otherwise = Double
x

doubleToWord64 :: Double -> Maybe Word64
doubleToWord64 :: Double -> Maybe Word64
doubleToWord64 Double
x
  | Double -> Bool
forall a. RealFloat a => a -> Bool
isNaN Double
x   = Maybe Word64
forall a. Maybe a
Nothing
  | Bool
otherwise = Word64 -> Maybe Word64
forall a. a -> Maybe a
Just (Double -> Word64
castDoubleToWord64 Double
x)

-- |Denotational equality for floating point numbers, checks bitwise equality.
--
--  NOTE: Denotational equality distinguishes NaNs, so its results may vary
--        depending on the architecture and compilation flags. Unfortunately,
--        this is a problem with floating-point numbers in general.
--
doubleDenotEq :: Double -> Double -> Bool
doubleDenotEq :: Double -> Double -> Bool
doubleDenotEq = Maybe Word64 -> Maybe Word64 -> Bool
forall a. Eq a => a -> a -> Bool
(==) (Maybe Word64 -> Maybe Word64 -> Bool)
-> (Double -> Maybe Word64) -> Double -> Double -> Bool
forall b c a. (b -> b -> c) -> (a -> b) -> a -> a -> c
`on` Double -> Maybe Word64
doubleToWord64

-- |I guess "denotational orderings" are now a thing? The point is that we need
--  an Ord instance which provides a total ordering, and is consistent with the
--  denotational equality.
--
--  NOTE: The ordering induced via `doubleToWord64` is total, and is consistent
--        with `doubleDenotEq`. However, it is *deeply* unintuitive. For one, it
--        considers all negative numbers to be larger than positive numbers.
--
doubleDenotOrd :: Double -> Double -> Ordering
doubleDenotOrd :: Double -> Double -> Ordering
doubleDenotOrd = Maybe Word64 -> Maybe Word64 -> Ordering
forall a. Ord a => a -> a -> Ordering
compare (Maybe Word64 -> Maybe Word64 -> Ordering)
-> (Double -> Maybe Word64) -> Double -> Double -> Ordering
forall b c a. (b -> b -> c) -> (a -> b) -> a -> a -> c
`on` Double -> Maybe Word64
doubleToWord64

-- |Return Just x if it's a finite number, otherwise return Nothing.
asFinite :: Double -> Maybe Double
asFinite :: Double -> Maybe Double
asFinite Double
x
  | Double -> Bool
forall a. RealFloat a => a -> Bool
isNaN      Double
x = Maybe Double
forall a. Maybe a
Nothing
  | Double -> Bool
forall a. RealFloat a => a -> Bool
isInfinite Double
x = Maybe Double
forall a. Maybe a
Nothing
  | Bool
otherwise    = Double -> Maybe Double
forall a. a -> Maybe a
Just Double
x

-- |Decode a Double to an integer ratio.
doubleToRatio :: Double -> (Integer, Integer)
doubleToRatio :: Double -> (Integer, Integer)
doubleToRatio Double
x
  | Double -> Bool
forall a. RealFloat a => a -> Bool
isNaN      Double
x = (Integer
0, Integer
0)
  | Double -> Bool
forall a. RealFloat a => a -> Bool
isInfinite Double
x = (Integer -> Integer
forall a. Num a => a -> a
signum (Double -> Integer
forall b. Integral b => Double -> b
forall a b. (RealFrac a, Integral b) => a -> b
floor Double
x), Integer
0)
  | Bool
otherwise    = let r :: Rational
r = Double -> Rational
forall a. Real a => a -> Rational
toRational Double
x in (Rational -> Integer
forall a. Ratio a -> a
numerator Rational
r, Rational -> Integer
forall a. Ratio a -> a
denominator Rational
r)

-- |Encode an integer ratio as a double.
ratioToDouble :: Integer -> Integer -> Double
ratioToDouble :: Integer -> Integer -> Double
ratioToDouble Integer
n Integer
d
  | Integer
d Integer -> Integer -> Bool
forall a. Eq a => a -> a -> Bool
== Integer
0 = case Integer -> Integer -> Ordering
forall a. Ord a => a -> a -> Ordering
compare Integer
n Integer
0 of
      Ordering
LT -> Double
negativeInfinity
      Ordering
EQ -> Double
nan
      Ordering
GT -> Double
positiveInfinity
  | Bool
otherwise = Rational -> Double
forall a. Fractional a => Rational -> a
fromRational (Integer
n Integer -> Integer -> Rational
forall a. Integral a => a -> a -> Ratio a
% Integer
d)

-- |Decode a Double to its mantissa and its exponent, normalised such that the
--  mantissa is the smallest possible number without loss of accuracy.
doubleDecode :: Double -> Maybe (Integer, Integer)
doubleDecode :: Double -> Maybe (Integer, Integer)
doubleDecode Double
x
  | Double -> Bool
forall a. RealFloat a => a -> Bool
isNaN      Double
x = Maybe (Integer, Integer)
forall a. Maybe a
Nothing
  | Double -> Bool
forall a. RealFloat a => a -> Bool
isInfinite Double
x = Maybe (Integer, Integer)
forall a. Maybe a
Nothing
  | Bool
otherwise    = (Integer, Integer) -> Maybe (Integer, Integer)
forall a. a -> Maybe a
Just ((Integer -> Integer -> (Integer, Integer))
-> (Integer, Integer) -> (Integer, Integer)
forall a b c. (a -> b -> c) -> (a, b) -> c
uncurry Integer -> Integer -> (Integer, Integer)
normalise ((Int -> Integer) -> (Integer, Int) -> (Integer, Integer)
forall b c a. (b -> c) -> (a, b) -> (a, c)
forall (p :: * -> * -> *) b c a.
Bifunctor p =>
(b -> c) -> p a b -> p a c
second Int -> Integer
forall a. Integral a => a -> Integer
toInteger (Double -> (Integer, Int)
forall a. RealFloat a => a -> (Integer, Int)
decodeFloat Double
x)))
  where
    normalise :: Integer -> Integer -> (Integer, Integer)
    normalise :: Integer -> Integer -> (Integer, Integer)
normalise Integer
mantissa Integer
exponent
      | Integer -> Bool
forall a. Integral a => a -> Bool
even Integer
mantissa = Integer -> Integer -> (Integer, Integer)
normalise (Integer
mantissa Integer -> Integer -> Integer
forall a. Integral a => a -> a -> a
`div` Integer
2) (Integer
exponent Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
+ Integer
1)
      | Bool
otherwise = (Integer
mantissa, Integer
exponent)

-- |Checks whether or not the Double is within a safe range of operation.
isSafeInteger :: Double -> Bool
isSafeInteger :: Double -> Bool
isSafeInteger Double
x = case Double -> (Integer, Double)
forall b. Integral b => Double -> (b, Double)
forall a b. (RealFrac a, Integral b) => a -> (b, a)
properFraction Double
x of
  (Integer
n, Double
f) -> Double
f Double -> Double -> Bool
forall a. Eq a => a -> a -> Bool
== Double
0.0 Bool -> Bool -> Bool
&& Integer
minMantissa Integer -> Integer -> Bool
forall a. Ord a => a -> a -> Bool
<= Integer
n Bool -> Bool -> Bool
&& Integer
n Integer -> Integer -> Bool
forall a. Ord a => a -> a -> Bool
<= Integer
maxMantissa

doubleRadix :: Integer
doubleRadix :: Integer
doubleRadix = Double -> Integer
forall a. RealFloat a => a -> Integer
floatRadix (Double
forall a. HasCallStack => a
undefined :: Double)

doubleDigits :: Int
doubleDigits :: Int
doubleDigits = Double -> Int
forall a. RealFloat a => a -> Int
floatDigits (Double
forall a. HasCallStack => a
undefined :: Double)

doubleRange :: (Int, Int)
doubleRange :: (Int, Int)
doubleRange = Double -> (Int, Int)
forall a. RealFloat a => a -> (Int, Int)
floatRange (Double
forall a. HasCallStack => a
undefined :: Double)

-- |The smallest representable mantissa. Simultaneously, the smallest integer which can be
--  represented as a Double without loss of precision.
minMantissa :: Integer
minMantissa :: Integer
minMantissa = - Integer
maxMantissa

-- |The largest representable mantissa. Simultaneously, the largest integer which can be
--  represented as a Double without loss of precision.
maxMantissa :: Integer
maxMantissa :: Integer
maxMantissa = (Integer
doubleRadix Integer -> Integer -> Integer
forall a b. (Num a, Integral b) => a -> b -> a
^ Int -> Integer
forall a. Integral a => a -> Integer
toInteger Int
doubleDigits) Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
- Integer
1

-- |The largest representable exponent.
minExponent :: Integer
minExponent :: Integer
minExponent = Int -> Integer
forall a. Integral a => a -> Integer
toInteger (Int -> Integer) -> Int -> Integer
forall a b. (a -> b) -> a -> b
$ ((Int, Int) -> Int
forall a b. (a, b) -> a
fst (Int, Int)
doubleRange Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
doubleDigits) Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1

-- |The smallest representable exponent.
maxExponent :: Integer
maxExponent :: Integer
maxExponent = Int -> Integer
forall a. Integral a => a -> Integer
toInteger (Int -> Integer) -> Int -> Integer
forall a b. (a -> b) -> a -> b
$ (Int, Int) -> Int
forall a b. (a, b) -> b
snd (Int, Int)
doubleRange Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
doubleDigits

-- |Encode a mantissa and an exponent as a Double.
doubleEncode :: Integer -> Integer -> Maybe Double
doubleEncode :: Integer -> Integer -> Maybe Double
doubleEncode Integer
mantissa Integer
exponent
  = if Integer
minMantissa Integer -> Integer -> Bool
forall a. Ord a => a -> a -> Bool
<= Integer
mantissa Bool -> Bool -> Bool
&& Integer
mantissa Integer -> Integer -> Bool
forall a. Ord a => a -> a -> Bool
<= Integer
maxMantissa Bool -> Bool -> Bool
&&
       Integer
minExponent Integer -> Integer -> Bool
forall a. Ord a => a -> a -> Bool
<= Integer
exponent Bool -> Bool -> Bool
&& Integer
exponent Integer -> Integer -> Bool
forall a. Ord a => a -> a -> Bool
<= Integer
maxExponent
    then Double -> Maybe Double
forall a. a -> Maybe a
Just (Integer -> Int -> Double
forall a. RealFloat a => Integer -> Int -> a
encodeFloat Integer
mantissa (Integer -> Int
forall a. Num a => Integer -> a
fromInteger Integer
exponent))
    else Maybe Double
forall a. Maybe a
Nothing