diff options
| author | Joseph Myers <joseph@codesourcery.com> | 2012-10-31 13:01:17 +0000 |
|---|---|---|
| committer | Joseph Myers <joseph@codesourcery.com> | 2012-10-31 13:01:17 +0000 |
| commit | ef82f4da79332e504e17729536d501add283719b (patch) | |
| tree | 02fa059ea66f20f529767cedbcb25a8be3d2fd8b /sysdeps | |
| parent | 954ef0d98d1fd8cc7a16a7289f286e885db36690 (diff) | |
| download | glibc-ef82f4da79332e504e17729536d501add283719b.zip glibc-ef82f4da79332e504e17729536d501add283719b.tar.gz glibc-ef82f4da79332e504e17729536d501add283719b.tar.bz2 | |
Fix fma underflow exceptions in after-rounding edge cases.
Diffstat (limited to 'sysdeps')
| -rw-r--r-- | sysdeps/ieee754/dbl-64/s_fma.c | 12 | ||||
| -rw-r--r-- | sysdeps/ieee754/ldbl-128/s_fmal.c | 12 | ||||
| -rw-r--r-- | sysdeps/ieee754/ldbl-96/s_fmal.c | 12 |
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 108fd66..155a773 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 7fb9856..201bff9 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 7310491..91333af 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. */ |