Is it possible to proof a-b+b = a for all double floating-point numbers?

Question

I want to know whether the equation : a-b+b = a is always true for a, b belongs to double precision floating-point number and |a|>=|b|.

If the equation is true, how can I proof it?

If not, what precondition can make the equation right?

A counterexample shows the equation is not true, so I now want to know, the inequality |((a-b)+b)-a|<=ulp(a) is true or not, how can I proof the inequality?

I have tested the python code below for half hours, the inequality |a-b+b-a|<=ulp(a) seems hold.

def getulp(x):
x = float(x)
k = frexp(x)[1]-1
if x == 0.0:
    return pow(2, -1074)
if (k<1023)&(k>-1022):
    return pow(2,k-52)
else:
    return pow(2,-1074)



while 1> 0:
k = np.random.randint(-100, 100)
ub = np.float_power(2,k-1)
db = np.float_power(2,k)
j = np.random.randint(-100, 100)
ub2 = np.float_power(2, j - 1)
db2 = np.float_power(2, j)
x = np.random.uniform(ub,db,1)
y = np.random.uniform(ub2,db2,1)
if x>y:
    temp = x-y+y
    if np.fabs(temp-x)>getulp(x):
        print x
        print y
        print "The inequality false"
else:
    temp = y-x+x
    if np.fabs(temp-y)>getulp(y):
        print x
        print y
        print "The inequality false"

Thanks for the counterexample by @Kirill, I will explain why I need to proof the equation. I want to proof, image the a, b are two end points of a line, then if I know the value b and a-b, I want to proof that I can get the value of a without error larger than one ULP(a) error. So I can give the condition |a|>=|b| or |a|<=|b|, but I can not limit their actual value.

Answer 1

You can sometimes prove such results (or get counterexamples) using an SMT solver such as Z3 that supports floating point arithmetic. Here is a proof of a version of your theorem that says $|((x+y)-y)-x| \leq 2^{-23}|x|$ when $x>y>1$ and $x+y\neq\infty_{32}$ in 32-bit floating point arithmetic:

λ> import Data.SBV
λ> :set -XScopedTypeVariables
λ> prove $ \(x::SFloat) (y::SFloat) -> fpIsNormal (x+y) &&& x .> y &&& y .> 1 ==> fpAbs (((x+y)-y)-x) .< fpAbs x * 1.1920928955078125e-7
Q.E.D.

Here I used the sbv package in Haskell. If it holds for 32-bit floats, it should hold for 64-bit floats also.