public inbox for libc-hacker@sourceware.org
 help / color / mirror / Atom feed
* [PATCH] Implement fmal, some fma bugfixes
@ 2010-10-15 19:08 Jakub Jelinek
  0 siblings, 0 replies; only message in thread
From: Jakub Jelinek @ 2010-10-15 19:08 UTC (permalink / raw)
  To: Ulrich Drepper; +Cc: Glibc hackers

[-- Attachment #1: Type: text/plain, Size: 26568 bytes --]

Hi!

This patch on top of
http://sources.redhat.com/ml/libc-alpha/2010-10/msg00022.html
implements extended long double and IEEE quad long double
fmal and fixes some bugs in fma, in particular
fma (DBL_MAX, DBL_MAX, -__builtin_inf ()) should be
__builtin_inf () rather than NaN, as the computation is supposed
to happen with infinite precision before rounding and thus x * y
will be still finite, albeit very large.
Also, unfortunately it seems two fetestexcept calls for adjust == -1
are unavoidable as shown by one of the added testcases -
if result is slightly above DBL_MIN, we compute (a1 + u.d) * 0x1p106
in round to nearest, but if the inexact bit from first addition
hasn't been ored into u.d, it is lost.

Attached is also greatly improved mpfr testing proglet, now handles
all of float (-DDO_FLT), double (-DDO_DBL) and long double (-DDO_LDBL,
both extended and quad).  Tested on x86_64-linux and i686-linux and
slightly on s390x-linux, unfortunately I've discovered mpfr_set_ld
is buggy there, so the tester is currently unusable there.

2010-10-15  Jakub Jelinek  <jakub@redhat.com>

	[BZ #3268]
 	* math/libm-test.inc (fma_test): Some new testcases.
	* sysdeps/ieee754/ldbl-128/s_fmal.c: New file.
	* sysdeps/ieee754/ldbl-96/s_fma.c (__fma): Fix fma with finite x and
	y and infinite z.  Do multiplication by C already in long double.
	* sysdeps/ieee754/ldbl-96/s_fmal.c: New file.
	* sysdeps/ieee754/dbl-64/s_fma.c (__fma): Fix fma with finite x and
	y and infinite z.  Do bitwise or of inexact bit into u.d.
	* sysdeps/ieee754/ldbl-64-128/s_fmal.c: New file.
	* sysdeps/i386/fpu/s_fmaf.S: Removed.
	* sysdeps/i386/fpu/s_fma.S: Removed.
	* sysdeps/i386/fpu/s_fmal.S: Removed.

--- libc/math/libm-test.inc.jj	2010-10-14 20:21:24.000000000 +0200
+++ libc/math/libm-test.inc	2010-10-15 17:51:16.000000000 +0200
@@ -2787,8 +2787,24 @@ fma_test (void)
   TEST_fff_f (fma, minus_infty, plus_infty, plus_infty, nan_value, INVALID_EXCEPTION);
   TEST_fff_f (fma, plus_infty, minus_infty, plus_infty, nan_value, INVALID_EXCEPTION);
   TEST_fff_f (fma, minus_infty, minus_infty, minus_infty, nan_value, INVALID_EXCEPTION);
+  TEST_fff_f (fma, plus_infty, 3.5L, minus_infty, nan_value, INVALID_EXCEPTION);
+  TEST_fff_f (fma, minus_infty, -7.5L, minus_infty, nan_value, INVALID_EXCEPTION);
+  TEST_fff_f (fma, -13.5L, plus_infty, plus_infty, nan_value, INVALID_EXCEPTION);
+  TEST_fff_f (fma, minus_infty, 7.5L, plus_infty, nan_value, INVALID_EXCEPTION);
 
   TEST_fff_f (fma, 1.25L, 0.75L, 0.0625L, 1.0L);
+
+  FLOAT fltmax = CHOOSE (LDBL_MAX, DBL_MAX, FLT_MAX,
+			 LDBL_MAX, DBL_MAX, FLT_MAX);
+  TEST_fff_f (fma, -fltmax, -fltmax, minus_infty, minus_infty);
+  TEST_fff_f (fma, fltmax / 2, fltmax / 2, minus_infty, minus_infty);
+  TEST_fff_f (fma, -fltmax, fltmax, plus_infty, plus_infty);
+  TEST_fff_f (fma, fltmax / 2, -fltmax / 4, plus_infty, plus_infty);
+  TEST_fff_f (fma, plus_infty, 4, plus_infty, plus_infty);
+  TEST_fff_f (fma, 2, minus_infty, minus_infty, minus_infty);
+  TEST_fff_f (fma, minus_infty, minus_infty, plus_infty, plus_infty);
+  TEST_fff_f (fma, plus_infty, minus_infty, minus_infty, minus_infty);
+
 #if defined (TEST_FLOAT) && FLT_MANT_DIG == 24
   TEST_fff_f (fma, 0x1.7ff8p+13, 0x1.000002p+0, 0x1.ffffp-24, 0x1.7ff802p+13);
   TEST_fff_f (fma, 0x1.fffp+0, 0x1.00001p+0, -0x1.fffp+0, 0x1.fffp-20);
@@ -2818,6 +2834,15 @@ fma_test (void)
   TEST_fff_f (fma, -0x1.19cab66d73e17p-959, 0x1.c7108a8c5ff51p-107, -0x0.80b0ad65d9b64p-1022, -0x0.80b0ad65d9d59p-1022);
   TEST_fff_f (fma, -0x1.d2eaed6e8e9d3p-979, -0x1.4e066c62ac9ddp-63, -0x0.9245e6b003454p-1022, -0x0.9245c09c5fb5dp-1022);
   TEST_fff_f (fma, 0x1.153d650bb9f06p-907, 0x1.2d01230d48407p-125, -0x0.b278d5acfc3cp-1022, -0x0.b22757123bbe9p-1022);
+  TEST_fff_f (fma, -0x1.fffffffffffffp-711, 0x1.fffffffffffffp-275, 0x1.fffffe00007ffp-983, 0x1.7ffffe00007ffp-983);
+#endif
+#if defined (TEST_LDOUBLE) && LDBL_MANT_DIG == 64
+  TEST_fff_f (fma, -0x8.03fcp+3696L, 0xf.fffffffffffffffp-6140L, 0x8.3ffffffffffffffp-2450L, -0x8.01ecp-2440L);
+  TEST_fff_f (fma, 0x9.fcp+2033L, -0x8.000e1f000ff800fp-3613L, -0xf.fffffffffffc0ffp-1579L, -0xd.fc119fb093ed092p-1577L);
+  TEST_fff_f (fma, 0xc.7fc000003ffffffp-1194L, 0x8.1e0003fffffffffp+15327L, -0x8.fffep+14072L, 0xc.ae9f164020effffp+14136L);
+  TEST_fff_f (fma, -0x8.0001fc000000003p+1798L, 0xcp-2230L, 0x8.f7e000000000007p-468L, -0xc.0002f9ffee10404p-429L);
+  TEST_fff_f (fma, 0xc.0000000000007ffp+10130L, -0x8.000000000000001p+4430L, 0xc.07000000001ffffp+14513L, -0xb.fffffffffffd7e4p+14563L);
+  TEST_fff_f (fma, 0xb.ffffp-4777L, 0x8.000000fffffffffp-11612L, -0x0.3800fff8p-16385L, 0x5.c7fe80c7ffeffffp-16385L);
 #endif
 
   END (fma);
--- libc/sysdeps/ieee754/ldbl-128/s_fmal.c.jj	2010-10-15 18:03:31.000000000 +0200
+++ libc/sysdeps/ieee754/ldbl-128/s_fmal.c	2010-10-15 18:08:43.000000000 +0200
@@ -0,0 +1,221 @@
+/* Compute x * y + z as ternary operation.
+   Copyright (C) 2010 Free Software Foundation, Inc.
+   This file is part of the GNU C Library.
+   Contributed by Jakub Jelinek <jakub@redhat.com>, 2010.
+
+   The GNU C Library is free software; you can redistribute it and/or
+   modify it under the terms of the GNU Lesser General Public
+   License as published by the Free Software Foundation; either
+   version 2.1 of the License, or (at your option) any later version.
+
+   The GNU C Library is distributed in the hope that it will be useful,
+   but WITHOUT ANY WARRANTY; without even the implied warranty of
+   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+   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.  */
+
+#include <float.h>
+#include <math.h>
+#include <fenv.h>
+#include <ieee754.h>
+
+/* This implementation uses rounding to odd to avoid problems with
+   double rounding.  See a paper by Boldo and Melquiond:
+   http://www.lri.fr/~melquion/doc/08-tc.pdf  */
+
+long double
+__fmal (long double x, long double y, long double z)
+{
+  union ieee854_long_double u, v, w;
+  int adjust = 0;
+  u.d = x;
+  v.d = y;
+  w.d = z;
+  if (__builtin_expect (u.ieee.exponent + v.ieee.exponent
+			>= 0x7fff + IEEE854_LONG_DOUBLE_BIAS
+			   - LDBL_MANT_DIG, 0)
+      || __builtin_expect (u.ieee.exponent >= 0x7fff - LDBL_MANT_DIG, 0)
+      || __builtin_expect (v.ieee.exponent >= 0x7fff - LDBL_MANT_DIG, 0)
+      || __builtin_expect (w.ieee.exponent >= 0x7fff - LDBL_MANT_DIG, 0)
+      || __builtin_expect (u.ieee.exponent + v.ieee.exponent
+			   <= IEEE854_LONG_DOUBLE_BIAS + LDBL_MANT_DIG, 0))
+    {
+      /* If z is Inf, but x and y are finite, the result should be
+	 z rather than NaN.  */
+      if (w.ieee.exponent == 0x7fff
+	  && u.ieee.exponent != 0x7fff
+          && v.ieee.exponent != 0x7fff)
+	return (z + 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 LDBL_DENORM_MIN,
+	 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_LONG_DOUBLE_BIAS
+	  || u.ieee.exponent + v.ieee.exponent
+	     < IEEE854_LONG_DOUBLE_BIAS - LDBL_MANT_DIG - 2)
+	return x * y + z;
+      if (u.ieee.exponent + v.ieee.exponent
+	  >= 0x7fff + IEEE854_LONG_DOUBLE_BIAS - LDBL_MANT_DIG)
+	{
+	  /* Compute 1p-113 times smaller result and multiply
+	     at the end.  */
+	  if (u.ieee.exponent > v.ieee.exponent)
+	    u.ieee.exponent -= LDBL_MANT_DIG;
+	  else
+	    v.ieee.exponent -= LDBL_MANT_DIG;
+	  /* If x + y exponent is very large and z exponent is very small,
+	     it doesn't matter if we don't adjust it.  */
+	  if (w.ieee.exponent > LDBL_MANT_DIG)
+	    w.ieee.exponent -= LDBL_MANT_DIG;
+	  adjust = 1;
+	}
+      else if (w.ieee.exponent >= 0x7fff - LDBL_MANT_DIG)
+	{
+	  /* 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)
+	    {
+	      if (u.ieee.exponent > LDBL_MANT_DIG)
+		u.ieee.exponent -= LDBL_MANT_DIG;
+	    }
+	  else if (v.ieee.exponent > LDBL_MANT_DIG)
+	    v.ieee.exponent -= LDBL_MANT_DIG;
+	  w.ieee.exponent -= LDBL_MANT_DIG;
+	  adjust = 1;
+	}
+      else if (u.ieee.exponent >= 0x7fff - LDBL_MANT_DIG)
+	{
+	  u.ieee.exponent -= LDBL_MANT_DIG;
+	  if (v.ieee.exponent)
+	    v.ieee.exponent += LDBL_MANT_DIG;
+	  else
+	    v.d *= 0x1p113L;
+	}
+      else if (v.ieee.exponent >= 0x7fff - LDBL_MANT_DIG)
+	{
+	  v.ieee.exponent -= LDBL_MANT_DIG;
+	  if (u.ieee.exponent)
+	    u.ieee.exponent += LDBL_MANT_DIG;
+	  else
+	    u.d *= 0x1p113L;
+	}
+      else /* if (u.ieee.exponent + v.ieee.exponent
+		  <= IEEE854_LONG_DOUBLE_BIAS + LDBL_MANT_DIG) */
+	{
+	  if (u.ieee.exponent > v.ieee.exponent)
+	    u.ieee.exponent += 2 * LDBL_MANT_DIG;
+	  else
+	    v.ieee.exponent += 2 * LDBL_MANT_DIG;
+	  if (w.ieee.exponent <= 4 * LDBL_MANT_DIG + 4)
+	    {
+	      if (w.ieee.exponent)
+		w.ieee.exponent += 2 * LDBL_MANT_DIG;
+	      else
+		w.d *= 0x1p226L;
+	      adjust = -1;
+	    }
+	  /* Otherwise x * y should just affect inexact
+	     and nothing else.  */
+	}
+      x = u.d;
+      y = v.d;
+      z = w.d;
+    }
+  /* Multiplication m1 + m2 = x * y using Dekker's algorithm.  */
+#define C ((1LL << (LDBL_MANT_DIG + 1) / 2) + 1)
+  long double x1 = x * C;
+  long double y1 = y * C;
+  long double m1 = x * y;
+  x1 = (x - x1) + x1;
+  y1 = (y - y1) + y1;
+  long double x2 = x - x1;
+  long double y2 = y - y1;
+  long double m2 = (((x1 * y1 - m1) + x1 * y2) + x2 * y1) + x2 * y2;
+
+  /* Addition a1 + a2 = z + m1 using Knuth's algorithm.  */
+  long double a1 = z + m1;
+  long double t1 = a1 - z;
+  long double t2 = a1 - t1;
+  t1 = m1 - t1;
+  t2 = z - t2;
+  long double a2 = t1 + t2;
+
+  fenv_t env;
+  feholdexcept (&env);
+  fesetround (FE_TOWARDZERO);
+  /* Perform m2 + a2 addition with round to odd.  */
+  u.d = a2 + m2;
+
+  if (__builtin_expect (adjust == 0, 1))
+    {
+      if ((u.ieee.mantissa3 & 1) == 0 && u.ieee.exponent != 0x7fff)
+	u.ieee.mantissa3 |= fetestexcept (FE_INEXACT) != 0;
+      feupdateenv (&env);
+      /* Result is a1 + u.d.  */
+      return a1 + u.d;
+    }
+  else if (__builtin_expect (adjust > 0, 1))
+    {
+      if ((u.ieee.mantissa3 & 1) == 0 && u.ieee.exponent != 0x7fff)
+	u.ieee.mantissa3 |= fetestexcept (FE_INEXACT) != 0;
+      feupdateenv (&env);
+      /* Result is a1 + u.d, scaled up.  */
+      return (a1 + u.d) * 0x1p113L;
+    }
+  else
+    {
+      if ((u.ieee.mantissa3 & 1) == 0)
+	u.ieee.mantissa3 |= fetestexcept (FE_INEXACT) != 0;
+      v.d = a1 + u.d;
+      int j = fetestexcept (FE_INEXACT) != 0;
+      feupdateenv (&env);
+      /* Ensure the following computations are performed in default rounding
+	 mode instead of just reusing the round to zero computation.  */
+      asm volatile ("" : "=m" (u) : "m" (u));
+      /* If a1 + u.d is exact, the only rounding happens during
+	 scaling down.  */
+      if (j == 0)
+	return v.d * 0x1p-226L;
+      /* If result rounded to zero is not subnormal, no double
+	 rounding will occur.  */
+      if (v.ieee.exponent > 226)
+	return (a1 + u.d) * 0x1p-226L;
+      /* If v.d * 0x1p-226L with round to zero is a subnormal above
+	 or equal to LDBL_MIN / 2, then v.d * 0x1p-226L shifts mantissa
+	 down just by 1 bit, which means v.ieee.mantissa3 |= j would
+	 change the round bit, not sticky or guard bit.
+	 v.d * 0x1p-226L never normalizes by shifting up,
+	 so round bit plus sticky bit should be already enough
+	 for proper rounding.  */
+      if (v.ieee.exponent == 226)
+	{
+	  /* 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.  In round-to-nearest 001 rounds down like 00,
+	     011 rounds up, even though 01 rounds down (thus we need
+	     to adjust), 101 rounds down like 10 and 111 rounds up
+	     like 11.  */
+	  if ((v.ieee.mantissa3 & 3) == 1)
+	    {
+	      v.d *= 0x1p-226L;
+	      if (v.ieee.negative)
+		return v.d - 0x1p-16493L /* __LDBL_DENORM_MIN__ */;
+	      else
+		return v.d + 0x1p-16493L /* __LDBL_DENORM_MIN__ */;
+	    }
+	  else
+	    return v.d * 0x1p-226L;
+	}
+      v.ieee.mantissa3 |= j;
+      return v.d * 0x1p-226L;
+    }
+}
+weak_alias (__fmal, fmal)
--- libc/sysdeps/ieee754/ldbl-96/s_fma.c.jj	2010-10-14 09:08:58.000000000 +0200
+++ libc/sysdeps/ieee754/ldbl-96/s_fma.c	2010-10-15 13:41:15.000000000 +0200
@@ -30,11 +30,20 @@
 double
 __fma (double x, double y, double z)
 {
+  if (__builtin_expect (isinf (z), 0))
+    {
+      /* If z is Inf, but x and y are finite, the result should be
+	 z rather than NaN.  */
+      if (finite (x) && finite (y))
+	return (z + x) + y;
+      return (x * y) + z;
+    }
+
   /* Multiplication m1 + m2 = x * y using Dekker's algorithm.  */
 #define C ((1ULL << (LDBL_MANT_DIG + 1) / 2) + 1)
-  long double x1 = x * C;
-  long double y1 = y * C;
-  long double m1 = x * y;
+  long double x1 = (long double) x * C;
+  long double y1 = (long double) y * C;
+  long double m1 = (long double) x * y;
   x1 = (x - x1) + x1;
   y1 = (y - y1) + y1;
   long double x2 = x - x1;
--- libc/sysdeps/ieee754/ldbl-96/s_fmal.c.jj	2010-10-14 22:39:23.000000000 +0200
+++ libc/sysdeps/ieee754/ldbl-96/s_fmal.c	2010-10-15 17:51:46.000000000 +0200
@@ -0,0 +1,221 @@
+/* Compute x * y + z as ternary operation.
+   Copyright (C) 2010 Free Software Foundation, Inc.
+   This file is part of the GNU C Library.
+   Contributed by Jakub Jelinek <jakub@redhat.com>, 2010.
+
+   The GNU C Library is free software; you can redistribute it and/or
+   modify it under the terms of the GNU Lesser General Public
+   License as published by the Free Software Foundation; either
+   version 2.1 of the License, or (at your option) any later version.
+
+   The GNU C Library is distributed in the hope that it will be useful,
+   but WITHOUT ANY WARRANTY; without even the implied warranty of
+   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+   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.  */
+
+#include <float.h>
+#include <math.h>
+#include <fenv.h>
+#include <ieee754.h>
+
+/* This implementation uses rounding to odd to avoid problems with
+   double rounding.  See a paper by Boldo and Melquiond:
+   http://www.lri.fr/~melquion/doc/08-tc.pdf  */
+
+long double
+__fmal (long double x, long double y, long double z)
+{
+  union ieee854_long_double u, v, w;
+  int adjust = 0;
+  u.d = x;
+  v.d = y;
+  w.d = z;
+  if (__builtin_expect (u.ieee.exponent + v.ieee.exponent
+			>= 0x7fff + IEEE854_LONG_DOUBLE_BIAS
+			   - LDBL_MANT_DIG, 0)
+      || __builtin_expect (u.ieee.exponent >= 0x7fff - LDBL_MANT_DIG, 0)
+      || __builtin_expect (v.ieee.exponent >= 0x7fff - LDBL_MANT_DIG, 0)
+      || __builtin_expect (w.ieee.exponent >= 0x7fff - LDBL_MANT_DIG, 0)
+      || __builtin_expect (u.ieee.exponent + v.ieee.exponent
+			   <= IEEE854_LONG_DOUBLE_BIAS + LDBL_MANT_DIG, 0))
+    {
+      /* If z is Inf, but x and y are finite, the result should be
+	 z rather than NaN.  */
+      if (w.ieee.exponent == 0x7fff
+	  && u.ieee.exponent != 0x7fff
+          && v.ieee.exponent != 0x7fff)
+	return (z + 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 LDBL_DENORM_MIN,
+	 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_LONG_DOUBLE_BIAS
+	  || u.ieee.exponent + v.ieee.exponent
+	     < IEEE854_LONG_DOUBLE_BIAS - LDBL_MANT_DIG - 2)
+	return x * y + z;
+      if (u.ieee.exponent + v.ieee.exponent
+	  >= 0x7fff + IEEE854_LONG_DOUBLE_BIAS - LDBL_MANT_DIG)
+	{
+	  /* Compute 1p-64 times smaller result and multiply
+	     at the end.  */
+	  if (u.ieee.exponent > v.ieee.exponent)
+	    u.ieee.exponent -= LDBL_MANT_DIG;
+	  else
+	    v.ieee.exponent -= LDBL_MANT_DIG;
+	  /* If x + y exponent is very large and z exponent is very small,
+	     it doesn't matter if we don't adjust it.  */
+	  if (w.ieee.exponent > LDBL_MANT_DIG)
+	    w.ieee.exponent -= LDBL_MANT_DIG;
+	  adjust = 1;
+	}
+      else if (w.ieee.exponent >= 0x7fff - LDBL_MANT_DIG)
+	{
+	  /* 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)
+	    {
+	      if (u.ieee.exponent > LDBL_MANT_DIG)
+		u.ieee.exponent -= LDBL_MANT_DIG;
+	    }
+	  else if (v.ieee.exponent > LDBL_MANT_DIG)
+	    v.ieee.exponent -= LDBL_MANT_DIG;
+	  w.ieee.exponent -= LDBL_MANT_DIG;
+	  adjust = 1;
+	}
+      else if (u.ieee.exponent >= 0x7fff - LDBL_MANT_DIG)
+	{
+	  u.ieee.exponent -= LDBL_MANT_DIG;
+	  if (v.ieee.exponent)
+	    v.ieee.exponent += LDBL_MANT_DIG;
+	  else
+	    v.d *= 0x1p64L;
+	}
+      else if (v.ieee.exponent >= 0x7fff - LDBL_MANT_DIG)
+	{
+	  v.ieee.exponent -= LDBL_MANT_DIG;
+	  if (u.ieee.exponent)
+	    u.ieee.exponent += LDBL_MANT_DIG;
+	  else
+	    u.d *= 0x1p64L;
+	}
+      else /* if (u.ieee.exponent + v.ieee.exponent
+		  <= IEEE854_LONG_DOUBLE_BIAS + LDBL_MANT_DIG) */
+	{
+	  if (u.ieee.exponent > v.ieee.exponent)
+	    u.ieee.exponent += 2 * LDBL_MANT_DIG;
+	  else
+	    v.ieee.exponent += 2 * LDBL_MANT_DIG;
+	  if (w.ieee.exponent <= 4 * LDBL_MANT_DIG + 4)
+	    {
+	      if (w.ieee.exponent)
+		w.ieee.exponent += 2 * LDBL_MANT_DIG;
+	      else
+		w.d *= 0x1p128L;
+	      adjust = -1;
+	    }
+	  /* Otherwise x * y should just affect inexact
+	     and nothing else.  */
+	}
+      x = u.d;
+      y = v.d;
+      z = w.d;
+    }
+  /* Multiplication m1 + m2 = x * y using Dekker's algorithm.  */
+#define C ((1LL << (LDBL_MANT_DIG + 1) / 2) + 1)
+  long double x1 = x * C;
+  long double y1 = y * C;
+  long double m1 = x * y;
+  x1 = (x - x1) + x1;
+  y1 = (y - y1) + y1;
+  long double x2 = x - x1;
+  long double y2 = y - y1;
+  long double m2 = (((x1 * y1 - m1) + x1 * y2) + x2 * y1) + x2 * y2;
+
+  /* Addition a1 + a2 = z + m1 using Knuth's algorithm.  */
+  long double a1 = z + m1;
+  long double t1 = a1 - z;
+  long double t2 = a1 - t1;
+  t1 = m1 - t1;
+  t2 = z - t2;
+  long double a2 = t1 + t2;
+
+  fenv_t env;
+  feholdexcept (&env);
+  fesetround (FE_TOWARDZERO);
+  /* Perform m2 + a2 addition with round to odd.  */
+  u.d = a2 + m2;
+
+  if (__builtin_expect (adjust == 0, 1))
+    {
+      if ((u.ieee.mantissa1 & 1) == 0 && u.ieee.exponent != 0x7fff)
+	u.ieee.mantissa1 |= fetestexcept (FE_INEXACT) != 0;
+      feupdateenv (&env);
+      /* Result is a1 + u.d.  */
+      return a1 + u.d;
+    }
+  else if (__builtin_expect (adjust > 0, 1))
+    {
+      if ((u.ieee.mantissa1 & 1) == 0 && u.ieee.exponent != 0x7fff)
+	u.ieee.mantissa1 |= fetestexcept (FE_INEXACT) != 0;
+      feupdateenv (&env);
+      /* Result is a1 + u.d, scaled up.  */
+      return (a1 + u.d) * 0x1p64L;
+    }
+  else
+    {
+      if ((u.ieee.mantissa1 & 1) == 0)
+	u.ieee.mantissa1 |= fetestexcept (FE_INEXACT) != 0;
+      v.d = a1 + u.d;
+      int j = fetestexcept (FE_INEXACT) != 0;
+      feupdateenv (&env);
+      /* Ensure the following computations are performed in default rounding
+	 mode instead of just reusing the round to zero computation.  */
+      asm volatile ("" : "=m" (u) : "m" (u));
+      /* If a1 + u.d is exact, the only rounding happens during
+	 scaling down.  */
+      if (j == 0)
+	return v.d * 0x1p-128L;
+      /* If result rounded to zero is not subnormal, no double
+	 rounding will occur.  */
+      if (v.ieee.exponent > 128)
+	return (a1 + u.d) * 0x1p-128L;
+      /* If v.d * 0x1p-128L with round to zero is a subnormal above
+	 or equal to LDBL_MIN / 2, then v.d * 0x1p-128L shifts mantissa
+	 down just by 1 bit, which means v.ieee.mantissa1 |= j would
+	 change the round bit, not sticky or guard bit.
+	 v.d * 0x1p-128L never normalizes by shifting up,
+	 so round bit plus sticky bit should be already enough
+	 for proper rounding.  */
+      if (v.ieee.exponent == 128)
+	{
+	  /* 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.  In round-to-nearest 001 rounds down like 00,
+	     011 rounds up, even though 01 rounds down (thus we need
+	     to adjust), 101 rounds down like 10 and 111 rounds up
+	     like 11.  */
+	  if ((v.ieee.mantissa1 & 3) == 1)
+	    {
+	      v.d *= 0x1p-128L;
+	      if (v.ieee.negative)
+		return v.d - 0x1p-16445L /* __LDBL_DENORM_MIN__ */;
+	      else
+		return v.d + 0x1p-16445L /* __LDBL_DENORM_MIN__ */;
+	    }
+	  else
+	    return v.d * 0x1p-128L;
+	}
+      v.ieee.mantissa1 |= j;
+      return v.d * 0x1p-128L;
+    }
+}
+weak_alias (__fmal, fmal)
--- libc/sysdeps/ieee754/dbl-64/s_fma.c.jj	2010-10-14 20:26:05.000000000 +0200
+++ libc/sysdeps/ieee754/dbl-64/s_fma.c	2010-10-15 17:56:30.000000000 +0200
@@ -43,6 +43,12 @@ __fma (double x, double y, double z)
       || __builtin_expect (u.ieee.exponent + v.ieee.exponent
 			   <= IEEE754_DOUBLE_BIAS + DBL_MANT_DIG, 0))
     {
+      /* If z is Inf, but x and y are finite, the result should be
+	 z rather than NaN.  */
+      if (w.ieee.exponent == 0x7ff
+	  && u.ieee.exponent != 0x7ff
+          && v.ieee.exponent != 0x7ff)
+	return (z + 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 DBL_DENORM_MIN,
 	 compute as x * y + z.  */
@@ -165,6 +171,8 @@ __fma (double x, double y, double z)
     }
   else
     {
+      if ((u.ieee.mantissa1 & 1) == 0)
+	u.ieee.mantissa1 |= fetestexcept (FE_INEXACT) != 0;
       v.d = a1 + u.d;
       int j = fetestexcept (FE_INEXACT) != 0;
       feupdateenv (&env);
--- libc/sysdeps/ieee754/ldbl-64-128/s_fmal.c.jj	2010-10-15 18:05:04.000000000 +0200
+++ libc/sysdeps/ieee754/ldbl-64-128/s_fmal.c	2010-10-15 18:05:17.000000000 +0200
@@ -0,0 +1,5 @@
+#include <math_ldbl_opt.h>
+#undef weak_alias
+#define weak_alias(n,a)
+#include <sysdeps/ieee754/ldbl-128/s_fmal.c>
+long_double_symbol (libm, __fmal, fmal);
--- libc/sysdeps/i386/fpu/s_fmaf.S.jj	2009-05-16 19:23:38.000000000 +0200
+++ libc/sysdeps/i386/fpu/s_fmaf.S	2010-09-14 11:29:01.831482015 +0200
@@ -1,31 +0,0 @@
-/* Compute (X * Y) + Z as ternary operation.
-   Copyright (C) 1997 Free Software Foundation, Inc.
-   This file is part of the GNU C Library.
-   Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997.
-
-   The GNU C Library is free software; you can redistribute it and/or
-   modify it under the terms of the GNU Lesser General Public
-   License as published by the Free Software Foundation; either
-   version 2.1 of the License, or (at your option) any later version.
-
-   The GNU C Library is distributed in the hope that it will be useful,
-   but WITHOUT ANY WARRANTY; without even the implied warranty of
-   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
-   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.  */
-
-#include <sysdep.h>
-
-	.text
-ENTRY(__fmaf)
-	flds	4(%esp)		// x
-	fmuls	8(%esp)		// x * y
-	flds	12(%esp)	// z : x * y
-	faddp			// (x * y) + z
-	ret
-END(__fmaf)
-weak_alias (__fmaf, fmaf)
--- libc/sysdeps/i386/fpu/s_fma.S.jj	2009-05-16 19:23:38.000000000 +0200
+++ libc/sysdeps/i386/fpu/s_fma.S	2010-09-14 11:29:01.831482015 +0200
@@ -1,31 +0,0 @@
-/* Compute (X * Y) + Z as ternary operation.
-   Copyright (C) 1997, 1998 Free Software Foundation, Inc.
-   This file is part of the GNU C Library.
-   Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997.
-
-   The GNU C Library is free software; you can redistribute it and/or
-   modify it under the terms of the GNU Lesser General Public
-   License as published by the Free Software Foundation; either
-   version 2.1 of the License, or (at your option) any later version.
-
-   The GNU C Library is distributed in the hope that it will be useful,
-   but WITHOUT ANY WARRANTY; without even the implied warranty of
-   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
-   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.  */
-
-#include <sysdep.h>
-
-	.text
-ENTRY(__fma)
-	fldl	4(%esp)		// x
-	fmull	12(%esp)	// x * y
-	fldl	20(%esp)	// z : x * y
-	faddp			// (x * y) + z
-	ret
-END(__fma)
-weak_alias (__fma, fma)
--- libc/sysdeps/i386/fpu/s_fmal.S.jj	2009-05-16 19:23:38.000000000 +0200
+++ libc/sysdeps/i386/fpu/s_fmal.S	2010-09-14 11:29:01.831482015 +0200
@@ -1,32 +0,0 @@
-/* Compute (X * Y) + Z as ternary operation.
-   Copyright (C) 1997 Free Software Foundation, Inc.
-   This file is part of the GNU C Library.
-   Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997.
-
-   The GNU C Library is free software; you can redistribute it and/or
-   modify it under the terms of the GNU Lesser General Public
-   License as published by the Free Software Foundation; either
-   version 2.1 of the License, or (at your option) any later version.
-
-   The GNU C Library is distributed in the hope that it will be useful,
-   but WITHOUT ANY WARRANTY; without even the implied warranty of
-   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
-   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.  */
-
-#include <sysdep.h>
-
-	.text
-ENTRY(__fmal)
-	fldt	4(%esp)		// x
-	fldt	16(%esp)	// x : y
-	fmulp			// x * y
-	fldt	28(%esp)	// z : x * y
-	faddp			// (x * y) + z
-	ret
-END(__fmal)
-weak_alias (__fmal, fmal)


	Jakub

[-- Attachment #2: testfmas.c --]
[-- Type: text/plain, Size: 4556 bytes --]

#define _GNU_SOURCE
#include <stdlib.h>
#include <mpfr.h>
#include <ieee754.h>
#include <stdio.h>
#include <math.h>

#ifdef DO_FLT
#define MANT_DIG __FLT_MANT_DIG__
#define MAX_EXP __FLT_MAX_EXP__
#define union_type union ieee754_float
#define suff(x) x##d
#define fma_suff fmaf
#define BIAS 0x7f
#define FMT "%a"
#define fltfield f
#elif defined DO_DBL
#define MANT_DIG __DBL_MANT_DIG__
#define MAX_EXP __DBL_MAX_EXP__
#define union_type union ieee754_double
#define suff(x) x##d
#define fma_suff fma
#define BIAS 0x3ff
#define FMT "%a"
#define fltfield d
#elif defined DO_LDBL
#define MANT_DIG __LDBL_MANT_DIG__
#define MAX_EXP __LDBL_MAX_EXP__
#define union_type union ieee854_long_double
#define suff(x) x##ld
#define fma_suff fmal
#define BIAS 0x3fff
#if MANT_DIG == 64
#define FMT "%La"
#else
#define FMT "%.40La"
#endif
#define fltfield d
#else
#error Define DO_FLT, DO_DBL or DO_LDBL
#endif

unsigned int
randbits (int count)
{
  static unsigned int val, bits;
  unsigned int ret = 0;

  if (count > bits)
    {
      if (count == 32 && bits == 0)
	{
	  ret = ((unsigned int) random ()) << 1;
	  count -= 31;
	}
      else
	{
	  ret = val << (count - bits);
	  count -= bits;
	}
      val = (random () & 0x7fffffffU);
      bits = 31;
    }
  ret |= val & ((1U << count) - 1);
  val >>= count;
  bits -= count;
  return ret;
}

void
randflt (int mode, union_type *p, unsigned int exp)
{
  if (mode == 1 || (mode >= 2 && randbits (4) == 0))
    {
      unsigned int mantissa[4] = { 0, 0, 0, 0 };
      unsigned int bits = 0, cur = randbits (1);
      do
	{
	  unsigned int curbits = randbits (8) % (MANT_DIG - bits + 6);
	  if (curbits > MANT_DIG - bits)
	    curbits = MANT_DIG - bits;
	  if (cur && curbits)
	    {
	      do
		{
		  unsigned int thisbits = 32 - (bits & 31);
		  if (curbits < thisbits)
		    thisbits = curbits;
		  mantissa[bits / 32] |= (~0U << (32 - thisbits))
					 >> (32 - thisbits - (bits & 31));
		  bits += thisbits;
		  curbits -= thisbits;
		}
	      while (curbits);
	    }
	  else
	    bits += curbits;
	  cur ^= 1;
	}
      while (bits < MANT_DIG);
#ifdef DO_FLT
      p->ieee.mantissa = mantissa[0];
#elif defined DO_DBL
      p->ieee.mantissa1 = mantissa[0];
      p->ieee.mantissa0 = mantissa[1];
#elif defined DO_LDBL
#if MANT_DIG == 64
      p->ieee.mantissa1 = mantissa[0];
      p->ieee.mantissa0 = mantissa[1] & 0x7fffffffU;
#elif MANT_DIG == 113
      p->ieee.mantissa3 = mantissa[0];
      p->ieee.mantissa2 = mantissa[1];
      p->ieee.mantissa1 = mantissa[2];
      p->ieee.mantissa0 = mantissa[3];
#endif
#endif
    }
  else
    {
#ifdef DO_FLT
      p->ieee.mantissa = randbits (23);
#elif defined DO_DBL
      p->ieee.mantissa1 = randbits (32);
      p->ieee.mantissa0 = randbits (20);
#elif defined DO_LDBL
#if MANT_DIG == 64
      p->ieee.mantissa1 = randbits (32);
      p->ieee.mantissa0 = randbits (31);
#elif MANT_DIG == 113
      p->ieee.mantissa3 = randbits (32);
      p->ieee.mantissa2 = randbits (32);
      p->ieee.mantissa1 = randbits (32);
      p->ieee.mantissa0 = randbits (16);
#endif
#endif
    }
  p->ieee.negative = randbits (1);
  p->ieee.exponent = exp;
#if defined (DO_LDBL) && MANT_DIG == 64
  p->ieee.empty = 0;
  if (exp != 0)
    p->ieee.mantissa0 |= 0x80000000;
#endif
}

int
main (int argc, const char **argv)
{
  long int i, count = 100;

  if (argc >= 2)
    count = strtoul (argv[1], NULL, 0);
  mpfr_set_default_prec (MANT_DIG);
  mpfr_set_emin (-MAX_EXP-MANT_DIG+4);
  mpfr_set_emax (MAX_EXP);
  for (i = 0; i < count; i++)
    {
      union_type u, v, w, x, y;
      int mode = randbits (2);
      randflt (mode, &u, randbits (15));
      randflt (mode, &v, randbits (15));
      if (randbits (4))
	randflt (mode, &w, ((int) u.ieee.exponent + (int) v.ieee.exponent)
			   - BIAS - MANT_DIG + randbits (18) % (2 * MANT_DIG));
      else
	randflt (mode, &w, randbits (15));
      mpfr_t xx, xy, xz;
      suff (mpfr_init_set_) (xx, u.fltfield, GMP_RNDN);
      suff (mpfr_init_set_) (xy, v.fltfield, GMP_RNDN);
      suff (mpfr_init_set_) (xz, w.fltfield, GMP_RNDN);
      int i = mpfr_fma (xx, xx, xy, xz, GMP_RNDN);
      mpfr_subnormalize (xx, i, GMP_RNDN);
      x.fltfield = suff (mpfr_get_) (xx, GMP_RNDN);
      mpfr_clear (xx); mpfr_clear (xy); mpfr_clear (xz);
      y.fltfield = fma_suff (u.fltfield, v.fltfield, w.fltfield);
      if (x.fltfield != y.fltfield
	  && !(isnan (x.fltfield) && isnan (y.fltfield)))
	printf ("fma ("FMT", "FMT", "FMT") = "FMT" mpfr_fma "FMT"\n",
		u.fltfield, v.fltfield, w.fltfield, y.fltfield, x.fltfield);
    }
  return 0;
}

^ permalink raw reply	[flat|nested] only message in thread

only message in thread, other threads:[~2010-10-15 19:08 UTC | newest]

Thread overview: (only message) (download: mbox.gz / follow: Atom feed)
-- links below jump to the message on this page --
2010-10-15 19:08 [PATCH] Implement fmal, some fma bugfixes Jakub Jelinek

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).