This is the mail archive of the libc-alpha@sourceware.org mailing list for the glibc project.


Index Nav: [Date Index] [Subject Index] [Author Index] [Thread Index]
Message Nav: [Date Prev] [Date Next] [Thread Prev] [Thread Next]
Other format: [Raw text]

Re: XFAIL libm-test.inc tests as needed for ibm128 [committed]


On Wed, 11 Jan 2017, Tulio Magno Quites Machado Filho wrote:

> > (The case of patched libgcc was already clean up to ulps
> 
> Which patch are you referring to?

This one (only suitable for building GCC specifically for testing glibc, 
not for more general use, since the changes are unconditional and make 
performance significantly worse).  This patch doesn't fully support soft 
float (that would require an additional interface exported from glibc) and 
as a minimal patch doesn't include testcases for the various bugs fixed 
(although I have such testcases in an earlier patch version).  With this 
patch it should be possible to define TEST_COND_ibm128_libgcc to 0 in 
libm-test.inc, and also possible to do random test generation to find 
cases of large errors in the ldbl-128ibm functions without running into 
problems with libgcc issues (being able to use random test generation to 
find bugs in glibc is probably the clearest use of this patch).

Index: libgcc/config/rs6000/ibm-ldouble.c
===================================================================
--- libgcc/config/rs6000/ibm-ldouble.c	(revision 228280)
+++ libgcc/config/rs6000/ibm-ldouble.c	(working copy)
@@ -47,6 +47,10 @@ see the files COPYING3 and COPYING.RUNTIME respect
 
 #if defined (__MACH__) || defined (__powerpc__) || defined (_AIX)
 
+typedef _Bool bool;
+#define false 0
+#define true 1
+
 #define fabs(x) __builtin_fabs(x)
 #define isless(x, y) __builtin_isless (x, y)
 #define inf() __builtin_inf()
@@ -87,6 +91,58 @@ __asm__ (".symver __gcc_qadd,_xlqadd@GCC_3.4\n\t"
 	 ".symver .__gcc_qdiv,._xlqdiv@GCC_3.4");
 #endif
 
+#define FE_TONEAREST 0
+
+/* Set the rounding mode to MODE, returning the previous mode.  */
+
+static int
+swapround (int mode)
+{
+  int old_mode;
+#ifdef __NO_FPRS__
+# ifdef _SOFT_FLOAT
+  /* Not implemented.  */
+  old_mode = FE_TONEAREST;
+# else
+  /* SPE floating point.  */
+  int spefscr;
+  asm volatile ("mfspefscr %0" : "=r" (spefscr));
+  old_mode = spefscr & 3;
+  spefscr = (spefscr & ~3) | mode;
+  asm volatile ("mtspefscr %0" : : "r" (spefscr));
+# endif
+#else
+  /* Logic as in glibc.  */
+  asm volatile ("mcrfs 7,7\n\t"
+		"mfcr %0" : "=r" (old_mode) : : "cr7");
+  old_mode &= 3;
+  if ((mode ^ old_mode) & 2)
+    {
+      if (mode & 2)
+	asm volatile ("mtfsb1 30");
+      else
+	asm volatile ("mtfsb0 30");
+    }
+  if ((mode ^ old_mode) & 1)
+    {
+      if (mode & 1)
+	asm volatile ("mtfsb1 31");
+      else
+	asm volatile ("mtfsb0 31");
+    }
+#endif
+  return old_mode;
+}
+
+/* Overflowing values in each rounding mode (FE_TONEAREST,
+   FE_TOWARDZERO, FE_UPWARD, FE_DOWNWARD).  */
+
+static const long double overflow_values[2][4] =
+  {
+    { __builtin_infl (), __LDBL_MAX__, __builtin_infl (), __LDBL_MAX__ },
+    { -__builtin_infl (), -__LDBL_MAX__, -__LDBL_MAX__, -__builtin_infl () }
+  };
+
 /* Combine two 'double' values into one 'long double' and return the result.  */
 static inline long double
 pack_ldouble (double dh, double dl)
@@ -110,40 +166,74 @@ pack_ldouble (double dh, double dl)
 long double
 __gcc_qadd (double a, double aa, double c, double cc)
 {
-  double xh, xl, z, q, zz;
+  double xh, xhs, xl, z, q, zz;
+  bool scale = false;
+  int old_mode;
 
+  if (nonfinite (a) || nonfinite (c))
+    return a + c;
+
+  old_mode = swapround (FE_TONEAREST);
+
+  if (fabs (a) >= 0x1p1022 || fabs (c) >= 0x1p1022)
+    {
+      scale = true;
+      if (fabs (a) >= 0x1p914)
+	{
+	  a *= 0.5;
+	  aa *= 0.5;
+	}
+      if (fabs (c) >= 0x1p914)
+	{
+	  c *= 0.5;
+	  cc *= 0.5;
+	}
+    }
+
   z = a + c;
+  q = a - z;
+  zz = q + c + (a - (q + z)) + aa + cc;
 
-  if (nonfinite (z))
+  /* Keep -0 result.  */
+  if (zz == 0.0)
     {
-      if (fabs (z) != inf())
-	return z;
-      z = cc + aa + c + a;
-      if (nonfinite (z))
-	return z;
-      xh = z;  /* Will always be DBL_MAX.  */
-      zz = aa + cc;
-      if (fabs(a) > fabs(c))
-	xl = a - z + c + zz;
-      else
-	xl = c - z + a + zz;
+      double ret = scale ? 2.0 * z : z;
+      if (unlikely (old_mode != FE_TONEAREST))
+	{
+	  swapround (old_mode);
+	  if (fabs (ret) == inf ())
+	    return overflow_values[ret < 0][old_mode];
+	  if (ret == 0)
+	    {
+	      asm volatile ("" : "+m" (a));
+	      ret = a + c;
+	    }
+	}
+      return ret;
     }
-  else
+  xh = z + zz;
+  xhs = xh;
+  if (scale)
+    xhs *= 2.0;
+  if (nonfinite (xhs))
     {
-      q = a - z;
-      zz = q + c + (a - (q + z)) + aa + cc;
+      if (unlikely (old_mode != FE_TONEAREST))
+	{
+	  swapround (old_mode);
+	  if (fabs (xhs) == inf ())
+	    return overflow_values[xhs < 0][old_mode];
+	}
+      return xhs;
+    }
 
-      /* Keep -0 result.  */
-      if (zz == 0.0)
-	return z;
+  xl = z - xh + zz;
+  if (scale)
+    xl *= 2.0;
 
-      xh = z + zz;
-      if (nonfinite (xh))
-	return xh;
+  if (unlikely (old_mode != FE_TONEAREST))
+    swapround (old_mode);
 
-      xl = z - xh + zz;
-    }
-  return pack_ldouble (xh, xl);
+  return pack_ldouble (xhs, xl);
 }
 
 long double
