Bad error computing 1/z

[Stephen Montgomery-Smith]

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.