Bug 706 reports that the x86 pow function is inaccurate for large
integer exponents because it uses an iterative multiplication
algorithm that allows rounding errors to accumulate.
I propose this patch to fix the bug by using the log+exp real-exponent
algorithm for integer exponents with absolute value at least 1024.
The choice of that exponent is somewhat arbitrary. Both algorithms
take advantage of having 11 bits excess precision, which is sufficient
for all inputs for the log+exp algorithm (you need about 53 fractional
bits in y*log2(x), and get them unless there is overflow or underflow
because without overflow or underflow there will be at most 10 bits in
the integer part of y*log2(x)) but not for exponents much bigger than
1024 with the multiplication algorithm. In terms of accuracy there is
probably no real advantage in the multiplication algorithm once the
exponent is large enough that no cases other than powers of 2 are
exact (i.e. once power of 3 are no longer exactly representable), but
I don't know how the speed of the algorithms compares; choosing a
threshold as large as possible is conservative in that it avoids
changing the algorithm unnecessarily.
Given that change, the real-exponent code then needs to handle
negative x with odd integer exponent, which is done using code similar
to that used in powl which already needed to handle that issue.
This bug does not affect powf since 40 bits excess precision is enough
for any integer exponent that uses the multiplication algorithm in x86
powf. It does affect powl, but for many inputs the real-exponent
algorithm is even worse there (see bug 13881) so this patch does not
change powl (but a fix for bug 13881 will need to stop the
integer-exponent algorithm being used for all but very small
exponents, as well as improving the real-exponent code).
Tested x86 and x86_64 and ulps updated accordingly (x86 didn't need
any new ulps, x86_64 did need ulps for some new tests for powf).
2012-04-08 Joseph Myers<joseph@codesourcery.com>
[BZ #706]
* sysdeps/i386/fpu/e_pow.S (p10): New object.
(__ieee754_pow): Use iterative multiplication algorithm only for
integer exponents with absolute value below 1024. Check for odd
integer exponents when using algorithm for real exponents.
* math/libm-test.inc (pow_test): Add more tests.
* sysdeps/x86_64/fpu/libm-test-ulps: Update.