Hide Forgot
From Bugzilla Helper: User-Agent: Mozilla/4.77 [en] (X11; U; Linux 2.4.2-2 i586) Description of problem: fma() is currently implemented as two binary operations, which may result in a loss of precision as x*y is rounded before adding z. How reproducible: Always Steps to Reproduce: Compile and run the following program: #include <stdio.h> #define __USE_ISOC99 1 #include <math.h> main () { const double x = 4660.3377777777768642408773303032; const double y = 8173.7955555555554383317939937115; const double z = -38092648.215387649834156036376953125; double a; a = fma (x, y, z); printf ("%f * %f + %f = %.9g\n", x, y, z, a); printf ("True result rounds to %.9g\n", -3.526534563229407686121644747209e-9); } Actual Results: 4660.337778 * 8173.795556 + -38092648.215388 = -3.52520146e-09 True result rounds to -3.52653456e-09 Expected Results: 4660.337778 * 8173.795556 + -38092648.215388 = -3.52653456e-09 True result rounds to -3.52653456e-09 Additional info: The fma() function is new in ISO C99. According to WG14/N869 7.12.13.1: 2 "The fma functions compute the sum z plus the product x times y, rounded as one ternary operation: they computes the sum z plus the product x times y (as if) to infinite precision and round once to the result format, according to the rounding mode characterized by the value of FLT_ROUNDS." A correct implementation requires the intermediate product x*y to have at least double the number of significant bits of the operands. This is easily achieved in the float version of fmaf by computing the intermediate result using double precision, but with double and long double, you need greater precision than most hardware can handle.
On ia64 it works correctly, the rest of architectures need better fma implementations. I think the current one can well stay in sysdeps/generic, plus there should be fmaf which does the computation on double (thus is FP_FAST_FMA) - 53 is more than twice as large as 24, so it will be correct. For fma and fmal somebody (perhaps me, but not this week) will need to write a 0.5 ulp precise implementation.
I've finished coding and testing an i386-specific implementation that passes all of the tests I gave it. It's written in assembly, and uses the Micros**t AsSeMbler syntax (which is what my employer's custom O/S is written in), but I think it shouldn't be hard to port to other architectures or generic C code. Feel free to study the code and borrow any ideas you like.
Created attachment 20593 [details] Fused multiply and add - long double version for i386, written for MASM
Created attachment 20594 [details] Test cases for fmal() (i386-specific). If anyone has ideas for additional tests, let me know. The above assembly file passes all these tests.
We have already a framework for this in glibc (soft-fp), I was just trying to think of an implementation which would work without the need to decompose the numbers, do the math in and compose them again. But you're right, it is probably best to code the generic fma and fmal using soft-fp (after checking for exceptional values) and arches can then optimize if they want.
Red Hat Linux is no longer supported by Red Hat, Inc. If you are still running Red Hat Linux, you are strongly advised to upgrade to a current Fedora Core release or Red Hat Enterprise Linux or comparable. Some information on which option may be right for you is available at http://www.redhat.com/rhel/migrate/redhatlinux/. Red Hat apologizes that these issues have not been resolved yet. We do want to make sure that no important bugs slip through the cracks. Please check if this issue is still present in a current Fedora Core release. If so, please change the product and version to match, and check the box indicating that the requested information has been provided. Note that any bug still open against Red Hat Linux on will be closed as 'CANTFIX' on September 30, 2006. Thanks again for your help.
I've verified that the bug is still present in Fedora Core 6 RC2 (glibc-2.4.90-17).
Fedora Core 5 and Fedora Core 6 are, as we're sure you've noticed, no longer test releases. We're cleaning up the bug database and making sure important bug reports filed against these test releases don't get lost. It would be helpful if you could test this issue with a released version of Fedora or with the latest development / test release. Thanks for your help and for your patience. [This is a bulk message for all open FC5/FC6 test release bugs. I'm adding myself to the CC list for each bug, so I'll see any comments you make after this and do my best to make sure every issue gets proper attention.]
Using the same program, I get different results this time. It's worse: [trevin@hydra src]$ gcc -lm -o bug-43358{,.c} [trevin@hydra src]$ ./bug-43358 4660.337778 * 8173.795556 + -38092648.215388 = 0 True result rounds to -3.52653456e-09 It may make a difference that my current computer is an AMD64 system, rather than the Athlon I had before (it died) or the K6-2 I had when I first reported the bug. So I compiled and ran the program again on a virtual machine running the 32-bit version of FC6: [trevin@fc6-i386 src]$ gcc -lm -o bug-43358{,.c} [trevin@fc6-i386 src]$ ./bug-43358 4660.337778 * 8173.795556 + -38092648.215388 = -3.52520146e-09 True result rounds to -3.52653456e-09 Well, what do you know? No change from the original report. So it looks like this bug exists in two forms -- one for the i386, and another for the x86_64.
This shouldn't be marked as "verified", since that means a *fix* has been verified. Putting to "assigned".
I can confirm this bug on Fedora 8, glibc-2.7-2, x86_64 architecture. $ gcc -lm -o bug-43358{,.c} $ ./bug-43358 4660.337778 * 8173.795556 + -38092648.215388 = 0 True result rounds to -3.52653456e-09
this bug is still present in Fedora 9. glibc-2.8-8 $ gcc -lm -o bug-43358{,.c} $ ./bug-43358 4660.337778 * 8173.795556 + -38092648.215388 = 0 True result rounds to -3.52653456e-09
this bug is still present in glibc-2.8.90-12.x86_64
The compiler now expands fma inline when optimizing.
glibc-2.12.90-16 has been submitted as an update for Fedora 14. https://admin.fedoraproject.org/updates/glibc-2.12.90-16
glibc-2.12.90-16 has been pushed to the Fedora 14 testing repository. If problems still persist, please make note of it in this bug report. If you want to test the update, you can install it with su -c 'yum --enablerepo=updates-testing update glibc'. You can provide feedback for this update here: https://admin.fedoraproject.org/updates/glibc-2.12.90-16
Confirmed fix with glibc-2.12.90-16.fc14.x86_64.
glibc-2.12.90-17 has been pushed to the Fedora 14 stable repository. If problems still persist, please make note of it in this bug report.