Associativity in floating point arithmetic failing by two values

Question

Cross-posting from math.stackexchange, since there might be people here familiar with this topic.

---

Assume working in floating point arithmetic with finite precision, bounded exponent and rounding to nearest.

Let $x,y$ be positive. It is not hard to find examples such that

$$(x+y)-x > y$$

In order to construct a test for a software I wanted to find an example or prove that such example doesn't exist of following:

Let $s(x)$ be the successor of $x$.

Is it possible to have both $$\begin{align}(x+y)-x&>y\\(x+y)-s(x)&>y\end{align}$$

It is easy to write a program that searches for such an example, but also it is unfeasible to test all possibilities and show that the example doesn't exist in this manner. So far my code hasn't got any example.

Example: In case seeing an example of $(x+y)-x>y$ helps somehow, take $$ \begin{align} x&=1.1234567891234568\\ y&=1e-5\text{ ( denoting }10^{-5}) \end{align} $$ Then $(x+y)-x=1.0000000000065512e-05 > y$. Examples of the first inequality there are many.

Answer 1

These kinds of floating point counterexamples are notoriously tricky to discover by hand. Usually you do expect random brute-force search to work, especially on 32-bit floats, but sometimes it fails.

These days, you have SMT solvers, like Z3 available to you, so you can try formulating your formula in Z3 to see if it can find a counterexample. This is the "high-tech" method of solving your problem. This approach is nice because it doesn't require you to know very much about floating point arithmetic to be able to use it.

I tried this, using Haskell, with the SBV package using Z3 as the backend. It looks only at the range $x\in(1,2)$, where $s(x) = x + \epsilon$

import Data.SBV

ε = 1.1920928955078125e-07

f :: SFloat -> SFloat -> SBool
f x y = 1 .< x &&& x .< 2 &&& fpIsNormal y &&& (x + y) - (x + ε) .> y

and it immediately found this counterexample:

λ> sat f
Satisfiable. Model:
  s0 =       1.5 :: Float
  s1 = 3.9999998 :: Float

This satisfies $(x+y)-s(x) > y$:

julia> x, y = 1.5f0, 3.9999998f0
julia> ((x + y) - nextfloat(x)) - y
2.3841858f-7
julia> (x + y) - nextfloat(x) > y
true

Note Looking at your example, it's unclear to me if you have restrictions on the magnitudes on $x$ and $y$. For example, if you restrict it to $x\in(1,2)$, $y\in(0,1)$, then there is no way to satisfy $(x+y)-s(x)>y$. But restricting only $x\in(1,2)$, but not $y$, it seems you get very large $n$ with $(x+y)-(x+n\epsilon) > y$.

Looking at some small examples, it seems the worst case can be reached by taking (in binary, single precision) $$ x = 1.5 = (1.1)_2, \qquad y = 2^k\times (0.111111111111111111111111)_2 = 2^k\times (1-2^{-24}), $$ which gives $n=2^{k-1}$, so $(x+y)-(x+2^{k-1}\epsilon)=y$. In double precision it works to take $x=1.5$ and $y=(1-2^{-53})\times 2^k$ (which is just all ones in binary), and the relationship only starts to break down from $k\geq51$.