* (patch,committed) libquadmath: Update math/fmaq.c
@ 2012-11-15 17:25 Tobias Burnus
0 siblings, 0 replies; only message in thread
From: Tobias Burnus @ 2012-11-15 17:25 UTC (permalink / raw)
To: gcc patches
[-- Attachment #1: Type: text/plain, Size: 167 bytes --]
Dear all,
I have committed (Rev. 193538) attached patch, which does an other
update from GLIBC.
Tobias
PS: I still want to update libquadmath's strtod and printf.
[-- Attachment #2: committed.diff --]
[-- Type: text/x-patch, Size: 4084 bytes --]
Index: libquadmath/ChangeLog
===================================================================
--- libquadmath/ChangeLog (Revision 193537)
+++ libquadmath/ChangeLog (Arbeitskopie)
@@ -1,3 +1,11 @@
+2012-11-15 Tobias Burnus <burnus@net-b.de>
+ Joseph Myers <joseph@codesourcery.com>
+
+ * math/fmaq.c (fmaq): Merge from GLIBC. Fix fma
+ underflows with small x * y; Fix overflow results
+ outside round-to-nearest mode; make use of Dekker
+ and Knuth algorithms use round-to-nearest.
+
2012-11-01 Tobias Burnus <burnus@net-b.de>
* math/fmaq.c (fmaq): Fix build.
Index: libquadmath/math/fmaq.c
===================================================================
--- libquadmath/math/fmaq.c (Revision 193537)
+++ libquadmath/math/fmaq.c (Arbeitskopie)
@@ -14,9 +14,8 @@
Lesser General Public License for more details.
You should have received a copy of the GNU Lesser General Public
- License along with the GNU C Library; if not, write to the Free
- Software Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA
- 02111-1307 USA. */
+ License along with the GNU C Library; if not, see
+ <http://www.gnu.org/licenses/>. */
#include "quadmath-imp.h"
#include <math.h>
@@ -62,17 +61,18 @@ fmaq (__float128 x, __float128 y, __float128 z)
underflows to 0. */
if (z == 0 && x != 0 && y != 0)
return x * y;
- /* If x or y or z is Inf/NaN, or if fma will certainly overflow,
- or if x * y is less than half of FLT128_DENORM_MIN,
- compute as x * y + z. */
+ /* If x or y or z is Inf/NaN, or if x * y is zero, compute as
+ x * y + z. */
if (u.ieee.exponent == 0x7fff
|| v.ieee.exponent == 0x7fff
|| w.ieee.exponent == 0x7fff
- || u.ieee.exponent + v.ieee.exponent
- > 0x7fff + IEEE854_FLOAT128_BIAS
- || u.ieee.exponent + v.ieee.exponent
- < IEEE854_FLOAT128_BIAS - FLT128_MANT_DIG - 2)
+ || x == 0
+ || y == 0)
return x * y + z;
+ /* If fma will certainly overflow, compute as x * y. */
+ if (u.ieee.exponent + v.ieee.exponent
+ > 0x7fff + IEEE854_FLOAT128_BIAS)
+ return x * y;
/* If x * y is less than 1/4 of FLT128_DENORM_MIN, neither the
result nor whether there is underflow depends on its exact
value, only on its sign. */
@@ -121,9 +121,18 @@ fmaq (__float128 x, __float128 y, __float128 z)
{
/* Similarly.
If z exponent is very large and x and y exponents are
- very small, it doesn't matter if we don't adjust it. */
- if (u.ieee.exponent > v.ieee.exponent)
+ very small, adjust them up to avoid spurious underflows,
+ rather than down. */
+ if (u.ieee.exponent + v.ieee.exponent
+ <= IEEE854_FLOAT128_BIAS + FLT128_MANT_DIG)
{
+ if (u.ieee.exponent > v.ieee.exponent)
+ u.ieee.exponent += 2 * FLT128_MANT_DIG + 2;
+ else
+ v.ieee.exponent += 2 * FLT128_MANT_DIG + 2;
+ }
+ else if (u.ieee.exponent > v.ieee.exponent)
+ {
if (u.ieee.exponent > FLT128_MANT_DIG)
u.ieee.exponent -= FLT128_MANT_DIG;
}
@@ -175,6 +184,12 @@ fmaq (__float128 x, __float128 y, __float128 z)
if (__builtin_expect ((x == 0 || y == 0) && z == 0, 0))
return x * y + z;
+#ifdef USE_FENV_H
+ fenv_t env;
+ feholdexcept (&env);
+ fesetround (FE_TONEAREST);
+#endif
+
/* Multiplication m1 + m2 = x * y using Dekker's algorithm. */
#define C ((1LL << (FLT128_MANT_DIG + 1) / 2) + 1)
__float128 x1 = x * C;
@@ -193,10 +208,25 @@ fmaq (__float128 x, __float128 y, __float128 z)
t1 = m1 - t1;
t2 = z - t2;
__float128 a2 = t1 + t2;
+#ifdef USE_FENV_H
+ feclearexcept (FE_INEXACT);
+#endif
+ /* If the result is an exact zero, ensure it has the correct
+ sign. */
+ if (a1 == 0 && m2 == 0)
+ {
#ifdef USE_FENV_H
- fenv_t env;
- feholdexcept (&env);
+ feupdateenv (&env);
+#endif
+ /* Ensure that round-to-nearest value of z + m1 is not
+ reused. */
+ asm volatile ("" : "=m" (z) : "m" (z));
+ return z + m1;
+ }
+
+
+#ifdef USE_FENV_H
fesetround (FE_TOWARDZERO);
#endif
/* Perform m2 + a2 addition with round to odd. */
^ permalink raw reply [flat|nested] only message in thread
only message in thread, other threads:[~2012-11-15 17:25 UTC | newest]
Thread overview: (only message) (download: mbox.gz / follow: Atom feed)
-- links below jump to the message on this page --
2012-11-15 17:25 (patch,committed) libquadmath: Update math/fmaq.c Tobias Burnus
This is a public inbox, see mirroring instructions
for how to clone and mirror all data and code used for this inbox;
as well as URLs for read-only IMAP folder(s) and NNTP newsgroup(s).