Created attachment 6327 [details] 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!
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.
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.