Fix IBM long double division inaccuracy (glibc bug 15396)

[Joseph S. Myers]

This patch fixes inaccuracy of IBM long double division that causes
large errors in glibc testing
<https://sourceware.org/bugzilla/show_bug.cgi?id=15396>.  The
implementation computes a first approximation to the result as t = a /
c, dividing the high parts of the arguments, then does adjustments
involving calculating c * t as the exact sum of two double values.  In
the problem cases, the low part of c * t underflows, so resulting in
loss of precision compared to the precision available for a number
with the exponent of the final result.  This patch scales the
arguments up in the problem case.

Tested with no regressions for cross to powerpc-linux-gnu.  OK to
commit?

(Note that there remain other bugs in the IBM long double code, some
causing glibc test failures, at least (a) invalid results in rounding
modes other than FE_TONEAREST, (b) spurious overflow and underflow
exceptions, mainly but not entirely where discontiguous mantissa bits
are involved.)

libgcc:
2014-01-02  Joseph Myers  <joseph@codesourcery.com>

	* config/rs6000/ibm-ldouble.c (__gcc_qdiv): Scale up arguments in
	case of small numerator and finite nonzero result.

gcc/testsuite:
2014-01-02  Joseph Myers  <joseph@codesourcery.com>

	* gcc.target/powerpc/rs6000-ldouble-3.c: New test.

Index: gcc/testsuite/gcc.target/powerpc/rs6000-ldouble-3.c
===================================================================
--- gcc/testsuite/gcc.target/powerpc/rs6000-ldouble-3.c	(revision 0)
+++ gcc/testsuite/gcc.target/powerpc/rs6000-ldouble-3.c	(revision 0)
@@ -0,0 +1,21 @@
+/* Test accuracy of long double division (glibc bug 15396).  */
+/* { dg-do run { target powerpc*-*-linux* powerpc*-*-darwin* powerpc*-*-aix* rs6000-*-* } } */
+/* { dg-options "-mlong-double-128" } */
+
+extern void exit (int);
+extern void abort (void);
+
+volatile long double a = 0x1p-1024L;
+volatile long double b = 0x3p-53L;
+volatile long double r;
+volatile long double expected = 0x1.55555555555555555555555555p-973L;
+
+int
+main (void)
+{
+  r = a / b;
+  /* Allow error up to 2ulp.  */
+  if (__builtin_fabsl (r - expected) > 0x1p-1073L)
+    abort ();
+  exit (0);
+}
Index: libgcc/config/rs6000/ibm-ldouble.c
===================================================================
--- libgcc/config/rs6000/ibm-ldouble.c	(revision 206280)
+++ libgcc/config/rs6000/ibm-ldouble.c	(working copy)
@@ -190,7 +190,16 @@ __gcc_qdiv (double a, double b, double c, double d
       || nonfinite (t))
     return t;
 
-  /* Finite nonzero result requires corrections to the highest order term.  */
+  /* Finite nonzero result requires corrections to the highest order
+     term.  These corrections require the low part of c * t to be
+     exactly represented in double.  */
+  if (fabs (a) <= 0x1p-969)
+    {
+      a *= 0x1p106;
+      b *= 0x1p106;
+      c *= 0x1p106;
+      d *= 0x1p106;
+    }
 
   s = c * t;                    /* (s,sigma) = c*t exactly.  */
   w = -(-b + d * t);	/* Written to get fnmsub for speed, but not

-- 
Joseph S. Myers
joseph@codesourcery.com