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 -- To UNSUBSCRIBE, email to [EMAIL PROTECTED] with a subject of "unsubscribe". Trouble? Contact [EMAIL PROTECTED]