Bug 13881 - powl inaccurate for x86/x86_64
Summary: powl inaccurate for x86/x86_64
Status: RESOLVED FIXED
Alias: None
Product: glibc
Classification: Unclassified
Component: math (show other bugs)
Version: 2.15
: P2 normal
Target Milestone: ---
Assignee: Not yet assigned to anyone
URL:
Keywords:
Depends on:
Blocks:
 
Reported: 2012-03-20 18:02 UTC by Joseph Myers
Modified: 2019-05-26 12:38 UTC (History)
4 users (show)

See Also:
Host:
Target:
Build:
Last reconfirmed:
fweimer: security-


Attachments
test case (311 bytes, text/x-csrc)
2012-04-06 15:11 UTC, Bruno Haible
Details
Test case for powl accuracy (476 bytes, text/plain)
2019-05-25 06:04 UTC, agner@agner.org
Details

Note You need to log in before you can comment on or make changes to this bug.
Description Joseph Myers 2012-03-20 18:02:17 UTC
The powl implementations for x86 and x86_64 work essentially through exp2 (y * log2 (x)).  This does not produce accurate results for large outputs; you need about 15 extra bits of precision in the evaluation of y * log2 (x) to do that.  As an example, powl (1e4932L, 0.75L) is 4840ulp off from the expected 1e3699L on x86_64.

For integer exponents a different algorithm is used; it's still not very accurate, e.g. 45ulp error for powl (10, 4932); that has the same underlying cause as the inaccuracy described for pow in bug 706.
Comment 1 Bruno Haible 2012-04-06 15:11:30 UTC
Created attachment 6328 [details]
test case
Comment 2 Bruno Haible 2012-04-06 15:24:30 UTC
Find attached a test case. The approximations to 1e4932L and 1e3699L
are computed through an (exact) ldexpl call, from a mantissa that was
computed with GNU clisp's high-precision long-float numbers.

$ gcc -Wall foo.c -lm
$ ./a.out

Expected result (without any rounding errors):

x = 1e+4932
y = 0.75
x^y = 1e+3699
1

Actual result (on x86 and x86_64):

x = 1.000000000000000000006018949387963891324e+4932
y = 0.75
x^y = 1.000000000000000298896762326719536039811e+3699
1.000000000000000298914538954253572455855

As you can see:
1) x^y has only 15 valid digits after the decimal point.
2) The printf of x has 22 correct digits after the decimal point; this matches
the 64 mantissa bits in an x86 'long double'.
Comment 3 Joseph Myers 2012-11-28 13:56:53 UTC
Fixed for 2.17 by:

commit 1bead169c32a3a688de863709b863207b7aafddd
Author: Joseph Myers <joseph@codesourcery.com>
Date:   Wed Nov 28 13:40:54 2012 +0000

    Fix powl inaccuracy for x86_64 and x86 (bug 13881).
Comment 4 agner@agner.org 2019-05-25 06:04:48 UTC
Created attachment 11802 [details]
Test case for powl accuracy
Comment 5 agner@agner.org 2019-05-25 06:07:36 UTC
Please reopen.

powl is still less accurate than pow in version 2.27 (both gcc and clang).

The attached test case powl.cpp has powl 5 times more inaccurate than pow. I have seen worse cases, but I have chosen a test case where it is easy to calculate the correct value. I have not seen any case where powl is more accurate than pow.
Comment 6 agner@agner.org 2019-05-25 07:28:09 UTC
Some more test cases:

powl(9.4401197044207130E-01,9.0396338480980503E+03) relative error 7.9E-13

powl(9.6205229632635447E-01,9.0396338480980503E+03) relative error 5.1E-13

powl(1.0567906354279399E+00,9.0396338480980503E+03) relative error 5.7E-13

powl(9.5952957579581444E-01,-1.6994966845030955E+04) relative error 1.8E-12

powl(9.4249308583367331E-01,5.1811882606772378E+03) relative error 6.4E-13

powl(1.0537202715568794E+00,5.1811882606772378E+03) relative error 4.2E-13
Comment 7 Florian Weimer 2019-05-26 12:38:09 UTC
(In reply to agner@agner.org from comment #5)
> Please reopen.
> 
> powl is still less accurate than pow in version 2.27 (both gcc and clang).
> 
> The attached test case powl.cpp has powl 5 times more inaccurate than pow. I
> have seen worse cases, but I have chosen a test case where it is easy to
> calculate the correct value. I have not seen any case where powl is more
> accurate than pow.

Would you please file a separate bug for that?  Thanks.