Skip site navigation (1)Skip section navigation (2)
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>