Fix x86 pow inaccuracy for large integer exponents (bug 706)

[Andreas Jaeger]

On 04/08/2012 07:59 PM, Joseph S. Myers wrote:
> 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.

Thanks, this is fine.

Andreas
-- 
  Andreas Jaeger aj@{suse.com,opensuse.org} Twitter/Identica: jaegerandi
   SUSE LINUX Products GmbH, Maxfeldstr. 5, 90409 Nürnberg, Germany
    GF: Jeff Hawn,Jennifer Guild,Felix Imendörffer,HRB16746 (AG Nürnberg)
     GPG fingerprint = 93A3 365E CE47 B889 DF7F  FED1 389A 563C C272 A126