about summary refs log tree commit diff
path: root/sysdeps/ieee754
diff options
context:
space:
mode:
authorJoseph Myers <joseph@codesourcery.com>2012-10-31 13:01:17 +0000
committerJoseph Myers <joseph@codesourcery.com>2012-10-31 13:01:17 +0000
commitef82f4da79332e504e17729536d501add283719b (patch)
tree02fa059ea66f20f529767cedbcb25a8be3d2fd8b /sysdeps/ieee754
parent954ef0d98d1fd8cc7a16a7289f286e885db36690 (diff)
downloadglibc-ef82f4da79332e504e17729536d501add283719b.tar.gz
glibc-ef82f4da79332e504e17729536d501add283719b.tar.xz
glibc-ef82f4da79332e504e17729536d501add283719b.zip
Fix fma underflow exceptions in after-rounding edge cases.
Diffstat (limited to 'sysdeps/ieee754')
-rw-r--r--sysdeps/ieee754/dbl-64/s_fma.c12
-rw-r--r--sysdeps/ieee754/ldbl-128/s_fmal.c12
-rw-r--r--sysdeps/ieee754/ldbl-96/s_fmal.c12
3 files changed, 36 insertions, 0 deletions
diff --git a/sysdeps/ieee754/dbl-64/s_fma.c b/sysdeps/ieee754/dbl-64/s_fma.c
index 108fd66192..155a773e36 100644
--- a/sysdeps/ieee754/dbl-64/s_fma.c
+++ b/sysdeps/ieee754/dbl-64/s_fma.c
@@ -22,6 +22,7 @@
 #include <fenv.h>
 #include <ieee754.h>
 #include <math_private.h>
+#include <tininess.h>
 
 /* This implementation uses rounding to odd to avoid problems with
    double rounding.  See a paper by Boldo and Melquiond:
@@ -208,6 +209,17 @@ __fma (double x, double y, double z)
 	 for proper rounding.  */
       if (v.ieee.exponent == 106)
 	{
+	  /* If the exponent would be in the normal range when
+	     rounding to normal precision with unbounded exponent
+	     range, the exact result is known and spurious underflows
+	     must be avoided on systems detecting tininess after
+	     rounding.  */
+	  if (TININESS_AFTER_ROUNDING)
+	    {
+	      w.d = a1 + u.d;
+	      if (w.ieee.exponent == 107)
+		return w.d * 0x1p-106;
+	    }
 	  /* v.ieee.mantissa1 & 2 is LSB bit of the result before rounding,
 	     v.ieee.mantissa1 & 1 is the round bit and j is our sticky
 	     bit.  */
diff --git a/sysdeps/ieee754/ldbl-128/s_fmal.c b/sysdeps/ieee754/ldbl-128/s_fmal.c
index 7fb9856ef1..201bff94ca 100644
--- a/sysdeps/ieee754/ldbl-128/s_fmal.c
+++ b/sysdeps/ieee754/ldbl-128/s_fmal.c
@@ -22,6 +22,7 @@
 #include <fenv.h>
 #include <ieee754.h>
 #include <math_private.h>
+#include <tininess.h>
 
 /* This implementation uses rounding to odd to avoid problems with
    double rounding.  See a paper by Boldo and Melquiond:
@@ -209,6 +210,17 @@ __fmal (long double x, long double y, long double z)
 	 for proper rounding.  */
       if (v.ieee.exponent == 226)
 	{
+	  /* If the exponent would be in the normal range when
+	     rounding to normal precision with unbounded exponent
+	     range, the exact result is known and spurious underflows
+	     must be avoided on systems detecting tininess after
+	     rounding.  */
+	  if (TININESS_AFTER_ROUNDING)
+	    {
+	      w.d = a1 + u.d;
+	      if (w.ieee.exponent == 227)
+		return w.d * 0x1p-226L;
+	    }
 	  /* v.ieee.mantissa3 & 2 is LSB bit of the result before rounding,
 	     v.ieee.mantissa3 & 1 is the round bit and j is our sticky
 	     bit.  */
diff --git a/sysdeps/ieee754/ldbl-96/s_fmal.c b/sysdeps/ieee754/ldbl-96/s_fmal.c
index 73104914b3..91333af353 100644
--- a/sysdeps/ieee754/ldbl-96/s_fmal.c
+++ b/sysdeps/ieee754/ldbl-96/s_fmal.c
@@ -22,6 +22,7 @@
 #include <fenv.h>
 #include <ieee754.h>
 #include <math_private.h>
+#include <tininess.h>
 
 /* This implementation uses rounding to odd to avoid problems with
    double rounding.  See a paper by Boldo and Melquiond:
@@ -209,6 +210,17 @@ __fmal (long double x, long double y, long double z)
 	 for proper rounding.  */
       if (v.ieee.exponent == 128)
 	{
+	  /* If the exponent would be in the normal range when
+	     rounding to normal precision with unbounded exponent
+	     range, the exact result is known and spurious underflows
+	     must be avoided on systems detecting tininess after
+	     rounding.  */
+	  if (TININESS_AFTER_ROUNDING)
+	    {
+	      w.d = a1 + u.d;
+	      if (w.ieee.exponent == 129)
+		return w.d * 0x1p-128L;
+	    }
 	  /* v.ieee.mantissa1 & 2 is LSB bit of the result before rounding,
 	     v.ieee.mantissa1 & 1 is the round bit and j is our sticky
 	     bit.  */