[Bug math/706] pow() produces inaccurate results for base ~ 1.0, and large exponent on 32-bit x86

["vincent-srcware at vinc17 dot net" <sourceware-bugzilla at sourceware dot org>]

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

--- Comment #13 from Vincent LefÃvre <vincent-srcware at vinc17 dot net> 2012-02-28 14:57:42 UTC ---
(In reply to comment #10)
> Note that if you want to compute pow(x,y) where y is not necessarily an
> integer, you can write y = yi + yf, where yi is an integer and |yf| < 1
> (e.g. |yf| <= 0.5 or 0 <= yf < 1). Then mathematically,
> pow(x,y) = pow(pow(x,yi),yf), and you can use that to get a more accurate
> result than the current algorithm (I haven't checked in detail...).

This is incorrect. I should have re-read my message. There are 2 possible
simple range reductions that wouldn't require much rewriting, in order to
reduce the problem to integer exponents (if one is able to fix that easily).

1. y = yi + yf, where yi is an integer and |yf| < 1. Then mathematically,
pow(x,y) = pow(x,yi) * pow(x,yf), where pow(x,yi) can be computed with an
iterated exponentiation. Note that yi and yf are computed exactly.

2. y = y0 * 2^e, where 0.5 <= |y0| < 1 for instance (which can be obtained with
frexp()). Then mathematically, pow(x,y) = pow(pow(x,y0),2^e), where the ^(2^e)
can be computed with repeated squarings.

Now, I fear that if the iterated exponentiation is not done in multiple
precision, this won't solve the problem for very large values of n. Indeed,
pow(x,n) evaluated with such an algorithm (whatever iteration is used) is
approximated with a factor of about[*] (1+epsilon)^(n-1) > 1 + (n-1).epsilon,
where 1+epsilon is the error bound on the multiplication in Higham's notation.
And this would be slow.

[*] Actually this is an upper bound, but in average, you would probably get the
same kind of problems with large values of n.

And Calixte, I've looked at your comment #9 and your code more closely, and I
agree that in the case of a negative exponent, the inversion should be done
after the exponentiation if such an algorithm is used.

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