Skip site navigation (1)Skip section navigation (2)
Date:      Wed, 1 Aug 2012 18:12:59 +1000 (EST)
From:      Bruce Evans <brde@optusnet.com.au>
To:        Stephen Montgomery-Smith <stephen@missouri.edu>
Cc:        freebsd-numerics@FreeBSD.org, Bruce Evans <brde@optusnet.com.au>
Subject:   Re: Bad error computing 1/z
Message-ID:  <20120801180703.B1153@besplex.bde.org>
In-Reply-To: <5018441C.5040909@missouri.edu>
References:  <5015DA8A.10200@missouri.edu> <20120730025118.GB7958@server.rulingia.com> <5015F9E9.9070408@missouri.edu> <20120730131929.Y876@besplex.bde.org> <5018441C.5040909@missouri.edu>

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



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