Date: Wed, 15 Aug 2012 15:56:24 -0500 From: Stephen Montgomery-Smith <stephen@missouri.edu> To: Bruce Evans <brde@optusnet.com.au> Cc: freebsd-numerics@freebsd.org Subject: Re: Complex arg-trig functions Message-ID: <502C0CF8.8040003@missouri.edu> In-Reply-To: <20120815223631.N1751@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>
next in thread | previous in thread | raw e-mail | index | archive | help
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.
Want to link to this message? Use this URL: <https://mail-archive.FreeBSD.org/cgi/mid.cgi?502C0CF8.8040003>