Is the addition property of equality valid for floating point numbers?

user_2341285 0
Viewed 58

That is to say, if I have two floating point numbers (IEEE-754 binary floating-point) a and b, and I calculate c as follows:

c = a - b

Is this always true?

c + b == a

The answer to this simple question has been surprisingly difficult to find an answer to on the interwebs. The closest thing I found was this paper. Maybe I am simply not using the correct terminology.

algebra precision floating-point
2 Answers

You definitely need to add some conditions excluding overflow. This java program prints "false".

public strictfp class Test {
    public static void main(String[] args) {
        double a = Double.MAX_VALUE;
        double b = -Double.MAX_VALUE;
        double c = a-b;
        System.out.println(c+b == a);
    }
}

a-b in real number arithmetic is twice the maximum finite double. In floating point arithmetic it overflows to infinity.

user_17985930

My assumption is that the question pertains to computations with IEEE-754 binary floating-point arithmetic.

While tests cannot prove a hypothesized property, they can easily disprove one by finding counter examples. In this case we would not even need to go that far, as we can pick out a failing case by hand. Because IEEE-754 floating formats have finite precision, for non-zero a and b sufficiently different in magnitude with |a| < |b|, and in the absence of overflow and underflow, c == -b, then c+b == 0 != a.

The following simplistic ISO-C99 test program finds this and other cases:

#include <stdio.h>
#include <stdlib.h>
#include <stdint.h>
#include <string.h>
#include <float.h>
#include <math.h>

// Fixes via: Greg Rose, KISS: A Bit Too Simple. http://eprint.iacr.org/2011/007
static unsigned int z=362436069,w=521288629,jsr=362436069,jcong=123456789;
#define znew (z=36969*(z&0xffff)+(z>>16))
#define wnew (w=18000*(w&0xffff)+(w>>16))
#define MWC  ((znew<<16)+wnew)
#define SHR3 (jsr^=(jsr<<13),jsr^=(jsr>>17),jsr^=(jsr<<5)) /* 2^32-1 */
#define CONG (jcong=69069*jcong+13579)                     /* 2^32 */
#define KISS ((MWC^CONG)+SHR3)

float __uint32_as_float (uint32_t a)
{
    float r;
    memcpy (&r, &a, sizeof r);
    return r;
}

int main (void)
{
    const float ULMT = sqrtf (FLT_MAX) / 2; // avoid overflow
    const float LLMT = sqrtf (FLT_MIN) * 2; // avoid underflow
    const uint64_t N = 1ULL << 10;    
    uint64_t count = 0LL;
    uint32_t ai, bi;
    float af, bf, cf, sum;
   
    do {
        do {
            ai = KISS;
            af = __uint32_as_float (ai);
        } while (!isfinite(af) || (fabsf (af) > ULMT) || (fabsf (af) < LLMT));
        do {
            bi = KISS;
            bf = __uint32_as_float (bi);
        } while (!isfinite(bf) || (fabsf (bf) > ULMT) || (fabsf (bf) < LLMT));
        
        cf = af - bf;
        sum = cf + bf;
        if (sum != af) {
            printf ("a!=c+b:  a=% 15.6a  b=% 15.6a c=% 15.6a  b+c=% 15.6a\n", 
                    af, bf, cf, sum);
        }
        count++;
    } while (count < N);
    return EXIT_SUCCESS;
}
user_7807170
Related