Copyright | (c) Levent Erkok |
---|---|

License | BSD3 |

Maintainer | erkokl@gmail.com |

Stability | experimental |

Safe Haskell | None |

Language | Haskell2010 |

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.

- assocPlus :: SFloat -> SFloat -> SFloat -> SBool
- assocPlusRegular :: IO ThmResult
- nonZeroAddition :: IO ThmResult
- multInverse :: IO ThmResult
- roundingAdd :: IO SatResult

# FP addition is not associative

assocPlus :: SFloat -> SFloat -> SFloat -> SBool Source #

Prove that floating point addition is not associative. For illustration purposes,
we will require one of the inputs to be a `NaN`

. We have:

`>>>`

Falsifiable. Counter-example: s0 = 0.0 :: Float s1 = 0.0 :: Float`prove $ assocPlus (0/0)`

Indeed:

`>>>`

`let i = 0/0 :: Float`

`>>>`

NaN`i + (0.0 + 0.0)`

`>>>`

NaN`((i + 0.0) + 0.0)`

But keep in mind that `NaN`

does not equal itself in the floating point world! We have:

`>>>`

False`let nan = 0/0 :: Float in nan == nan`

assocPlusRegular :: IO ThmResult Source #

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:

`>>>`

Falsifiable. Counter-example: x = 5.5511693e-15 :: Float y = -2.2204438e-16 :: Float z = -3.1763606e-21 :: Float`assocPlusRegular`

Indeed, we have:

`>>>`

5.3291218e-15`((5.5511693e-15) + ((-2.2204438e-16) + (-3.1763606e-21))) :: Float`

`>>>`

5.329122e-15`(((5.5511693e-15) + (-2.2204438e-16)) + (-3.1763606e-21)) :: Float`

Note the difference between two additions!

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

nonZeroAddition :: IO ThmResult Source #

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:

`>>>`

Falsifiable. Counter-example: a = 5.1705105e-26 :: Float b = -3.8518597e-34 :: Float`nonZeroAddition`

Indeed, we have:

`>>>`

True`(5.1705105e-26 + (-3.8518597e-34)) == (5.1705105e-26 :: Float)`

But:

`>>>`

False`-3.8518597e-34 == (0::Float)`

# FP multiplicative inverses may not exist

multInverse :: IO ThmResult Source #

This example illustrates that `a * (1/a)`

does not necessarily equal `1`

. Again,
we protect against division by `0`

and `NaN`

/`Infinity`

.

We have:

`>>>`

Falsifiable. Counter-example: a = 8.988465674311586e307 :: Double`multInverse`

Indeed, we have:

`>>>`

`let a = 8.988465674311586e307 :: Double`

`>>>`

1.0000000000000002`a * (1/a)`

# Effect of rounding modes

roundingAdd :: IO SatResult Source #

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:

`>>>`

Satisfiable. Model: rm = RoundTowardPositive :: RoundingMode x = 1.9531249 :: Float y = 6.2499996e-2 :: Float`roundingAdd`

(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.0156248`(1.9531249 + 6.2499996e-2) :: Float`

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:

`>>>`

Satisfiable. Model: s0 = 2.015625 :: Float`sat $ \z -> z .== fpAdd sRoundTowardPositive 1.9531249 (6.2499996e-2 :: SFloat)`

We can see why these two resuls are indeed different. But we can see that the `RoundTowardsPositive`

(which rounds towards positive-infinity) produces a larger result. Indeed, if we treat these numbers
as `Double`

values, we get:

`>>>`

2.015624896`(1.9531249 + 6.2499996e-2) :: Double`

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 https://ece.uwaterloo.ca/~dwharder/NumericalAnalysis/02Numerics/Double/paper.pdf as
an excellent guide.