tree 01b518e9b2eba1adc2c9982cac6b378201e0d685
parent 49e9a673b36a1241747bf3ea95d040f127753f81
author Liu Hao <lh_mouse@126.com> 1576155844 +0800
committer Liu Hao <lh_mouse@126.com> 1576745552 +0800

crt/fma{,f}.c: Implement FMA for `double` and `float` properly

FMA of `double` shall not be implemented using `long double`, as
the result has 106 bits and cannot be stored accurately in a
`long double`, whose mantissa only has 64 bits.

Let's calculate `FMA(0x1.0000000000001p0, 1, 0x0.00000000000007FFp0)`:
First, we multiply `0x1.0000000000001p0` and `1`, which yields
`0x1.0000000000001p0`. Then we add `0x0.00000000000007FFp0` to it:

      0x1.0000 0000 0000 1     p0
  +)  0x0.0000 0000 0000 07FF  p0
  ---------------------------------
      0x1.0000 0000 0000 17FF  p0   (long double)

This result has 65 significant bits. When it is stored into a
`long double`, the last bit is rounded to even, as follows:

      0x1.0000 0000 0000 17FF  p0
                            1       <= rounded UP as this bit is set
  ---------------------------------
      0x1.0000 0000 0000 1800  p0   (long double)

If we attempt to round this value into `double` again, we get:

      0x1.0000 0000 0000 1800  p0
                          8         <= rounded UP as this bit is set
  ---------------------------------
      0x1.0000 0000 0000 2000  p0   (double)

This is wrong. Because FMA shall round the result exactly once. If we
round the first result to double we get a different result:

      0x0.0000 0000 0000 17FF  p0
                          8         <= rounded DOWN as this bit is clear
  ---------------------------------
      0x1.0000 0000 0000 1000  p0   (double)

`float` suffers from a similar problem.

This patch fixes such issues.

Testcase for `double`:

  #include <assert.h>
  #include <math.h>

  int main(void)
    {
      volatile double a = 0x1.0000000000001p0;
      volatile double b = 1;
      volatile double c = 0x0.00000000000007FFp0;

      assert(a * b + c    == 0x1.0000000000001p0);
      assert(fma(a, b, c) == 0x1.0000000000001p0);
      assert((double)((long double)a * b + c)
                          == 0x1.0000000000002p0);
    }

Testcase for `float`:

  #include <assert.h>
  #include <math.h>

  int main(void)
    {
      volatile float a = 0x0.800001p0f;
      volatile float b = 0x0.5FFFFFp0f;;
      volatile float c = 0x0.000000000000FFFFFFp0f;

      assert(a * b + c     == 0x0.2FFFFFCp0);
      assert(fmaf(a, b, c) == 0x0.2FFFFFCp0);
      assert((float)((double)a * b + c)
                           == 0x0.3000000p0);
    }

Reference: https://www.exploringbinary.com/double-rounding-errors-in-floating
Signed-off-by: Liu Hao <lh_mouse@126.com>