[Bug math/13957] New: powl very inaccurate on powerpc

["bruno at clisp dot org" <sourceware-bugzilla at sourceware dot org>]

http://sourceware.org/bugzilla/show_bug.cgi?id=13957

             Bug #: 13957
           Summary: powl very inaccurate on powerpc
           Product: glibc
           Version: 2.11
            Status: NEW
          Severity: normal
          Priority: P2
         Component: math
        AssignedTo: unassigned@sourceware.org
        ReportedBy: bruno@clisp.org
    Classification: Unclassified


Created attachment 6327
  --> http://sourceware.org/bugzilla/attachment.cgi?id=6327
test case

While on x86, x86_64 platforms the results of powl() generally have an error
of less than 100 ulps, on a PowerPC platform I'm seeing errors of more than
50000 ulps.

Recall that on PowerPC, 'long double' numbers have 106 mantissa bits,
i.e. ca. 32 decimal digits after the decimal point.

How to reproduce:
================================ foo.c ================================
#include <math.h>
#include <stdio.h>
long double x = 0.9500175466493529918258722439486448724157L;
long double y = 13499.12711242096089256937988377554090961L;
int main ()
{
  long double p = powl (x, y);
  long double q = powl (p, 1.0L / y);
  printf ("x = %.63Lg\n", x);
  printf ("y = %.63Lg\n", y);
  printf ("x^y = %.63Lg\n", p);
  printf ("x^y*2^1000 = %.63Lg\n", ldexpl (p, 1000));
  printf ("(x^y)^(1/y) = %.63Lg\n", q);
  printf ("(x^y)^(1/y)/x = %.63Lg\n", q / x);
  return 0;
}
=======================================================================
$ gcc -Wall foo.c -lm
$ ./a.out

Expected results (assuming no rounding errors at all):

x = 0.950017546649352991825872243948644872415686653396830865564350738
y = 13499.1271124209608925693798837755409096128741713277563070427778
x^y = 2.49114074548890053222918160895394492812126078716433839487765876e-301
x^y*2^1000 = 2.66927875050377145601812886620084351250727601302586070667964568
(x^y)^(1/y) = 0.950017546649352991825872243948644872415686653396830865564350738
(x^y)^(1/y)/x = 1

Actual results:

x = 0.950017546649352991825872243948644872415686653396830865564350738
y = 13499.1271124209608925693798837755409096128741713277563070427778
x^y = 2.49114074548890047185112002893359654155432307599856022554444831e-301
x^y*2^1000 = 2.66927875050377145601815149518866790434579172597295837476849556
(x^y)^(1/y) = 0.950017546649352991825872244545270235795653144477288826094787998
(x^y)^(1/y)/x = 1.0000000000000000000000000006280072362657898669644932875180931

You can see:
1) x^y*2^1000 has only 22 correct digits after the decimal point. This means
that x^y has an error of more than 900000000 ulp!
2) (x^y)^(1/y)/x has only 27 correct digits after the decimal point. This means
that (x^y)^(1/y)/x has an error of more than 50000 ulp!

-- 
Configure bugmail: http://sourceware.org/bugzilla/userprefs.cgi?tab=email
------- You are receiving this mail because: -------
You are on the CC list for the bug.