    -----------------------------------------------------------------------------
    -- |
    -- Module    : Documentation.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 #-}

    {-# OPTIONS_GHC -Wall -Werror #-}

    module Documentation.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 :: SFloat -> SFloat -> SFloat -> SBool
    assocPlus x :: SFloat
    x y :: SFloat
    y z :: SFloat
    z = SFloat
    x SFloat -> SFloat -> SFloat
    forall a. Num a => a -> a -> a
    + (SFloat
    y SFloat -> SFloat -> SFloat
    forall a. Num a => a -> a -> a
    + SFloat
    z) SFloat -> SFloat -> SBool
    forall a. EqSymbolic a => a -> a -> SBool
    .== (SFloat
    x SFloat -> SFloat -> SFloat
    forall a. Num a => a -> a -> a
    + SFloat
    y) SFloat -> SFloat -> SFloat
    forall a. Num a => a -> a -> a
    + SFloat
    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 =    128.00029 :: Float
    --   y =  -7.27236e-4 :: Float
    --   z = -6.875994e-3 :: Float
    --
    -- Indeed, we have:
    --
    -- >>> let x =    128.00029 :: Float
    -- >>> let y =  -7.27236e-4 :: Float
    -- >>> let z = -6.875994e-3 :: Float
    -- >>> x + (y + z)
    -- 127.99268
    -- >>> (x + y) + z
    -- 127.99269
    --
    -- Note the difference in the results!
    assocPlusRegular :: IO ThmResult
    assocPlusRegular :: IO ThmResult
    assocPlusRegular = SymbolicT IO SBool -> IO ThmResult
    forall a. Provable a => a -> IO ThmResult
    prove (SymbolicT IO SBool -> IO ThmResult)
    -> SymbolicT IO SBool -> IO ThmResult
    forall a b. (a -> b) -> a -> b
    $ do [x :: SFloat
    x, y :: SFloat
    y, z :: SFloat
    z] <- [String] -> Symbolic [SFloat]
    sFloats ["x", "y", "z"]
                                  let lhs :: SFloat
    lhs = SFloat
    xSFloat -> SFloat -> SFloat
    forall a. Num a => a -> a -> a
    +(SFloat
    ySFloat -> SFloat -> SFloat
    forall a. Num a => a -> a -> a
    +SFloat
    z)
                                      rhs :: SFloat
    rhs = (SFloat
    xSFloat -> SFloat -> SFloat
    forall a. Num a => a -> a -> a
    +SFloat
    y)SFloat -> SFloat -> SFloat
    forall a. Num a => a -> a -> a
    +SFloat
    z
                                  -- make sure we do not overflow at the intermediate points
                                  SBool -> SymbolicT IO ()
    forall (m :: * -> *). SolverContext m => SBool -> m ()
    constrain (SBool -> SymbolicT IO ()) -> SBool -> SymbolicT IO ()
    forall a b. (a -> b) -> a -> b
    $ SFloat -> SBool
    forall a. IEEEFloating a => SBV a -> SBool
    fpIsPoint SFloat
    lhs
                                  SBool -> SymbolicT IO ()
    forall (m :: * -> *). SolverContext m => SBool -> m ()
    constrain (SBool -> SymbolicT IO ()) -> SBool -> SymbolicT IO ()
    forall a b. (a -> b) -> a -> b
    $ SFloat -> SBool
    forall a. IEEEFloating a => SBV a -> SBool
    fpIsPoint SFloat
    rhs
                                  SBool -> SymbolicT IO SBool
    forall (m :: * -> *) a. Monad m => a -> m a
    return (SBool -> SymbolicT IO SBool) -> SBool -> SymbolicT IO SBool
    forall a b. (a -> b) -> a -> b
    $ SFloat
    lhs SFloat -> SFloat -> SBool
    forall a. EqSymbolic a => a -> a -> SBool
    .== SFloat
    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.060287e28 :: Float
    --   b = 3.6780381e19 :: Float
    --
    -- Indeed, we have:
    --
    -- >>> let a =  5.060287e28 :: Float
    -- >>> let b = 3.6780381e19 :: Float
    -- >>> a + b == a
    -- True
    -- >>> b == 0
    -- False
    nonZeroAddition :: IO ThmResult
    nonZeroAddition :: IO ThmResult
    nonZeroAddition = SymbolicT IO SBool -> IO ThmResult
    forall a. Provable a => a -> IO ThmResult
    prove (SymbolicT IO SBool -> IO ThmResult)
    -> SymbolicT IO SBool -> IO ThmResult
    forall a b. (a -> b) -> a -> b
    $ do [a :: SFloat
    a, b :: SFloat
    b] <- [String] -> Symbolic [SFloat]
    sFloats ["a", "b"]
                                 SBool -> SymbolicT IO ()
    forall (m :: * -> *). SolverContext m => SBool -> m ()
    constrain (SBool -> SymbolicT IO ()) -> SBool -> SymbolicT IO ()
    forall a b. (a -> b) -> a -> b
    $ SFloat -> SBool
    forall a. IEEEFloating a => SBV a -> SBool
    fpIsPoint SFloat
    a
                                 SBool -> SymbolicT IO ()
    forall (m :: * -> *). SolverContext m => SBool -> m ()
    constrain (SBool -> SymbolicT IO ()) -> SBool -> SymbolicT IO ()
    forall a b. (a -> b) -> a -> b
    $ SFloat -> SBool
    forall a. IEEEFloating a => SBV a -> SBool
    fpIsPoint SFloat
    b
                                 SBool -> SymbolicT IO ()
    forall (m :: * -> *). SolverContext m => SBool -> m ()
    constrain (SBool -> SymbolicT IO ()) -> SBool -> SymbolicT IO ()
    forall a b. (a -> b) -> a -> b
    $ SFloat
    a SFloat -> SFloat -> SFloat
    forall a. Num a => a -> a -> a
    + SFloat
    b SFloat -> SFloat -> SBool
    forall a. EqSymbolic a => a -> a -> SBool
    .== SFloat
    a
                                 SBool -> SymbolicT IO SBool
    forall (m :: * -> *) a. Monad m => a -> m a
    return (SBool -> SymbolicT IO SBool) -> SBool -> SymbolicT IO SBool
    forall a b. (a -> b) -> a -> b
    $ SFloat
    b SFloat -> SFloat -> SBool
    forall a. EqSymbolic a => a -> a -> SBool
    .== 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 = 2.4907063e38 :: Float
    --
    -- Indeed, we have:
    --
    -- >>> let a = 2.4907063e38 :: Float
    -- >>> a * (1/a)
    -- 1.0000001
    multInverse :: IO ThmResult
    multInverse :: IO ThmResult
    multInverse = SymbolicT IO SBool -> IO ThmResult
    forall a. Provable a => a -> IO ThmResult
    prove (SymbolicT IO SBool -> IO ThmResult)
    -> SymbolicT IO SBool -> IO ThmResult
    forall a b. (a -> b) -> a -> b
    $ do SFloat
    a <- String -> Symbolic SFloat
    sFloat "a"
                             SBool -> SymbolicT IO ()
    forall (m :: * -> *). SolverContext m => SBool -> m ()
    constrain (SBool -> SymbolicT IO ()) -> SBool -> SymbolicT IO ()
    forall a b. (a -> b) -> a -> b
    $ SFloat -> SBool
    forall a. IEEEFloating a => SBV a -> SBool
    fpIsPoint SFloat
    a
                             SBool -> SymbolicT IO ()
    forall (m :: * -> *). SolverContext m => SBool -> m ()
    constrain (SBool -> SymbolicT IO ()) -> SBool -> SymbolicT IO ()
    forall a b. (a -> b) -> a -> b
    $ SFloat -> SBool
    forall a. IEEEFloating a => SBV a -> SBool
    fpIsPoint (1SFloat -> SFloat -> SFloat
    forall a. Fractional a => a -> a -> a
    /SFloat
    a)
                             SBool -> SymbolicT IO SBool
    forall (m :: * -> *) a. Monad m => a -> m a
    return (SBool -> SymbolicT IO SBool) -> SBool -> SymbolicT IO SBool
    forall a b. (a -> b) -> a -> b
    $ SFloat
    a SFloat -> SFloat -> SFloat
    forall a. Num a => a -> a -> a
    * (1SFloat -> SFloat -> SFloat
    forall a. Fractional a => a -> a -> a
    /SFloat
    a) SFloat -> SFloat -> SBool
    forall a. EqSymbolic a => a -> a -> SBool
    .== 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. 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 = RoundTowardPositive :: RoundingMode
    --   x  =      -2.3509886e-38 :: Float
    --   y  =            -6.0e-45 :: 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:
    --
    -- >>> -2.3509886e-38 + (-6.0e-45) :: Float
    -- -2.3509893e-38
    --
    -- While we cannot directly see the result when the mode is 'RoundTowardPositive' in Haskell, we can use
    -- SBV to provide us with that result thusly:
    --
    -- >>> sat $ \z -> z .== fpAdd sRoundTowardPositive (-2.3509886e-38) (-6.0e-45 :: SFloat)
    -- Satisfiable. Model:
    --   s0 = -2.350989e-38 :: Float
    --
    -- We can see why these two resuls are indeed different: The 'RoundTowardPositive'
    -- (which rounds towards positive infinity from zero) produces a larger result. Indeed, if we treat these numbers
    -- as 'Double' values, we get:
    --
    -- >> -2.3509886e-38 + (-6.0e-45) :: Double
    -- -2.3509892e-38
    --
    -- we see that the "more precise" result is larger than what the 'Float' value is, justifying the
    -- larger value with 'RoundTowardPositive'. A more detailed study is beyond our current scope, so we'll
    -- merely note that floating point representation and semantics is indeed a thorny
    -- subject, and point to <http://ece.uwaterloo.ca/~dwharder/NumericalAnalysis/02Numerics/Double/paper.pdf> as
    -- an excellent guide.
    roundingAdd :: IO SatResult
    roundingAdd :: IO SatResult
    roundingAdd = SymbolicT IO SBool -> IO SatResult
    forall a. Provable a => a -> IO SatResult
    sat (SymbolicT IO SBool -> IO SatResult)
    -> SymbolicT IO SBool -> IO SatResult
    forall a b. (a -> b) -> a -> b
    $ do SRoundingMode
    m :: SRoundingMode <- String -> Symbolic SRoundingMode
    forall a. SymVal a => String -> Symbolic (SBV a)
    free "rm"
                           SBool -> SymbolicT IO ()
    forall (m :: * -> *). SolverContext m => SBool -> m ()
    constrain (SBool -> SymbolicT IO ()) -> SBool -> SymbolicT IO ()
    forall a b. (a -> b) -> a -> b
    $ SRoundingMode
    m SRoundingMode -> SRoundingMode -> SBool
    forall a. EqSymbolic a => a -> a -> SBool
    ./= RoundingMode -> SRoundingMode
    forall a. SymVal a => a -> SBV a
    literal RoundingMode
    RoundNearestTiesToEven
                           SFloat
    x <- String -> Symbolic SFloat
    sFloat "x"
                           SFloat
    y <- String -> Symbolic SFloat
    sFloat "y"
                           let lhs :: SFloat
    lhs = SRoundingMode -> SFloat -> SFloat -> SFloat
    forall a.
    IEEEFloating a =>
    SRoundingMode -> SBV a -> SBV a -> SBV a
    fpAdd SRoundingMode
    m SFloat
    x SFloat
    y
                           let rhs :: SFloat
    rhs = SFloat
    x SFloat -> SFloat -> SFloat
    forall a. Num a => a -> a -> a
    + SFloat
    y
                           SBool -> SymbolicT IO ()
    forall (m :: * -> *). SolverContext m => SBool -> m ()
    constrain (SBool -> SymbolicT IO ()) -> SBool -> SymbolicT IO ()
    forall a b. (a -> b) -> a -> b
    $ SFloat -> SBool
    forall a. IEEEFloating a => SBV a -> SBool
    fpIsPoint SFloat
    lhs
                           SBool -> SymbolicT IO ()
    forall (m :: * -> *). SolverContext m => SBool -> m ()
    constrain (SBool -> SymbolicT IO ()) -> SBool -> SymbolicT IO ()
    forall a b. (a -> b) -> a -> b
    $ SFloat -> SBool
    forall a. IEEEFloating a => SBV a -> SBool
    fpIsPoint SFloat
    rhs
                           SBool -> SymbolicT IO SBool
    forall (m :: * -> *) a. Monad m => a -> m a
    return (SBool -> SymbolicT IO SBool) -> SBool -> SymbolicT IO SBool
    forall a b. (a -> b) -> a -> b
    $ SFloat
    lhs SFloat -> SFloat -> SBool
    forall a. EqSymbolic a => a -> a -> SBool
    ./= SFloat
    rhs