From owner-freebsd-numerics@FreeBSD.ORG Thu Sep 6 13:51:08 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 9DB77106566C for ; Thu, 6 Sep 2012 13:51:08 +0000 (UTC) (envelope-from brde@optusnet.com.au) Received: from mail07.syd.optusnet.com.au (mail07.syd.optusnet.com.au [211.29.132.188]) by mx1.freebsd.org (Postfix) with ESMTP id 20BB88FC1E for ; Thu, 6 Sep 2012 13:51:07 +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 mail07.syd.optusnet.com.au (8.13.1/8.13.1) with ESMTP id q86DowsF008413 (version=TLSv1/SSLv3 cipher=DHE-RSA-AES256-SHA bits=256 verify=NO); Thu, 6 Sep 2012 23:51:00 +1000 Date: Thu, 6 Sep 2012 23:50:58 +1000 (EST) From: Bruce Evans X-X-Sender: bde@besplex.bde.org To: Stephen Montgomery-Smith In-Reply-To: <502C0CF8.8040003@missouri.edu> Message-ID: <20120906221028.O1542@besplex.bde.org> References: <5017111E.6060003@missouri.edu> <501C361D.4010807@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> 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: Thu, 06 Sep 2012 13:51:08 -0000 On Wed, 15 Aug 2012, Stephen Montgomery-Smith wrote: > On 08/15/2012 08:35 AM, Bruce Evans wrote: >> An up-front check may still be simpler, and gives more control. In >> csqrt*(), I needed an explicit check and special expressions to get >> uniform behaviour. > > I still like it the way I have it. There is a definite logic in the way infs > and nans come out of casinh, etc. Changed now :-). > There is only one place I disagree with C99: > catanh(1) = Inf + 0*I > > I think mpc gets it correct: > atanh(1) = Inf + nan*I C99 is trying be compatible with the real function to a fault. Real atanh(1) must be Inf, and any imaginary part would make catanh(1) different. C99 requires similar beahaviour for many other functions. >> I added this to the NaN mixing in catan[h]*(), >> and now all my tests pass: >> >> % diff -c2 catrig.c~ catrig.c >> % *** catrig.c~ Sun Aug 12 17:29:18 2012 >> % --- catrig.c Wed Aug 15 11:57:02 2012 >> % *************** >> % *** 605,609 **** >> % */ >> % if (ISNAN(x) || ISNAN(y)) >> % ! return (cpack(x+y, x+y)); >> % % /* (x,inf) and (inf,y) and (inf,inf) -> (0,PI/2) */ >> % --- 609,613 ---- >> % */ >> % if (ISNAN(x) || ISNAN(y)) >> % ! return (cpack((x+0.0L)+(y+0), (x+0.0L)+(y+0))); >> % % /* (x,inf) and (inf,y) and (inf,inf) -> (0,PI/2) */ >> >> Use this expression in all precisions. > > Would this work? > > if (ISNAN(x) || ISNAN(y)) > return (cpack((x+x)+(y+y), (x+x)+(y+y))); No, since this would be evaluated by a different ALU than my version, if the ALU for long doubles is different from the ALU for doubles, as happens on amd64. However, I could use y+y instead of y+0 in some cases. The point of using y+0 is to use the same code and not change the result much if y happens to be finite (only -0 gets changed, to 0). ALUs: (x+0.0L)+(y+0): x+0.0L: long double + long double, after promotion of x. Always in the i387. y+0 (or y+y): typeof(y) + typeof(y), after promotion of x. In the i387 on i386, but in SSE on amd64. y+0 ends up quasi-promoted on i386. final addition: long double + long double, after promotion of y+0. Always in the i387. Always long double result. Returning this then converts to double. (x+x)+(y+y): x+x: as above for y+y. Only in the i387 on i386. y+y: as above final addition: still only in the i387 on i386. So the result of casin() on 2 NaNs is quite different from that of casinl() on the same 2 NaNs (initially), since the i387 ALU has different semantics for NaNs than the SSE ALU despite the ALUs sharing the FPU hardware. Only quieting of NaNs works perfectly in this version. >> I forgot to comment it. Adding 0 quietens signaling NaNs before mixing >> NaNs. I should have tried y+y. Adding 0.0L promotes part of the >> expression to long double together with quietening signaling NaNs. >> The rest of the expression is promoted to match. I should try the >> old way again: of (long double)x+x. Another explanation. Without testing it, I think it has the same speed with the x+0[.0L] versions, but is slightly better for the x+x version, since x+x can be done without another load for 0 in SSE only. But the extra load is likely to free on x86 since it can be hidden in scheduling latencies. >> ... >> % Index: ../s_csqrt.c >> % =================================================================== >> % RCS file: /home/ncvs/src/lib/msun/src/s_csqrt.c,v >> % retrieving revision 1.4 >> % diff -u -2 -r1.4 s_csqrt.c >> % --- ../s_csqrt.c 8 Aug 2008 00:15:16 -0000 1.4 >> % +++ ../s_csqrt.c 14 Aug 2012 20:34:07 -0000 >> % @@ -34,14 +34,5 @@ >> % #include "math_private.h" >> % % -/* >> % - * gcc doesn't implement complex multiplication or division correctly, >> % - * so we need to handle infinities specially. We turn on this pragma to >> % - * notify conforming c99 compilers that the fast-but-incorrect code that >> % - * gcc generates is acceptable, since the special cases have already >> been >> % - * handled. >> % - */ >> % -#pragma STDC CX_LIMITED_RANGE ON >> >> Remove this. There was only 1 complex expression, and it depended on the >> negation of this pragma to work. Since gcc doesn't support this pragma, >> the expression only worked accidentally when it was optimized. > > I removed it. (I copied it verbatim from csqrt without really understanding > it.) Good. > The part that follows - is this all referencing csqrt? Depends on whether you clipped an "Index" header :-). >> >> % - >> % -/* We risk spurious overflow for components >= DBL_MAX / (1 + >> sqrt(2)). */ >> % +/* For avoiding overflow for components >= DBL_MAX / (1 + sqrt(2)). */ >> % #define THRESH 0x1.a827999fcef32p+1022 >> >> .............. snip Looks like csqrt. I didn't check that THRESH is actually correct. It seems to have more bits than are strictly necessary. >> This is like a fix in clog(). hypot() handles denormals OK, but >> necessarily loses accuracy when it returns a denormal result, so >> the expression (a + hypot(a, b)) is more inaccurate than necessary. > > Which code is being referenced here? I use expressions like this catrig. > Although I think when I use it, I am somewhat certain that neither a nor b > are denormal. Now "This" is the change in csqrt() in the snipped quote and is being compared to related code in clog() that wasn't quoted recently in this thread. Snipped quote: % + scale = 2; % } else { % - scale = 0; % + scale = 1; % + } % + % + /* Scale to reduce inaccuracies when both components are denormal. */ % + if (fabs(a) <= 0x1p-1023 && fabs(b) <= 0x1p-1023) { % + a *= 0x1p54; % + b *= 0x1p54; % + scale = 0x1p-27; Related code in clog(): % /* Reduce inaccuracies and avoid underflow when ax is denormal. */ % if (kx <= MIN_EXP - 2) % return (cpack(log(hypot(x * 0x1p1023, y * 0x1p1023)) + % (MIN_EXP - 2) * ln2_lo + (MIN_EXP - 2) * ln2_hi, v)); % % /* Avoid remaining underflows (when ax is small but not denormal). */ % if (ky < (MIN_EXP - 1) / 2 + MANT_DIG) % return (cpack(log(hypot(x, y)), v)); This code in catrig.c seems to be related: @ /* @ * sum_squares(x,y) = x*x + y*y. @ * Assumes x*x and y*y will not overflow. @ * Assumes x and y are finite. @ * Assumes fabs(x) >= DBL_EPSILON. @ */ @ inline static double @ sum_squares(double x, double y) @ { @ /* @ * The following code seems to do better than @ * return (hypot(x, y) * hypot(x, y)); @ */ @ @ if (fabs(y) < MIN_4TH_ROOT) @ if ((int)y==0) /* raise inexact */ @ return (x*x); @ return (x*x + y*y); @ } but it isn't really -- it is more like the case of a large difference of exponents, so x*x+y*y reduces to x*x (at a smaller exponent difference than for hypot(x, y) reducing to |x|). hypot() sets inexact by returning |x|+|y| in that case, but the above avoids using y*y since it would give spurious underflow. Hmm, MIN_4TH_ROOT seems to be logically wrong here, and physically wrong for float precision. In double precision, MIN_4TH_ROOT/DBL_EPSILON is ~0.5e-59, but in float precision, MIN_4TH_ROOT/FLT_EPSILON is only ~1e-2, so when x = FLT_EPSILON, y = MIN_4TH_ROOT is very far from giving a discardable y*y. The logically correct threshold is something like y/x < 2**-MANT_DIG/2 or y < 2**-MANT_DIG/2 * EPSILON to avoid the division. clog*() uses essentially y/x < 2**-MANT_DIG for some reason (it actually uses a fuzzy check based on exponent differences: kx - ky > MANT_DIG; MANT_DIG/2 would be a little too small here, but doubling it seems wasteful). hypot() uses essentially y/x < 2**-60, where the magic 60 is DBL_MANT_DIG plys a safety margin. hypotf() uses a magic 30. The non-magic integer exponent thresholds work very well in clog*(), but depend on having the exponents handy for efficiency, which may be possible via optimization for an inline function like the above but takes more source code. I don't see how the above avoids cancelation problems like the ones in clog(). You never subtract 1 from sum_squares(), but there are open-coded expressions like '(1-ax)*(1+ax) - ay*ay' to reduce the cancelation error. I found that for clog() it was worth using extra precision out to about |z| = 4 (or |z| < 1/4) to reduce the cancelation error. But somehow my tests didn't show large errors. Perhaps because I'm happy with < 4 ulps instead of < 1. Bruce