    -----------------------------------------------------------------------------
    -- |
    -- Module      :  Data.SBV.Examples.Misc.Floating
    -- Copyright   :  (c) Levent Erkok
    -- License     :  BSD3
    -- Maintainer  :  erkokl@gmail.com
    -- Stability   :  experimental
    --
    -- Several examples involving IEEE-754 floating point numbers, i.e., single
    -- precision 'Float' ('SFloat') and double precision 'Double' ('SDouble') types.
    --
    -- Note that arithmetic with floating point is full of surprises; due to precision
    -- issues associativity of arithmetic operations typically do not hold. Also,
    -- the presence of @NaN@ is always something to look out for.
    -----------------------------------------------------------------------------

    {-# LANGUAGE ScopedTypeVariables #-}

    module Data.SBV.Examples.Misc.Floating where

    import Data.SBV

    -----------------------------------------------------------------------------
    -- * FP addition is not associative
    -----------------------------------------------------------------------------

    -- | Prove that floating point addition is not associative. For illustration purposes,
    -- we will require one of the inputs to be a @NaN@. We have:
    --
    -- >>> prove $ assocPlus (0/0)
    -- Falsifiable. Counter-example:
    --   s0 = 0.0 :: Float
    --   s1 = 0.0 :: Float
    --
    -- Indeed:
    --
    -- >>> let i = 0/0 :: Float
    -- >>> i + (0.0 + 0.0)
    -- NaN
    -- >>> ((i + 0.0) + 0.0)
    -- NaN
    --
    -- But keep in mind that @NaN@ does not equal itself in the floating point world! We have:
    --
    -- >>> let nan = 0/0 :: Float in nan == nan
    -- False
    assocPlus :: SFloat -> SFloat -> SFloat -> SBool
    assocPlus x y z = x + (y + z) .== (x + y) + z

    -- | Prove that addition is not associative, even if we ignore @NaN@/@Infinity@ values.
    -- To do this, we use the predicate 'fpIsPoint', which is true of a floating point
    -- number ('SFloat' or 'SDouble') if it is neither @NaN@ nor @Infinity@. (That is, it's a
    -- representable point in the real-number line.)
    --
    -- We have:
    --
    -- >>> assocPlusRegular
    -- Falsifiable. Counter-example:
    --   x = -1.5991211e-2 :: Float
    --   y =     131071.99 :: Float
    --   z =    -131069.99 :: Float
    --
    -- Indeed, we have:
    --
    -- >>> ((-1.5991211e-2) + (131071.99 + (-131069.99))) :: Float
    -- 1.9840088
    -- >>> ((-1.5991211e-2) + 131071.99) + (-131069.99) :: Float
    -- 1.984375
    --
    -- Note the significant difference between two additions!
    assocPlusRegular :: IO ThmResult
    assocPlusRegular = prove $ do [x, y, z] <- sFloats ["x", "y", "z"]
                                  let lhs = x+(y+z)
                                      rhs = (x+y)+z
                                  -- make sure we do not overflow at the intermediate points
                                  constrain $ fpIsPoint lhs
                                  constrain $ fpIsPoint rhs
                                  return $ lhs .== rhs

    -----------------------------------------------------------------------------
    -- * FP addition by non-zero can result in no change
    -----------------------------------------------------------------------------

    -- | Demonstrate that @a+b = a@ does not necessarily mean @b@ is @0@ in the floating point world,
    -- even when we disallow the obvious solution when @a@ and @b@ are @Infinity.@
    -- We have:
    --
    -- >>> nonZeroAddition
    -- Falsifiable. Counter-example:
    --   a =  5.1705105e-26 :: Float
    --   b = -3.8518597e-34 :: Float
    --
    -- Indeed, we have:
    --
    -- >>> (5.1705105e-26 + (-3.8518597e-34)) == (5.1705105e-26 :: Float)
    -- True
    --
    -- But:
    --
    -- >>> -3.8518597e-34 == (0::Float)
    -- False
    --
    nonZeroAddition :: IO ThmResult
    nonZeroAddition = prove $ do [a, b] <- sFloats ["a", "b"]
                                 constrain $ fpIsPoint a
                                 constrain $ fpIsPoint b
                                 constrain $ a + b .== a
                                 return $ b .== 0

    -----------------------------------------------------------------------------
    -- * FP multiplicative inverses may not exist
    -----------------------------------------------------------------------------

    -- | This example illustrates that @a * (1/a)@ does not necessarily equal @1@. Again,
    -- we protect against division by @0@ and @NaN@/@Infinity@.
    --
    -- We have:
    --
    -- >>> multInverse
    -- Falsifiable. Counter-example:
    --   a = 1.1058928764217435e308 :: Double
    --
    -- Indeed, we have:
    --
    -- >>> let a = 1.1058928764217435e308 :: Double
    -- >>> a * (1/a)
    -- 0.9999999999999998
    multInverse :: IO ThmResult
    multInverse = prove $ do a <- sDouble "a"
                             constrain $ fpIsPoint a
                             constrain $ fpIsPoint (1/a)
                             return $ a * (1/a) .== 1

    -----------------------------------------------------------------------------
    -- * Effect of rounding modes
    -----------------------------------------------------------------------------

    -- | One interesting aspect of floating-point is that the chosen rounding-mode
    -- can effect the results of a computation if the exact result cannot be precisely
    -- represented. SBV exports the functions 'fpAdd', 'fpSub', 'fpMul', 'fpDiv', 'fpFMA'
    -- and 'fpSqrt' which allows users to specify the IEEE supported 'RoundingMode' for
    -- the operation. (Also see the class 'RoundingFloat'.) This example illustrates how SBV
    -- can be used to find rounding-modes where, for instance, addition can produce different
    -- results. We have:
    --
    -- >>> roundingAdd
    -- Satisfiable. Model:
    --   rm = RoundTowardZero :: RoundingMode
    --   x  =        7.984373 :: Float
    --   y  =  -2.1684042e-19 :: Float
    --
    -- (Note that depending on your version of Z3, you might get a different result.)
    -- Unfortunately we can't directly validate this result at the Haskell level, as Haskell only supports
    -- 'RoundNearestTiesToEven'. We have:
    --
    -- >>> (7.984373 + (-2.1684042e-19)) :: Float
    -- 7.984373
    --
    -- While we cannot directly see the result when the mode is 'RoundTowardZero' in Haskell, we can use
    -- SBV to provide us with that result thusly:
    --
    -- >>> sat $ \z -> z .== fpAdd sRoundTowardZero 7.984373 (-2.1684042e-19 :: SFloat)
    -- Satisfiable. Model:
    --   s0 = 7.9843726 :: Float
    --
    -- We can see why these two resuls are indeed different. To see why, one would have to convert the
    -- individual numbers to Float's, which would induce rounding-errors, add them up, and round-back;
    -- a tedious operation, but one that might prove illimunating for the interested reader. We'll merely
    -- note that floating point representation and semantics is indeed a thorny
    -- subject, and point to <https://ece.uwaterloo.ca/~dwharder/NumericalAnalysis/02Numerics/Double/paper.pdf> as
    -- an excellent guide.
    roundingAdd :: IO SatResult
    roundingAdd = sat $ do m :: SRoundingMode <- free "rm"
                           constrain $ m ./= literal RoundNearestTiesToEven
                           x <- sFloat "x"
                           y <- sFloat "y"
                           let lhs = fpAdd m x y
                           let rhs = x + y
                           constrain $ fpIsPoint lhs
                           constrain $ fpIsPoint rhs
                           return $ lhs ./= rhs