From owner-freebsd-numerics@FreeBSD.ORG Tue Jul 31 20:46:23 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 522791065670 for ; Tue, 31 Jul 2012 20:46:23 +0000 (UTC) (envelope-from stephen@missouri.edu) Received: from wilberforce.math.missouri.edu (wilberforce.math.missouri.edu [128.206.184.213]) by mx1.freebsd.org (Postfix) with ESMTP id 1A1F48FC19 for ; Tue, 31 Jul 2012 20:46:23 +0000 (UTC) Received: from [128.206.184.213] (wilberforce.math.missouri.edu [128.206.184.213]) by wilberforce.math.missouri.edu (8.14.5/8.14.5) with ESMTP id q6VKkK7b024207; Tue, 31 Jul 2012 15:46:22 -0500 (CDT) (envelope-from stephen@missouri.edu) Message-ID: <5018441C.5040909@missouri.edu> Date: Tue, 31 Jul 2012 15:46:20 -0500 From: Stephen Montgomery-Smith User-Agent: Mozilla/5.0 (X11; FreeBSD amd64; rv:14.0) Gecko/20120728 Thunderbird/14.0 MIME-Version: 1.0 To: Bruce Evans References: <5015DA8A.10200@missouri.edu> <20120730025118.GB7958@server.rulingia.com> <5015F9E9.9070408@missouri.edu> <20120730131929.Y876@besplex.bde.org> In-Reply-To: <20120730131929.Y876@besplex.bde.org> Content-Type: text/plain; charset=ISO-8859-1; format=flowed Content-Transfer-Encoding: 7bit Cc: freebsd-numerics@FreeBSD.org 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: Tue, 31 Jul 2012 20:46:23 -0000 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.