From owner-freebsd-numerics@FreeBSD.ORG Fri Sep 21 19:05:14 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 C4E40106566B for ; Fri, 21 Sep 2012 19:05:14 +0000 (UTC) (envelope-from brde@optusnet.com.au) Received: from mail03.syd.optusnet.com.au (mail03.syd.optusnet.com.au [211.29.132.184]) by mx1.freebsd.org (Postfix) with ESMTP id 4B3068FC0C for ; Fri, 21 Sep 2012 19:05:14 +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 mail03.syd.optusnet.com.au (8.13.1/8.13.1) with ESMTP id q8LJ5APD011591 (version=TLSv1/SSLv3 cipher=DHE-RSA-AES256-SHA bits=256 verify=NO); Sat, 22 Sep 2012 05:05:12 +1000 Date: Sat, 22 Sep 2012 05:05:10 +1000 (EST) From: Bruce Evans X-X-Sender: bde@besplex.bde.org To: Stephen Montgomery-Smith In-Reply-To: <505C7490.90600@missouri.edu> Message-ID: <20120922042112.E3044@besplex.bde.org> References: <5017111E.6060003@missouri.edu> <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> <20120919010613.T2493@besplex.bde.org> <505BD9B4.8020801@missouri.edu> <20120921172402.W945@besplex.bde.org> <20120921212525.W1732@besplex.bde.org> <505C7490.90600@missouri.edu> MIME-Version: 1.0 Content-Type: TEXT/PLAIN; charset=US-ASCII; format=flowed Cc: 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: Fri, 21 Sep 2012 19:05:15 -0000 On Fri, 21 Sep 2012, Stephen Montgomery-Smith wrote: > I will keep the inexact raising up-front (mostly because I forgot how I did > it earlier). It is working very well. I have minor cleanups to it using a 1-line raise_inexact() macro. > I will still use atanh(fl), and rely on someone else to fix it. (If it is > inexact near 0, it is only a few ULP, and that is good enough for me.) I see you simplified the SQRT_[34]_EPSILON thresholds and fixed their initialization by removing them (the function calls in the static initializers didn't compile). I spent more time on them and found the best values in practice (take the sqrt's accurately and divide them by 2 or 4 instead of your 8), but they are painful to initialize, especiallu for long doubles. A few points of more general interest turned up while debugging this, - atan()'s series is alternating, while the others are not. Alternation causes more cancelation errors. - FOO_EPSILON only applies to 1 side of an addition. E.g., it applies to 1.0 + x where x > 0, for 1.0 + x where x < 0 the size of the corresponding epsilon is half as much. Non-alternation means that the FOO_EPSILON side applies. - the general approximation in cacos(z) and casin(z) is quite good for small z, so larger thresholds for using the special approximations don't affect accuracy much. However, the general approximation in catanh(z) is not so good for small z, so using larger threshold for the special approximation affects accuracy significantly. I found the approximate best point to switch the approximations. - some combination of the previous 3 points means that the switching point is about twice as large relative to the SQRT_N_EPSILON threshold for catan() as for the others (divide by 2 instead of 4). > I'll go ahead and see about the pio2h and pio2l. Please wait for my patch for this. It has all the details for all precisions including pretty-printing the declarations. Or you can do a nearly-global substition of m_pi_2 by pio2_hi + pio2_lo. My other mostly-complete changes: - avoid all the scalb() and related calls. This makes do_hard_work() a bit faster and simpler and real_value_reciprocal() much faster and a bit more complex. - make real_value_reciprocal() handle signs (everything is automatic except for x = inf), and avoid a copysign() after it - a few improvements in comments My other unfinished changes: - figure out if the up-front things in catanh() are best placed there. - decide whether to handle pure real and pure imaginary args specially like I do for both in catanhf(). This interacts with the previous point. - decide whether my old change to remove unnecessary accuracy for the case where ax == 1, ay < FLT_EPSILON in catanh() is correct (you didn't accept it, and maybe other accuracy changes make it extra accuracy more interesting). Bruce