Complex arg-trig functions

[Stephen Montgomery-Smith]

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.