@@ -159,13 +249,41 @@ static double fmsub (double, double, double);
 long double
 __gcc_qmul (double a, double b, double c, double d)
 {
-  double xh, xl, t, tau, u, v, w;
+  double xh, xl, t, tau, u, us, v, w;
+  bool scale = false;
+  int old_mode;
+
+  if (nonfinite (a) || nonfinite (c))
+    return a * c;
+
+  old_mode = swapround (FE_TONEAREST);
+
+  if (fabs (a) >= 0x1p512)
+    {
+      scale = true;
+      a *= 0.5;
+      b *= 0.5;
+    }
+  else if (fabs (c) >= 0x1p512)
+    {
+      scale = true;
+      c *= 0.5;
+      d *= 0.5;
+    }
   
   t = a * c;			/* Highest order double term.  */
 
   if (unlikely (t == 0)		/* Preserve -0.  */
       || nonfinite (t))
-    return t;
+    {
+      if (unlikely (old_mode != FE_TONEAREST))
+	{
+	  swapround (old_mode);
+	  if (fabs (t) == inf ())
+	    return overflow_values[t < 0][old_mode];
+	}
+      return t;
+    }
 
   /* Sum terms of two highest orders. */
   
@@ -179,12 +297,29 @@ __gcc_qmul (double a, double b, double c, double d
   w = b*c;
   tau += v + w;	    /* Add in other second-order terms.	 */
   u = t + tau;
+  us = u;
+  if (scale)
+    us *= 2.0;
 
   /* Construct long double result.  */
-  if (nonfinite (u))
-    return u;
-  xh = u;
+  if (nonfinite (us))
+    {
+      if (unlikely (old_mode != FE_TONEAREST))
+	{
+	  swapround (old_mode);
+	  if (fabs (us) == inf ())
+	    return overflow_values[us < 0][old_mode];
+	}
+      return us;
+    }
+  xh = us;
   xl = (t - u) + tau;
+  if (scale)
+    xl *= 2.0;
+
+  if (unlikely (old_mode != FE_TONEAREST))
+    swapround (old_mode);
+
   return pack_ldouble (xh, xl);
 }
 
@@ -191,13 +326,41 @@ __gcc_qmul (double a, double b, double c, double d
 long double
 __gcc_qdiv (double a, double b, double c, double d)
 {
-  double xh, xl, s, sigma, t, tau, u, v, w;
+  double xh, xl, s, sigma, t, tau, u, us, v, w;
+  bool scale = false;
+  int old_mode;
+
+  if (nonfinite (a) || nonfinite (c) || c == 0)
+    return a / c;
+
+  old_mode = swapround (FE_TONEAREST);
+
+  if (fabs (a) >= 0x1p512)
+    {
+      scale = true;
+      a *= 0.5;
+      b *= 0.5;
+    }
+  else if (fabs (c) <= 0x1p-511)
+    {
+      scale = true;
+      c *= 2.0;
+      d *= 2.0;
+    }
   
   t = a / c;                    /* highest order double term */
   
   if (unlikely (t == 0)		/* Preserve -0.  */
       || nonfinite (t))
-    return t;
+    {
+      if (unlikely (old_mode != FE_TONEAREST))
+	{
+	  swapround (old_mode);
+	  if (fabs (t) == inf ())
+	    return overflow_values[t < 0][old_mode];
+	}
+      return t;
+    }
 
   /* Finite nonzero result requires corrections to the highest order
      term.  These corrections require the low part of c * t to be
@@ -224,12 +387,29 @@ __gcc_qdiv (double a, double b, double c, double d
   
   tau = ((v-sigma)+w)/c;   /* Correction to t.  */
   u = t + tau;
+  us = u;
+  if (scale)
+    us *= 2.0;
 
   /* Construct long double result.  */
-  if (nonfinite (u))
-    return u;
-  xh = u;
+  if (nonfinite (us))
+    {
+      if (unlikely (old_mode != FE_TONEAREST))
+	{
+	  swapround (old_mode);
+	  if (fabs (us) == inf ())
+	    return overflow_values[us < 0][old_mode];
+	}
+      return us;
+    }
+  xh = us;
   xl = (t - u) + tau;
+  if (scale)
+    xl *= 2.0;
+
+  if (unlikely (old_mode != FE_TONEAREST))
+    swapround (old_mode);
+
   return pack_ldouble (xh, xl);
 }
 

-- 
Joseph S. Myers
joseph@codesourcery.com


Index Nav: [Date Index] [Subject Index] [Author Index] [Thread Index]
Message Nav: [Date Prev] [Date Next] [Thread Prev] [Thread Next]