Bug#372544: libc6-dev: fma() is incorrect (inaccurate), not conform to C99
Note: this is upstream bug 3268: http://sourceware.org/bugzilla/show_bug.cgi?id=3268 -- Vincent Lefèvre [EMAIL PROTECTED] - Web: http://www.vinc17.org/ 100% accessible validated (X)HTML - Blog: http://www.vinc17.org/blog/ Work: CR INRIA - computer arithmetic / Arenaire project (LIP, ENS-Lyon)
Bug#372544: libc6-dev: fma() is incorrect (inaccurate), not conform to C99
On 2006-06-10 03:11:51 +0200, Vincent Lefevre wrote:
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).
The overflow problem occurs on x86_64 (since SSE2 is used); see
attached C file.
I haven't tested it, but I assume that it also occurs in fmal
on all x86 machines.
--
Vincent Lefèvre [EMAIL PROTECTED] - Web: http://www.vinc17.org/
100% accessible validated (X)HTML - Blog: http://www.vinc17.org/blog/
Work: CR INRIA - computer arithmetic / Arenaire project (LIP, ENS-Lyon)
/* $Id: fma-tests.c 17487 2007-05-21 11:03:03Z lefevre $
*
* Miscellaneous fma() test.
*/
#include stdio.h
#include stdlib.h
#include float.h
#include math.h
#define WRONG(S,X,Y,Z,F,C)\
fprintf (stderr, ERROR in test ' S '!\nfma (%a, %a, %a)\n \
returned %a instead of %a.\n, X, Y, Z, F, C)
/* Modified Nelson H. F. Beebe's fma() test from stds-754 list, 2006-06-09 */
static int beebe (void)
{
volatile 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)
{
WRONG (beebe, x, x, z, f, e2);
return 1;
}
return 0;
}
static int overflowed_mul (void)
{
volatile double x, f;
x = DBL_MAX;
f = fma (x, 2.0, -x);
if (f != x)
{
WRONG (overflowed_mul, x, 2.0, -x, f, x);
return 1;
}
return 0;
}
int main (void)
{
int err = 0;
err += beebe ();
err += overflowed_mul ();
return err;
}
Bug#372544: libc6-dev: fma() is incorrect (inaccurate), not conform to C99
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 libc62.3.6-15GNU 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-1The GNU C compiler
ii gcc-4.0 [c-compiler] 4.0.3-3The GNU C compiler
ii gcc-4.1 [c-compiler] 4.1.1-2The GNU C compiler
-- no debconf information
--
To UNSUBSCRIBE, email to [EMAIL PROTECTED]
with a subject of unsubscribe. Trouble? Contact [EMAIL PROTECTED]
