*To*: 372544@bugs.debian.org*Subject*: Bug#372544: libc6-dev: fma() is incorrect (inaccurate), not conform to C99*From*: Vincent Lefevre <vincent@vinc17.org>*Date*: Mon, 21 May 2007 13:28:32 +0200*Message-id*: <[🔎] 20070521112832.GA31247@vin.lip.ens-lyon.fr>*Reply-to*: Vincent Lefevre <vincent@vinc17.org>, 372544@bugs.debian.org*In-reply-to*: <20060610011151.GA32497@dixsept.loria.fr>*References*: <20060610011151.GA32497@dixsept.loria.fr>

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 <vincent@vinc17.org> - 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; }

- Prev by Date:
**Bug#425383: [INTL:lt] Lithuanian debconf templates translation** - Next by Date:
**No video on earth will enlighten me.** - Previous by thread:
**Bug#425383: marked as done ([INTL:lt] Lithuanian debconf templates translation)** - Next by thread:
**No video on earth will enlighten me.** - Index(es):