about summary refs log tree commit diff
path: root/sysdeps/ieee754/dbl-64/s_fma.c
diff options
context:
space:
mode:
authorJoseph Myers <joseph@codesourcery.com>2012-11-03 19:48:53 +0000
committerJoseph Myers <joseph@codesourcery.com>2012-11-03 19:48:53 +0000
commit5b5b04d6282df0364424c6f2c0462e5c1a4394b0 (patch)
treebef8cb91fffbf78e56f18479234abb47ce3054b2 /sysdeps/ieee754/dbl-64/s_fma.c
parentfbeafedeea37e0af1984a6511018d159f5ceed6a (diff)
downloadglibc-5b5b04d6282df0364424c6f2c0462e5c1a4394b0.tar.gz
glibc-5b5b04d6282df0364424c6f2c0462e5c1a4394b0.tar.xz
glibc-5b5b04d6282df0364424c6f2c0462e5c1a4394b0.zip
Make fma use of Dekker and Knuth algorithms use round-to-nearest (bug 14796).
Diffstat (limited to 'sysdeps/ieee754/dbl-64/s_fma.c')
-rw-r--r--sysdeps/ieee754/dbl-64/s_fma.c18
1 files changed, 16 insertions, 2 deletions
diff --git a/sysdeps/ieee754/dbl-64/s_fma.c b/sysdeps/ieee754/dbl-64/s_fma.c
index 68c63a1987..07fe715617 100644
--- a/sysdeps/ieee754/dbl-64/s_fma.c
+++ b/sysdeps/ieee754/dbl-64/s_fma.c
@@ -167,6 +167,9 @@ __fma (double x, double y, double z)
   if (__builtin_expect ((x == 0 || y == 0) && z == 0, 0))
     return x * y + z;
 
+  fenv_t env;
+  libc_feholdexcept_setround (&env, FE_TONEAREST);
+
   /* Multiplication m1 + m2 = x * y using Dekker's algorithm.  */
 #define C ((1 << (DBL_MANT_DIG + 1) / 2) + 1)
   double x1 = x * C;
@@ -185,9 +188,20 @@ __fma (double x, double y, double z)
   t1 = m1 - t1;
   t2 = z - t2;
   double a2 = t1 + t2;
+  feclearexcept (FE_INEXACT);
 
-  fenv_t env;
-  libc_feholdexcept_setround (&env, FE_TOWARDZERO);
+  /* If the result is an exact zero, ensure it has the correct
+     sign.  */
+  if (a1 == 0 && m2 == 0)
+    {
+      libc_feupdateenv (&env);
+      /* Ensure that round-to-nearest value of z + m1 is not
+	 reused.  */
+      asm volatile ("" : "=m" (z) : "m" (z));
+      return z + m1;
+    }
+
+  libc_fesetround (FE_TOWARDZERO);
 
   /* Perform m2 + a2 addition with round to odd.  */
   u.d = a2 + m2;