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

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

[bruno at clisp 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.

[Bug math/13957] powl very inaccurate on powerpc

[bruno at clisp dot org]

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

Bruno Haible <bruno at clisp dot org> changed:

           What    |Removed                     |Added
----------------------------------------------------------------------------
           See Also|                            |http://sourceware.org/bugzi
                   |                            |lla/show_bug.cgi?id=13881

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

[Bug math/13957] powl very inaccurate on powerpc

[jsm28 at gcc dot gnu.org]

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

Joseph Myers <jsm28 at gcc dot gnu.org> changed:

           What    |Removed                     |Added
----------------------------------------------------------------------------
             Status|NEW                         |RESOLVED
         Resolution|                            |INVALID

--- Comment #1 from Joseph Myers <jsm28 at gcc dot gnu.org> 2012-11-06 16:38:42 UTC ---
Actually, I don't think this case demonstrates a significant error, although
maybe other cases do.

You give an expected result of
2.49114074548890053222918160895394492812126078716433839487765876e-301.  That
has binary exponent -999.  The least normal binary exponent for IBM long double
is -969.  So there are only 76 mantissa bits available for this number - that
is, 22 correct decimal digits after the point is the most you could expect.

Then you raise a number, with an error of about one part in 2^75 from the
mathematical value of x^y, to the power 1/13499, resulting in a number expected
to have an error of about one part in 2^89.  That corresponds pretty well to
the 27 digits after the point.

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

[Bug math/13957] powl very inaccurate on powerpc

[jsm28 at gcc dot gnu.org]

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

--- Comment #2 from Joseph Myers <jsm28 at gcc dot gnu.org> 2012-11-06 16:43:31 UTC ---
If you do find cases that show significant errors, I advise just giving the x
and y values as hex floats, and the observed and expected results as hex
floats, to avoid any issues of exactness of decimal representations, and to
make clear exactly what call is inaccurate rather than making the program do
internal consistency checks such as (x^y)^(1/y); we can easily confirm the
correct expected values using MPFR.

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

[Bug math/13957] powl very inaccurate on powerpc

[fweimer at redhat dot com]

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

Florian Weimer <fweimer at redhat dot com> changed:

           What    |Removed                     |Added
----------------------------------------------------------------------------
              Flags|                            |security-

-- 
You are receiving this mail because:
You are on the CC list for the bug.