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 <jos...@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 <jos...@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 jos...@codesourcery.com