From owner-freebsd-numerics@FreeBSD.ORG Wed Aug 15 20:56:25 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 E5383106566C for ; Wed, 15 Aug 2012 20:56:25 +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 B33788FC14 for ; Wed, 15 Aug 2012 20:56:25 +0000 (UTC) Received: from [127.0.0.1] (wilberforce.math.missouri.edu [128.206.184.213]) by wilberforce.math.missouri.edu (8.14.5/8.14.5) with ESMTP id q7FKuOwp062046; Wed, 15 Aug 2012 15:56:24 -0500 (CDT) (envelope-from stephen@missouri.edu) Message-ID: <502C0CF8.8040003@missouri.edu> Date: Wed, 15 Aug 2012 15:56:24 -0500 From: Stephen Montgomery-Smith User-Agent: Mozilla/5.0 (X11; Linux i686; rv:14.0) Gecko/20120714 Thunderbird/14.0 MIME-Version: 1.0 To: Bruce Evans 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> In-Reply-To: <20120815223631.N1751@besplex.bde.org> Content-Type: text/plain; charset=ISO-8859-1; format=flowed Content-Transfer-Encoding: 7bit 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, 15 Aug 2012 20:56:26 -0000 On 08/15/2012 08:35 AM, Bruce Evans wrote: > On Tue, 14 Aug 2012, Stephen Montgomery-Smith wrote: >> It seemed to me that there is a logic behind why the the infs and nans >> produce the results they do. I noticed that do_the_hard_work() >> already got the answers correct for the real part *rx. Getting the >> imaginary part to work as well seemed to me to be the cleanest way to >> make it work. (I added all the nan and inf checking after writing the >> rest of the code.) > > 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. 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 > 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))); > > 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. > > % diff -c2 catrigf.c~ catrigf.c > % *** catrigf.c~ Sun Aug 12 17:00:52 2012 > % --- catrigf.c Wed Aug 15 11:57:08 2012 > % *************** > % *** 349,353 **** > % % if (isnan(x) || isnan(y)) > % ! return (cpackf(x+y, x+y)); > % % if (isinf(x) || isinf(y)) > % --- 351,355 ---- > % % if (isnan(x) || isnan(y)) > % ! return (cpack((x+0.0L)+(y+0), (x+0.0L)+(y+0))); > % % if (isinf(x) || isinf(y)) > % diff -c2 catrigl.c~ catrigl.c > % *** catrigl.c~ Sun Aug 12 06:54:46 2012 > % --- catrigl.c Wed Aug 15 11:58:46 2012 > % *************** > % *** 323,327 **** > % % if (isnan(x) || isnan(y)) > % ! return (cpackl(x+y, x+y)); > % % if (isinf(x) || isinf(y)) > % --- 325,329 ---- > % % if (isnan(x) || isnan(y)) > % ! return (cpack((x+0.0L)+(y+0), (x+0.0L)+(y+0))); > % % if (isinf(x) || isinf(y)) > % 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.) The part that follows - is this all referencing csqrt? > > % - > % -/* 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 > 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.