From owner-freebsd-numerics@FreeBSD.ORG Wed Aug 1 08:13:04 2012 Return-Path: Delivered-To: freebsd-numerics@FreeBSD.org Received: from mx1.freebsd.org (mx1.freebsd.org [69.147.83.52]) by hub.freebsd.org (Postfix) with ESMTP id D3EFA10656AA for ; Wed, 1 Aug 2012 08:13:04 +0000 (UTC) (envelope-from brde@optusnet.com.au) Received: from mail06.syd.optusnet.com.au (mail06.syd.optusnet.com.au [211.29.132.187]) by mx1.freebsd.org (Postfix) with ESMTP id 520238FC0C for ; Wed, 1 Aug 2012 08:13:03 +0000 (UTC) Received: from c122-106-171-246.carlnfd1.nsw.optusnet.com.au (c122-106-171-246.carlnfd1.nsw.optusnet.com.au [122.106.171.246]) by mail06.syd.optusnet.com.au (8.13.1/8.13.1) with ESMTP id q718D0sh030560 (version=TLSv1/SSLv3 cipher=DHE-RSA-AES256-SHA bits=256 verify=NO); Wed, 1 Aug 2012 18:13:01 +1000 Date: Wed, 1 Aug 2012 18:12:59 +1000 (EST) From: Bruce Evans X-X-Sender: bde@besplex.bde.org To: Stephen Montgomery-Smith In-Reply-To: <5018441C.5040909@missouri.edu> Message-ID: <20120801180703.B1153@besplex.bde.org> References: <5015DA8A.10200@missouri.edu> <20120730025118.GB7958@server.rulingia.com> <5015F9E9.9070408@missouri.edu> <20120730131929.Y876@besplex.bde.org> <5018441C.5040909@missouri.edu> MIME-Version: 1.0 Content-Type: TEXT/PLAIN; charset=US-ASCII; format=flowed Cc: freebsd-numerics@FreeBSD.org, Bruce Evans Subject: Re: Bad error computing 1/z X-BeenThere: freebsd-numerics@freebsd.org X-Mailman-Version: 2.1.5 Precedence: list List-Id: "Discussions of high quality implementation of libm functions." List-Unsubscribe: , List-Archive: List-Post: List-Help: List-Subscribe: , X-List-Received-Date: Wed, 01 Aug 2012 08:13:04 -0000 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