Skip site navigation (1)Skip section navigation (2)
Date:      Tue, 31 Jul 2012 15:46:20 -0500
From:      Stephen Montgomery-Smith <stephen@missouri.edu>
To:        Bruce Evans <brde@optusnet.com.au>
Cc:        freebsd-numerics@FreeBSD.org
Subject:   Re: Bad error computing 1/z
Message-ID:  <5018441C.5040909@missouri.edu>
In-Reply-To: <20120730131929.Y876@besplex.bde.org>
References:  <5015DA8A.10200@missouri.edu> <20120730025118.GB7958@server.rulingia.com> <5015F9E9.9070408@missouri.edu> <20120730131929.Y876@besplex.bde.org>

next in thread | previous in thread | raw e-mail | index | archive | help
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.



Want to link to this message? Use this URL: <https://mail-archive.FreeBSD.org/cgi/mid.cgi?5018441C.5040909>