Bad error computing 1/z
Stephen Montgomery-Smith
stephen at missouri.edu
Tue Jul 31 20:46:23 UTC 2012
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?
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.
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.
More information about the freebsd-numerics
mailing list