#848 fma returns an incorrect result

[ARATA Mizuki]

The implementation of fma is not correct. I wrote a test program to check for some corner cases (attached; also available at https://gist.github.com/minoki/ff4412f2a385eab3c7e689acd589ea3e).

Here is an explanation for a part of the problem:

The current implementation of fma returns x * y + z if it's not finite.

https://sourceforge.net/p/mingw-w64/mingw-w64/ci/master/tree/mingw-w64-crt/math/fma.c#l77

  double ret = x * y + z;
  if(!isfinite(ret)) {
    return ret; /* If this naive check doesn't yield a finite value, the FMA isn't
                   likely to return one either. Forward the value as is. */
  }

But, there are cases where fma should return a finite value even if x * y + z is not finite; examples are

• x = 0x1p1000, y = 0x1p1000, z = -INFINITY: The product x * y overflows, but the exact value is still a finite value. Therefore, adding -INFINITY to it should yield -INFINITY, not a NaN.
• x = 0x0.ffffffep513, y = 0x1.0000002p511, z = -0x1p-1074: The exact value of x * y is 0x0.fffffffffffffcp1024 (54 bits of 1), which is 0.5 ulp above DBL_MAX. Therefore, adding a small negative number to it yields a number between DBL_MAX and DBL_MAX + 0.5ulp, and the rounded value would be DBL_MAX.

The rest of the problem is... well, the four statements after break_down should be reconsidered, I guess.

[Jameson]

You might want to pick up something like glibc (GNU) or openlibm (BSD) for the fma/fmaf/fmal algorithms from https://github.com/lattera/glibc/blob/master/sysdeps/ieee754/dbl-64/s_fma.c or
https://github.com/JuliaMath/openlibm/blob/master/src/s_fma.c
since you are missing a lot of edge cases for the fma algorithm that are critical to getting the rounding correct and avoiding double-rounding errors in the answer.

For example, this is another case we had come across in random testing, where mingw-w64 is giving the incorrect answer:
ASSERT(fmaf(-1.9369631e13f, 2.1513551e-7f, -1.7354427e-24f) == -4.1670958f6)