From owner-freebsd-numerics@FreeBSD.ORG Tue Sep 18 15:19:16 2012 Return-Path: Delivered-To: freebsd-numerics@FreeBSD.org Received: from mx1.freebsd.org (mx1.freebsd.org [IPv6:2001:4f8:fff6::34]) by hub.freebsd.org (Postfix) with ESMTP id 94BAC106566C for ; Tue, 18 Sep 2012 15:19:16 +0000 (UTC) (envelope-from brde@optusnet.com.au) Received: from mail02.syd.optusnet.com.au (mail02.syd.optusnet.com.au [211.29.132.183]) by mx1.freebsd.org (Postfix) with ESMTP id 2285B8FC12 for ; Tue, 18 Sep 2012 15:19:15 +0000 (UTC) Received: from c122-106-157-84.carlnfd1.nsw.optusnet.com.au (c122-106-157-84.carlnfd1.nsw.optusnet.com.au [122.106.157.84]) by mail02.syd.optusnet.com.au (8.13.1/8.13.1) with ESMTP id q8IFJC2E008798 (version=TLSv1/SSLv3 cipher=DHE-RSA-AES256-SHA bits=256 verify=NO); Wed, 19 Sep 2012 01:19:14 +1000 Date: Wed, 19 Sep 2012 01:19:12 +1000 (EST) From: Bruce Evans X-X-Sender: bde@besplex.bde.org To: Bruce Evans In-Reply-To: <20120918232850.N2144@besplex.bde.org> Message-ID: <20120919010613.T2493@besplex.bde.org> References: <5017111E.6060003@missouri.edu> <20120912225847.J1771@besplex.bde.org> <50511B40.3070009@missouri.edu> <20120913204808.T1964@besplex.bde.org> <5051F59C.6000603@missouri.edu> <20120914014208.I2862@besplex.bde.org> <50526050.2070303@missouri.edu> <20120914212403.H1983@besplex.bde.org> <50538E28.6050400@missouri.edu> <20120915231032.C2669@besplex.bde.org> <50548E15.3010405@missouri.edu> <5054C027.2040008@missouri.edu> <5054C200.7090307@missouri.edu> <20120916041132.D6344@besplex.bde.org> <50553424.2080902@missouri.edu> <20120916134730.Y957@besplex.bde.org> <5055ECA8.2080008@missouri.edu> <20120917022614.R2943@besplex.bde.org> <50562213.9020400@missouri.edu> <20120917060116.G3825@besplex.bde.org> <50563C57.60806@missouri.edu> <20120918012459.V5094@besplex.bde.org> <5057A932.3000603@missouri.edu> <5057F24B.7020605@missouri.edu> <20120918162105.U991@besplex.bde.org> <20120918232850.N2144@besplex.bde.org> MIME-Version: 1.0 Content-Type: TEXT/PLAIN; charset=US-ASCII; format=flowed Cc: Stephen Montgomery-Smith , freebsd-numerics@FreeBSD.org Subject: Re: Complex arg-trig functions 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, 18 Sep 2012 15:19:16 -0000 On Wed, 19 Sep 2012, I wrote: > @ + /* > @ + * Handle the annoying special case +-1 + I*+-0, and collaterally > @ + * handle the not-so-special case y == 0. C99 specifies that > @ + * catanh(+-1 + I*+-0) = +-Inf + I*+-0 instead of the limiting > @ + * value +-Inf + I*+-PI/2 since it wants y == 0 to give the same > @ + * result as the real atanh() (at least for y == +0). The special > @ + * behaviour for +-1 + I*+-0 begins with classifying it to avoid > @ + * raising inexact for it. Make the classification as simple and > @ + * short as possible (except for this comment about it) and ensure > @ + * identical results by calling the real atanh() for all non-NaN x > @ + * when y == 0. This turns out to be significantly more accurate. > @ + * > @ + * TODO: move this before the NaN classification and let atanh() > @ + * handle NaN x too. Make a similar special case for x == 0 to > @ + * improve accuracy; this takes no extra lines of code since it > @ + * removes the need to handle x == 0 under the NaN classification. > @ + */ > @ + if (y == 0) > @ + return (cpackf(atanh(x), y)); > > See the comment. Duh, this has to be under (y == 0 && ax <= 1) so that the real function actually applies. > @ @ - if (ax == 1 && ay < FLT_EPSILON) { > @ - if ((int)ay==0) { > @ - if ( ilogbf(ay) > FLT_MIN_EXP) > @ - rx = - logf(ay/2) / 2; > @ - else > @ - rx = - (logf(ay) - m_ln2) / 2; > @ - } > @ - } else > > Inexact was raised up front, and will also be raised by logf(). > > ilogbf() is rather slow, though it is now a builtin. There is no need > to use it here: > - the condition can be written as (ay > 2 * FLT_MIN_EXP) > - the expression with m_ln2 is accurate enough, so there is no need for > the condition. I thought I tested this assertion and found no difference > at all in accuracy, but now I can't see why it is true. This case is > fundamentally quite accurate -- within about 1 ulp for the logf() part -- > and subtracting m_ln2 will only lose about 0.5 ulps pf accuracy (since > |logf(FLT_EPSILON)| dominates m_ln2), so it is not near the worse case > for accuracy, but the loss of accuracy is not null. Now tested. The increase in inaccuracy is only from ~0.7 ulps to ~0.9 ulps. This is acceptable. > @ ... > @ + > @ + if (ax == 1 && ay < FLT_EPSILON) > @ + rx = - (logf(ay) - m_ln2) / 2; However, the outer FLT_EPSILON threshold for the above is too conservative. It can be increased to FLT_EPSILON**2 without expanding the error to above 0.7 ulps, provided this optimization is not used -- with the exanded threshold, this optimization expands the error by another 0.2 ulps, to ~1.1 ulps instead of to ~0.9 ulps. These errors are still in the noise compared with the worst case error of ~2.6 ulps, but it is good to keep errors nelow 1 ulp if this is easy. Bruce