-----------------------------------------------------------------------------
-- |
-- 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'), double precision 'Double' ('SDouble'), and
-- the generic 'SFloatingPoint' @eb@ @sb@ type where the user can specify the
-- exponent and significand bit-widths. (Note that there is always an extra
-- sign-bit, and the value of @sb@ includes the hidden bit.)
--
-- 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 DataKinds           #-}
{-# LANGUAGE ScopedTypeVariables #-}

{-# OPTIONS_GHC -Wall -Werror #-}

module Documentation.SBV.Examples.Misc.Floating where

import Data.SBV

-- $setup
-- >>> -- For doctest purposes only:
-- >>> 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 :: SBV Float -> SBV Float -> SBV Float -> SBool
assocPlus SBV Float
x SBV Float
y SBV Float
z = SBV Float
x SBV Float -> SBV Float -> SBV Float
forall a. Num a => a -> a -> a
+ (SBV Float
y SBV Float -> SBV Float -> SBV Float
forall a. Num a => a -> a -> a
+ SBV Float
z) SBV Float -> SBV Float -> SBool
forall a. EqSymbolic a => a -> a -> SBool
.== (SBV Float
x SBV Float -> SBV Float -> SBV Float
forall a. Num a => a -> a -> a
+ SBV Float
y) SBV Float -> SBV Float -> SBV Float
forall a. Num a => a -> a -> a
+ SBV Float
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 =  2097149.9 :: Float
--   y =  13.999817 :: Float
--   z = -1.9998167 :: Float
--
-- Indeed, we have:
--
-- >>> let x =  2097149.9 :: Float
-- >>> let y =  13.999817 :: Float
-- >>> let z = -1.9998167 :: Float
-- >>> x + (y + z)
-- 2097162.0
-- >>> (x + y) + z
-- 2097161.8
--
-- 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 [SBV Float
x, SBV Float
y, SBV Float
z] <- [String] -> Symbolic [SBV Float]
sFloats [String
"x", String
"y", String
"z"]
                              let lhs :: SBV Float
lhs = SBV Float
xSBV Float -> SBV Float -> SBV Float
forall a. Num a => a -> a -> a
+(SBV Float
ySBV Float -> SBV Float -> SBV Float
forall a. Num a => a -> a -> a
+SBV Float
z)
                                  rhs :: SBV Float
rhs = (SBV Float
xSBV Float -> SBV Float -> SBV Float
forall a. Num a => a -> a -> a
+SBV Float
y)SBV Float -> SBV Float -> SBV Float
forall a. Num a => a -> a -> a
+SBV Float
z
                              -- make sure we do not overflow at the intermediate points
                              SBool -> SymbolicT IO ()
forall a. QuantifiedBool a => a -> SymbolicT IO ()
forall (m :: * -> *) a.
(SolverContext m, QuantifiedBool a) =>
a -> m ()
constrain (SBool -> SymbolicT IO ()) -> SBool -> SymbolicT IO ()
forall a b. (a -> b) -> a -> b
$ SBV Float -> SBool
forall a. IEEEFloating a => SBV a -> SBool
fpIsPoint SBV Float
lhs
                              SBool -> SymbolicT IO ()
forall a. QuantifiedBool a => a -> SymbolicT IO ()
forall (m :: * -> *) a.
(SolverContext m, QuantifiedBool a) =>
a -> m ()
constrain (SBool -> SymbolicT IO ()) -> SBool -> SymbolicT IO ()
forall a b. (a -> b) -> a -> b
$ SBV Float -> SBool
forall a. IEEEFloating a => SBV a -> SBool
fpIsPoint SBV Float
rhs
                              SBool -> SymbolicT IO SBool
forall a. a -> SymbolicT IO a
forall (m :: * -> *) a. Monad m => a -> m a
return (SBool -> SymbolicT IO SBool) -> SBool -> SymbolicT IO SBool
forall a b. (a -> b) -> a -> b
$ SBV Float
lhs SBV Float -> SBV Float -> SBool
forall a. EqSymbolic a => a -> a -> SBool
.== SBV Float
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 =   -7.441692e23 :: Float
--   b = -7.5039685e-27 :: Float
--
-- Indeed, we have:
--
-- >>> let a =   -7.441692e23 :: Float
-- >>> let b = -7.5039685e-27 :: 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 [SBV Float
a, SBV Float
b] <- [String] -> Symbolic [SBV Float]
sFloats [String
"a", String
"b"]
                             SBool -> SymbolicT IO ()
forall a. QuantifiedBool a => a -> SymbolicT IO ()
forall (m :: * -> *) a.
(SolverContext m, QuantifiedBool a) =>
a -> m ()
constrain (SBool -> SymbolicT IO ()) -> SBool -> SymbolicT IO ()
forall a b. (a -> b) -> a -> b
$ SBV Float -> SBool
forall a. IEEEFloating a => SBV a -> SBool
fpIsPoint SBV Float
a
                             SBool -> SymbolicT IO ()
forall a. QuantifiedBool a => a -> SymbolicT IO ()
forall (m :: * -> *) a.
(SolverContext m, QuantifiedBool a) =>
a -> m ()
constrain (SBool -> SymbolicT IO ()) -> SBool -> SymbolicT IO ()
forall a b. (a -> b) -> a -> b
$ SBV Float -> SBool
forall a. IEEEFloating a => SBV a -> SBool
fpIsPoint SBV Float
b
                             SBool -> SymbolicT IO ()
forall a. QuantifiedBool a => a -> SymbolicT IO ()
forall (m :: * -> *) a.
(SolverContext m, QuantifiedBool a) =>
a -> m ()
constrain (SBool -> SymbolicT IO ()) -> SBool -> SymbolicT IO ()
forall a b. (a -> b) -> a -> b
$ SBV Float
a SBV Float -> SBV Float -> SBV Float
forall a. Num a => a -> a -> a
+ SBV Float
b SBV Float -> SBV Float -> SBool
forall a. EqSymbolic a => a -> a -> SBool
.== SBV Float
a
                             SBool -> SymbolicT IO SBool
forall a. a -> SymbolicT IO a
forall (m :: * -> *) a. Monad m => a -> m a
return (SBool -> SymbolicT IO SBool) -> SBool -> SymbolicT IO SBool
forall a b. (a -> b) -> a -> b
$ SBV Float
b SBV Float -> SBV Float -> SBool
forall a. EqSymbolic a => a -> a -> SBool
.== SBV Float
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 = -5.69379e-39 :: Float
--
-- Indeed, we have:
--
-- >>> let a = -5.69379e-39 :: Float
-- >>> a * (1/a)
-- 0.99999994
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 SBV Float
a <- String -> Symbolic (SBV Float)
sFloat String
"a"
                         SBool -> SymbolicT IO ()
forall a. QuantifiedBool a => a -> SymbolicT IO ()
forall (m :: * -> *) a.
(SolverContext m, QuantifiedBool a) =>
a -> m ()
constrain (SBool -> SymbolicT IO ()) -> SBool -> SymbolicT IO ()
forall a b. (a -> b) -> a -> b
$ SBV Float -> SBool
forall a. IEEEFloating a => SBV a -> SBool
fpIsPoint SBV Float
a
                         SBool -> SymbolicT IO ()
forall a. QuantifiedBool a => a -> SymbolicT IO ()
forall (m :: * -> *) a.
(SolverContext m, QuantifiedBool a) =>
a -> m ()
constrain (SBool -> SymbolicT IO ()) -> SBool -> SymbolicT IO ()
forall a b. (a -> b) -> a -> b
$ SBV Float -> SBool
forall a. IEEEFloating a => SBV a -> SBool
fpIsPoint (SBV Float
1SBV Float -> SBV Float -> SBV Float
forall a. Fractional a => a -> a -> a
/SBV Float
a)
                         SBool -> SymbolicT IO SBool
forall a. a -> SymbolicT IO a
forall (m :: * -> *) a. Monad m => a -> m a
return (SBool -> SymbolicT IO SBool) -> SBool -> SymbolicT IO SBool
forall a b. (a -> b) -> a -> b
$ SBV Float
a SBV Float -> SBV Float -> SBV Float
forall a. Num a => a -> a -> a
* (SBV Float
1SBV Float -> SBV Float -> SBV Float
forall a. Fractional a => a -> a -> a
/SBV Float
a) SBV Float -> SBV Float -> SBool
forall a. EqSymbolic a => a -> a -> SBool
.== SBV Float
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 = RoundTowardZero :: RoundingMode
--   x  =  -1.0579198e-37 :: Float
--   y  =   4.7017266e-38 :: Float
--
-- (Note that depending on your version of Z3, you might get a different result.)
-- Unfortunately Haskell floats do not allow computation with arbitrary rounding modes, but SBV's
-- 'SFloatingPoint' type does. We have:
--
-- >>> fpAdd sRoundTowardZero (-1.0579198e-37) (4.7017266e-38) :: SFPSingle
-- -5.87747119e-38 :: SFloatingPoint 8 24
-- >>> fpAdd sRoundNearestTiesToEven (-1.0579198e-37) (4.7017266e-38) :: SFPSingle
-- -5.87747175e-38 :: SFloatingPoint 8 24
--
-- We can see why these two results are indeed different: The 'RoundTowardZero'
-- (which rounds towards zero) produces a larger result, closer to 0. Indeed, if we treat these numbers
-- as 'Double' values, we get:
--
-- >>> (-1.0579198e-37) + (4.7017266e-38) :: Double
-- -5.8774714e-38
--
-- we see that the "more precise" result is larger than what the 'Float' value is, justifying the
-- larger value with 'RoundTowardZero'. 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. Satisfiable 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 String
"rm"
                       SBool -> SymbolicT IO ()
forall a. QuantifiedBool a => a -> SymbolicT IO ()
forall (m :: * -> *) a.
(SolverContext m, QuantifiedBool a) =>
a -> 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
                       SBV Float
x <- String -> Symbolic (SBV Float)
sFloat String
"x"
                       SBV Float
y <- String -> Symbolic (SBV Float)
sFloat String
"y"
                       let lhs :: SBV Float
lhs = SRoundingMode -> SBV Float -> SBV Float -> SBV Float
forall a.
IEEEFloating a =>
SRoundingMode -> SBV a -> SBV a -> SBV a
fpAdd SRoundingMode
m SBV Float
x SBV Float
y
                       let rhs :: SBV Float
rhs = SBV Float
x SBV Float -> SBV Float -> SBV Float
forall a. Num a => a -> a -> a
+ SBV Float
y
                       SBool -> SymbolicT IO ()
forall a. QuantifiedBool a => a -> SymbolicT IO ()
forall (m :: * -> *) a.
(SolverContext m, QuantifiedBool a) =>
a -> m ()
constrain (SBool -> SymbolicT IO ()) -> SBool -> SymbolicT IO ()
forall a b. (a -> b) -> a -> b
$ SBV Float -> SBool
forall a. IEEEFloating a => SBV a -> SBool
fpIsPoint SBV Float
lhs
                       SBool -> SymbolicT IO ()
forall a. QuantifiedBool a => a -> SymbolicT IO ()
forall (m :: * -> *) a.
(SolverContext m, QuantifiedBool a) =>
a -> m ()
constrain (SBool -> SymbolicT IO ()) -> SBool -> SymbolicT IO ()
forall a b. (a -> b) -> a -> b
$ SBV Float -> SBool
forall a. IEEEFloating a => SBV a -> SBool
fpIsPoint SBV Float
rhs
                       SBool -> SymbolicT IO SBool
forall a. a -> SymbolicT IO a
forall (m :: * -> *) a. Monad m => a -> m a
return (SBool -> SymbolicT IO SBool) -> SBool -> SymbolicT IO SBool
forall a b. (a -> b) -> a -> b
$ SBV Float
lhs SBV Float -> SBV Float -> SBool
forall a. EqSymbolic a => a -> a -> SBool
./= SBV Float
rhs

-- | Arbitrary precision floating-point numbers. SBV can talk about floating point numbers with arbitrary
-- exponent and significand sizes as well. Here is a simple example demonstrating the minimum non-zero positive
-- and maximum floating point values with exponent width 5 and significand width 4, which is actually 3
-- bits for the significand explicitly stored, includes the hidden bit. We have:
--
-- >>> fp54Bounds
-- Objective "max": Optimal model:
--   x   = 61400 :: FloatingPoint 5 4
--   max =   503 :: WordN 9
--   min =   503 :: WordN 9
-- Objective "min": Optimal model:
--   x   = 0.00000763 :: FloatingPoint 5 4
--   max =        257 :: WordN 9
--   min =        257 :: WordN 9
--
-- The careful reader will notice that the numbers @61400@ and @0.00000763@ are quite suspicious, but the metric
-- space equivalents are correct. The reason for this is due to the sparcity of floats. The "computed" value of
-- the maximum in this bound is actually @61440@, however in @FloatingPoint 5 4@ representation all numbers
-- between @57344@ and @61440@ collapse to the same bit-pattern, and the pretty-printer picks a string
-- representation in decimal that falls within range that it considers is the "simplest." (Printing
-- floats precisely is a thorny subject!) Likewise, the minimum value we're looking for is actually
-- 2^-17, but any number between 2^-16 and 2^-17 will map to this number. It turns out that 0.00000763
-- in decimal is one such value. Moral of the story is that when reading floating-point numbers in
-- decimal notation one should be very careful about the printed representation and the numeric value; while
-- they will match in vsalue (if there are no bugs!), they can print quite differently! (Also keep in
-- mind the rounding modes that impact how the conversion is done.)
fp54Bounds :: IO OptimizeResult
fp54Bounds :: IO OptimizeResult
fp54Bounds = OptimizeStyle -> SymbolicT IO SBool -> IO OptimizeResult
forall a. Satisfiable a => OptimizeStyle -> a -> IO OptimizeResult
optimize OptimizeStyle
Independent (SymbolicT IO SBool -> IO OptimizeResult)
-> SymbolicT IO SBool -> IO OptimizeResult
forall a b. (a -> b) -> a -> b
$ do SFloatingPoint 5 4
x :: SFloatingPoint 5 4 <- String -> Symbolic (SFloatingPoint 5 4)
forall (eb :: Nat) (sb :: Nat).
ValidFloat eb sb =>
String -> Symbolic (SFloatingPoint eb sb)
sFloatingPoint String
"x"

                                       SBool -> SymbolicT IO ()
forall a. QuantifiedBool a => a -> SymbolicT IO ()
forall (m :: * -> *) a.
(SolverContext m, QuantifiedBool a) =>
a -> m ()
constrain (SBool -> SymbolicT IO ()) -> SBool -> SymbolicT IO ()
forall a b. (a -> b) -> a -> b
$ SFloatingPoint 5 4 -> SBool
forall a. IEEEFloating a => SBV a -> SBool
fpIsPoint SFloatingPoint 5 4
x
                                       SBool -> SymbolicT IO ()
forall a. QuantifiedBool a => a -> SymbolicT IO ()
forall (m :: * -> *) a.
(SolverContext m, QuantifiedBool a) =>
a -> m ()
constrain (SBool -> SymbolicT IO ()) -> SBool -> SymbolicT IO ()
forall a b. (a -> b) -> a -> b
$ SFloatingPoint 5 4
x SFloatingPoint 5 4 -> SFloatingPoint 5 4 -> SBool
forall a. OrdSymbolic a => a -> a -> SBool
.> SFloatingPoint 5 4
0

                                       String -> SFloatingPoint 5 4 -> SymbolicT IO ()
forall a. Metric a => String -> SBV a -> SymbolicT IO ()
maximize String
"max" SFloatingPoint 5 4
x
                                       String -> SFloatingPoint 5 4 -> SymbolicT IO ()
forall a. Metric a => String -> SBV a -> SymbolicT IO ()
minimize String
"min" SFloatingPoint 5 4
x

                                       SBool -> SymbolicT IO SBool
forall a. a -> SymbolicT IO a
forall (f :: * -> *) a. Applicative f => a -> f a
pure SBool
sTrue