Bug#372544: libc6-dev: fma() is incorrect (inaccurate), not conform to C99

[Vincent Lefevre <vincent@vinc17.org>]

Package: libc6-dev
Version: 2.3.6-15
Severity: normal

The fma() function is incorrect. For instance (this example is based
on one given in the ieee754r mailing-list):

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

int main (void)
{
  double eps, e2, f, x, z;

  eps = DBL_EPSILON;
  e2 = eps * eps;
  x = 1.0 + eps;
  z = -1.0 - 2.0 * eps;
  f = fma (x, x, z);
  if (f != e2)
    {
      fprintf (stderr, "fma() is WRONG on this platform!\n"
               "Got %a instead of %a.\n", f, e2);
      exit (EXIT_FAILURE);
    }
  return 0;
}

compiled with -std=c99 gives:

fma() is WRONG on this platform!
Got 0x0p+0 instead of 0x1p-104.

The cause is that fma() is implemented on the x86 with (x * y) + z.
The ISO C99 standard requires:

    The fma functions compute (x × y) + z, rounded as one ternary
    operation: they compute the value (as if) to infinite precision
    and round once to the result format, according to the rounding
    mode characterized by the value of FLT_ROUNDS.

I currently don't have the time to fix it, but I think this should
be done like this:

1. Compute the product as an exact sum of two FP numbers with Dekker's
   algorithm.
2. Determine the rounded results (see algorithms on FP expansions).

Well, there would be 2 problems: First, the fact that values may be in
extended precision (GCC bug). Then, the possible intermediate overflow
or underflow should be avoided, probably by doing comparisons first
and scale the values if need be (this problem is much less important
than the current situation).

Note: on PowerPC, which has a hardware fma, this bug does not appear.

-- System Information:
Debian Release: testing/unstable
  APT prefers unstable
  APT policy: (500, 'unstable'), (500, 'stable')
Architecture: i386 (i686)
Shell:  /bin/sh linked to /bin/bash
Kernel: Linux 2.6.14.4-20051215
Locale: LANG=POSIX, LC_CTYPE=en_US.ISO8859-1 (charmap=ISO-8859-1)

Versions of packages libc6-dev depends on:
ii  libc6                    2.3.6-15        GNU C Library: Shared libraries
ii  linux-kernel-headers     2.6.13+0rc3-2.1 Linux Kernel Headers for developme

Versions of packages libc6-dev recommends:
ii  gcc [c-compiler]              4:4.1.1-1  The GNU C compiler
ii  gcc-3.3 [c-compiler]          1:3.3.6-13 The GNU C compiler
ii  gcc-3.4 [c-compiler]          3.4.6-1    The GNU C compiler
ii  gcc-4.0 [c-compiler]          4.0.3-3    The GNU C compiler
ii  gcc-4.1 [c-compiler]          4.1.1-2    The GNU C compiler

-- no debconf information