From owner-freebsd-numerics@FreeBSD.ORG Wed Sep 12 10:27:05 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 4674E106566B for ; Wed, 12 Sep 2012 10:27:05 +0000 (UTC) (envelope-from brde@optusnet.com.au) Received: from mail35.syd.optusnet.com.au (mail35.syd.optusnet.com.au [211.29.133.51]) by mx1.freebsd.org (Postfix) with ESMTP id C3C9E8FC08 for ; Wed, 12 Sep 2012 10:27:04 +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 mail35.syd.optusnet.com.au (8.13.1/8.13.1) with ESMTP id q8CAQtJ1019686 (version=TLSv1/SSLv3 cipher=DHE-RSA-AES256-SHA bits=256 verify=NO); Wed, 12 Sep 2012 20:26:56 +1000 Date: Wed, 12 Sep 2012 20:26:54 +1000 (EST) From: Bruce Evans X-X-Sender: bde@besplex.bde.org To: Stephen Montgomery-Smith In-Reply-To: <504D294F.1050901@missouri.edu> Message-ID: <20120912195727.H1078@besplex.bde.org> References: <5017111E.6060003@missouri.edu> <20120804165555.X1231@besplex.bde.org> <501D51D7.1020101@missouri.edu> <20120805030609.R3101@besplex.bde.org> <501D9C36.2040207@missouri.edu> <20120805175106.X3574@besplex.bde.org> <501EC015.3000808@missouri.edu> <20120805191954.GA50379@troutmask.apl.washington.edu> <20120807205725.GA10572@server.rulingia.com> <20120809025220.N4114@besplex.bde.org> <5027F07E.9060409@missouri.edu> <20120814003614.H3692@besplex.bde.org> <50295F5C.6010800@missouri.edu> <20120814072946.S5260@besplex.bde.org> <50297CA5.5010900@missouri.edu> <50297E43.7090309@missouri.edu> <20120814201105.T934@besplex.bde.org> <502A780B.2010106@missouri.edu> <20120815223631.N1751@besplex.bde.org> <502C0CF8.8040003@missouri.edu> <503265E8.3060101@missouri.edu> <5036EFDB.3010304@missouri.edu> <503A9A1A.1050004@missouri.edu> <503ABAC4.3060503@missouri.edu> <20120906200420.H1056@besplex.bde.org> <504D294F.1050901@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: Wed, 12 Sep 2012 10:27:05 -0000 On Sun, 9 Sep 2012, Stephen Montgomery-Smith wrote: > On 09/06/2012 06:42 AM, Bruce Evans wrote: >> I'm back from a holiday. >> >> On Sun, 26 Aug 2012, Stephen Montgomery-Smith wrote: > >> I tested your latest version and found regressions for NaNs. (Also, >> the float version is less accurate now that it doesn't call double >> functions, but that is a progression.) >> 1. casin*() abuses x and y for |x| and |y| and thus loses the sign of >> NaNs. The quick fix is to reload the original values. The correct >> fix is to use separate variables for the absolute values, as is >> done for cacosh*() and catanh*(). > > I now uniformly use ax and ay for |x| and |y| throughout. > >> 2. NaNs are now filtered out early > > It is not correct to filter out the NaNs early. This is because of the way > expressions like cacosh(NaN + I*INF) are defined in C99. The INFs must be > filtered out first. Oops. I still prefer filtering out NaNs first. They already have a special case for one NaN and one +-0 (with +-0 only special for y in cacos()), and can 4 more for one Nan and one Inf (2 combinations in each of 2 functions) just as easily as the Inf cases can have special subcases for NaNs. Here is the Inf case for casinhf(): % if (ax > 1/FLT_EPSILON || ay > 1/FLT_EPSILON) % if (isinf(x) || isinf(y) || (int)(1+tiny)==1) { % if (signbit(x) == 0) % w = clog_for_large_values(z) + M_LN2; Adding M_LN2_ risks problems for NaNs. These were fatal for cacos*() so we avoid them there. By filtering out NaNs first, we avoid them better. M_LN2 is a double constant, so using it pessimizes the float case (it happens not to break it). % else % w = clog_for_large_values(-z) + M_LN2; -z does arithmetic on a complex value so it risks doing bad things for NaNs. At best, it toggles the sign of a NaN. clog_for_large_values() always clears the sign bit of NaNs, and thus gives inconsistent sign handling to cases which don't go through here. Filtering out NaNs first would avoid these cases. clog_for_large_values() is pessimized a bit for the usual case by having to handle NaNs: @ static float complex @ clog_for_large_values(float complex z) @ { @ float x, y; @ float ax, ay, t; @ @ x = crealf(z); @ y = cimagf(z); @ @ @ if (x != x || y != y) @ return (cpackf(logf(hypotf(x, y)), atan2f(y, x))); This is `if (isnan(x) || isnan(y))' spelled fairly optimally. In the usual case, this will be repeated in the caller, with pessimal spelling for the float and long double cases and different fairly optimal spelling for the double case. By filtering out the NaN cases early, we can avoid repeating this. >> % @@ -289,7 +284,16 @@ >> % * the arguments is not NaN, so we opt not to raise it. >> % */ >> % - return (cpack(x+y, x+y)); >> % + return (cpack(x+0.0L+(y+0), x+0.0L+(y+0))); >> % } > > Why does the second way work, and the first way doesn't? Because the result may depend on the ALU, and only the second way is evaluated in the same ALU in all precisions. Bruce