{-# LANGUAGE CPP #-}

-- | Logically consistent comparison of floating point numbers.
module Agda.Utils.Float
  ( asFinite
  , isPosInf
  , isNegInf
  , isPosZero
  , isNegZero
  , isSafeInteger
  , doubleEq
  , doubleLe
  , doubleLt
  , intToDouble
  , doublePlus
  , doubleMinus
  , doubleTimes
  , doubleNegate
  , doubleDiv
  , doublePow
  , doubleSqrt
  , doubleExp
  , doubleLog
  , doubleSin
  , doubleCos
  , doubleTan
  , doubleASin
  , doubleACos
  , doubleATan
  , doubleATan2
  , doubleSinh
  , doubleCosh
  , doubleTanh
  , doubleASinh
  , doubleACosh
  , doubleATanh
  , doubleRound
  , doubleFloor
  , doubleCeiling
  , doubleDenotEq
  , doubleDenotOrd
  , doubleToWord64
  , doubleToRatio
  , ratioToDouble
  , doubleDecode
  , doubleEncode
  , toStringWithoutDotZero
  ) where

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

import Agda.Utils.List  ( stripSuffix )

#if __GLASGOW_HASKELL__ >= 804
import GHC.Float (castDoubleToWord64, castWord64ToDouble)
#else
import System.IO.Unsafe (unsafePerformIO)
import qualified Foreign          as F
import qualified Foreign.Storable as F
#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

-- |Remove suffix @.0@ from printed floating point number.
toStringWithoutDotZero :: Double -> String
toStringWithoutDotZero :: Double -> Suffix Char
toStringWithoutDotZero Double
d = Suffix Char -> Maybe (Suffix Char) -> Suffix Char
forall a. a -> Maybe a -> a
fromMaybe Suffix Char
s (Maybe (Suffix Char) -> Suffix Char)
-> Maybe (Suffix Char) -> Suffix Char
forall a b. (a -> b) -> a -> b
$ Suffix Char -> Suffix Char -> Maybe (Suffix Char)
forall a. Eq a => Suffix a -> Suffix a -> Maybe (Suffix a)
stripSuffix Suffix Char
".0" Suffix Char
s
  where s :: Suffix Char
s = Double -> Suffix Char
forall a. Show a => a -> Suffix Char
show Double
d

-- |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