Bad error computing 1/z
Bruce Evans
brde at optusnet.com.au
Wed Aug 1 08:13:04 UTC 2012
On Tue, 31 Jul 2012, Stephen Montgomery-Smith wrote:
> On 07/29/12 22:35, Bruce Evans wrote:
>
>> 1/z calculated by gcc is much more accurate than the singular function
>> z*z :-). The real part of the latter is x*x-y*y. This may have an
>> error of almost 2**MANT_DIG ulps due to cancelation.
>
> Could one avoid the cancellation errors in computing x*x-y*y by using the
> formula (x-y)*(x+y) instead?
Yes, but there seems to be no similar trick for general complex
multiplication or higher powers of z.
> Here is my proposed argument: Case 1: x and y have the same sign. If they
> are nearly equal, then the computation (x-y) should be exact. The error in
> computing (x+y) shouldn't be worse than 0.5 ULP. So the whole calculation
> could probably be performed with an error less than 1 ULP, maybe even 0.5 ULP
> if one took a lot of care. Case 2: x and y have different signs - just
> switch the roles of (x+y) and (x-y) in the preceding argument.
1.5 or even 2.0 ulps I think (since there are 3 operations).
> Do you know for a fact that gcc uses x*x-y*y? To use (x-y)*(x+y) would be no
> worse performance wise - perhaps better if addition is cheaper than
> multiplication. It would be the simplest programming switch in the world, as
> long as one knew where to find the code fragment that needed to be changed.
I checked, and it doesn't seem to.
Bruce
More information about the freebsd-numerics
mailing list