Bug 43358

Summary: fma(x,y,z) does not compute x*y+z as a single ternary operation
Product: [Fedora] Fedora Reporter: Trevin Beattie <trevin>
Component: glibcAssignee: Andreas Schwab <schwab>
Status: CLOSED ERRATA QA Contact:
Severity: medium Docs Contact:
Priority: low    
Version: rawhideCC: fweimer, jakub, kparal, mattdm, notting, petrosyan
Target Milestone: ---Keywords: FutureFeature
Target Release: ---   
Hardware: All   
OS: Linux   
Whiteboard:
Fixed In Version: glibc-2.12.90-17 Doc Type: Enhancement
Doc Text:
Story Points: ---
Clone Of: Environment:
Last Closed: 2010-10-19 22:22:42 UTC Type: ---
Regression: --- Mount Type: ---
Documentation: --- CRM:
Verified Versions: Category: ---
oVirt Team: --- RHEL 7.3 requirements from Atomic Host:
Cloudforms Team: --- Target Upstream Version:
Embargoed:
Attachments:
Description Flags
Fused multiply and add - long double version for i386, written for MASM
none
Test cases for fmal() (i386-specific). If anyone has ideas for additional tests, let me know. The above assembly file passes all these tests. none

Description Trevin Beattie 2001-06-03 20:11:59 UTC
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.

Comment 1 Jakub Jelinek 2001-06-05 08:06:48 UTC
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.

Comment 2 Trevin Beattie 2001-06-08 00:53:26 UTC
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.


Comment 3 Trevin Beattie 2001-06-08 00:55:13 UTC
Created attachment 20593 [details]
Fused multiply and add - long double version for i386, written for MASM

Comment 4 Trevin Beattie 2001-06-08 00:58:57 UTC
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.

Comment 5 Jakub Jelinek 2001-06-08 08:28:02 UTC
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.

Comment 6 Bill Nottingham 2006-08-07 18:31:17 UTC
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.


Comment 7 Trevin Beattie 2006-08-08 16:49:15 UTC
I've verified that the bug is still present in Fedora Core 6 RC2 (glibc-2.4.90-17).

Comment 8 Matthew Miller 2007-04-06 16:11:47 UTC
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.]


Comment 9 Trevin Beattie 2007-04-07 02:23:23 UTC
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.


Comment 10 Matthew Miller 2007-04-07 03:55:00 UTC
This shouldn't be marked as "verified", since that means a *fix* has been
verified. Putting to "assigned".

Comment 11 petrosyan 2008-02-14 00:57:10 UTC
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

Comment 12 petrosyan 2008-08-15 05:28:44 UTC
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

Comment 13 petrosyan 2008-10-01 17:51:11 UTC
this bug is still present in glibc-2.8.90-12.x86_64

Comment 14 Andreas Schwab 2009-06-25 14:40:42 UTC
The compiler now expands fma inline when optimizing.

Comment 15 Fedora Update System 2010-10-14 12:03:47 UTC
glibc-2.12.90-16 has been submitted as an update for Fedora 14.
https://admin.fedoraproject.org/updates/glibc-2.12.90-16

Comment 16 Fedora Update System 2010-10-14 13:51:38 UTC
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

Comment 17 Kamil Páral 2010-10-19 17:44:37 UTC
Confirmed fix with glibc-2.12.90-16.fc14.x86_64.

Comment 18 Fedora Update System 2010-10-19 22:22:20 UTC
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.