{-# 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 = forall a. Eq a => a -> a -> Bool
(==)

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

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

truncateDouble :: Double -> Double
truncateDouble :: Double -> Double
truncateDouble = Word64 -> Double
castWord64ToDouble 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 forall b c a. (b -> c) -> (a -> b) -> a -> c
. 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 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 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 forall a. Num a => a -> a -> a
* Double
y)

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

{-# INLINE doubleDiv #-}
doubleDiv :: Double -> Double -> Double
doubleDiv :: Double -> Double -> Double
doubleDiv = 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 forall a. Floating a => a -> a -> a
** Double
y)

{-# INLINE doubleSqrt #-}
doubleSqrt :: Double -> Double
doubleSqrt :: Double -> Double
doubleSqrt = 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 (forall a. Floating a => a -> a
exp Double
x)

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

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

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

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

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

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

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

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

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

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

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

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

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

{-# INLINE doubleATanh #-}
doubleATanh :: Double -> Double
doubleATanh :: Double -> Double
doubleATanh = 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 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 forall a. Fractional a => a -> a -> a
/ Double
0.0

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

isNegInf :: Double -> Bool
isNegInf :: Double -> Bool
isNegInf Double
x = Double
x forall a. Ord a => a -> a -> Bool
< Double
0.0 Bool -> Bool -> 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 = forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
fmap forall a b. (RealFrac a, Integral b) => a -> b
round forall b c a. (b -> c) -> (a -> b) -> a -> c
. Double -> Maybe Double
asFinite

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

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

normaliseNaN :: Double -> Double
normaliseNaN :: Double -> Double
normaliseNaN Double
x
  | 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
  | forall a. RealFloat a => a -> Bool
isNaN Double
x   = forall a. Maybe a
Nothing
  | Bool
otherwise = 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 = forall a. Eq a => a -> a -> 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 = forall a. Ord a => a -> a -> Ordering
compare 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
  | forall a. RealFloat a => a -> Bool
isNaN      Double
x = forall a. Maybe a
Nothing
  | forall a. RealFloat a => a -> Bool
isInfinite Double
x = forall a. Maybe a
Nothing
  | Bool
otherwise    = 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 = forall a. a -> Maybe a -> a
fromMaybe Suffix Char
s forall a b. (a -> b) -> a -> b
$ forall a. Eq a => Suffix a -> Suffix a -> Maybe (Suffix a)
stripSuffix Suffix Char
".0" Suffix Char
s
  where s :: Suffix Char
s = 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
  | forall a. RealFloat a => a -> Bool
isNaN      Double
x = (Integer
0, Integer
0)
  | forall a. RealFloat a => a -> Bool
isInfinite Double
x = (forall a. Num a => a -> a
signum (forall a b. (RealFrac a, Integral b) => a -> b
floor Double
x), Integer
0)
  | Bool
otherwise    = let r :: Rational
r = forall a. Real a => a -> Rational
toRational Double
x in (forall a. Ratio a -> a
numerator Rational
r, 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 forall a. Eq a => a -> a -> Bool
== Integer
0 = case 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 = forall a. Fractional a => Rational -> a
fromRational (Integer
n 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
  | forall a. RealFloat a => a -> Bool
isNaN      Double
x = forall a. Maybe a
Nothing
  | forall a. RealFloat a => a -> Bool
isInfinite Double
x = forall a. Maybe a
Nothing
  | Bool
otherwise    = forall a. a -> Maybe a
Just (forall a b c. (a -> b -> c) -> (a, b) -> c
uncurry Integer -> Integer -> (Integer, Integer)
normalise (forall (p :: * -> * -> *) b c a.
Bifunctor p =>
(b -> c) -> p a b -> p a c
second forall a. Integral a => a -> Integer
toInteger (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
      | forall a. Integral a => a -> Bool
even Integer
mantissa = Integer -> Integer -> (Integer, Integer)
normalise (Integer
mantissa forall a. Integral a => a -> a -> a
`div` Integer
2) (Integer
exponent 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 forall a b. (RealFrac a, Integral b) => a -> (b, a)
properFraction Double
x of
  (Integer
n, Double
f) -> Double
f forall a. Eq a => a -> a -> Bool
== Double
0.0 Bool -> Bool -> Bool
&& Integer
minMantissa forall a. Ord a => a -> a -> Bool
<= Integer
n Bool -> Bool -> Bool
&& Integer
n forall a. Ord a => a -> a -> Bool
<= Integer
maxMantissa

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

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

doubleRange :: (Int, Int)
doubleRange :: (Int, Int)
doubleRange = forall a. RealFloat a => a -> (Int, Int)
floatRange (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 forall a b. (Num a, Integral b) => a -> b -> a
^ forall a. Integral a => a -> Integer
toInteger Int
doubleDigits) forall a. Num a => a -> a -> a
- Integer
1

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

-- |The smallest representable exponent.
maxExponent :: Integer
maxExponent :: Integer
maxExponent = forall a. Integral a => a -> Integer
toInteger forall a b. (a -> b) -> a -> b
$ forall a b. (a, b) -> b
snd (Int, Int)
doubleRange 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 forall a. Ord a => a -> a -> Bool
<= Integer
mantissa Bool -> Bool -> Bool
&& Integer
mantissa forall a. Ord a => a -> a -> Bool
<= Integer
maxMantissa Bool -> Bool -> Bool
&&
       Integer
minExponent forall a. Ord a => a -> a -> Bool
<= Integer
exponent Bool -> Bool -> Bool
&& Integer
exponent forall a. Ord a => a -> a -> Bool
<= Integer
maxExponent
    then forall a. a -> Maybe a
Just (forall a. RealFloat a => Integer -> Int -> a
encodeFloat Integer
mantissa (forall a. Num a => Integer -> a
fromInteger Integer
exponent))
    else forall a. Maybe a
Nothing