weird behavior with casinh()

[Leo Filippini]

Hello,


I noticed while writing code for am ARM project that the casinh() 
function in newlib behaves oddly for real numbers larger than about 10.


I don't have a working newlib for my desktop machine so I copied the 
code from the source, but I see the same behavior on my ARM project. 
This is a minimum working example:


#include <stdio.h>
#include <complex.h>
#include <math.h>

void printComplex(char * name, double complex num) {
     printf("%s: ", name);
     printf("%.16f %.16f * j\n", creal(num), cimag(num));
}

double complex newlib_casin(double complex z)
{
     double complex w;
     double complex ca, ct, zz, z2;
     double x, y;

     x = creal(z);
     y = cimag(z);

#if 0
     if (y == 0.0) {
         if (fabs(x) > 1.0) {
             w = M_PI_2 + 0.0 * I;
#if 0
             mtherr ("casin", DOMAIN);
#endif
         } else {
             w = asin(x) + 0.0 * I;
         }
         return w;
     }
#endif

/* Power series expansion */
/*
b = cabs(z);
if( b < 0.125 )
{
z2.r = (x - y) * (x + y);
z2.i = 2.0 * x * y;

cn = 1.0;
n = 1.0;
ca.r = x;
ca.i = y;
sum.r = x;
sum.i = y;
do
     {
     ct.r = z2.r * ca.r  -  z2.i * ca.i;
     ct.i = z2.r * ca.i  +  z2.i * ca.r;
     ca.r = ct.r;
     ca.i = ct.i;

     cn *= n;
     n += 1.0;
     cn /= n;
     n += 1.0;
     b = cn/n;

     ct.r *= b;
     ct.i *= b;
     sum.r += ct.r;
     sum.i += ct.i;
     b = fabs(ct.r) + fabs(ct.i);
     }
while( b > MACHEP );
w->r = sum.r;
w->i = sum.i;
return;
}
*/


     ca = x + y * I;
     ct = ca * I;
     /* sqrt( 1 - z*z) */
     /* cmul( &ca, &ca, &zz ) */
     /*x * x  -  y * y */
     zz = (x - y) * (x + y) + (2.0 * x * y) * I;

     zz = 1.0 - creal(zz) - cimag(zz) * I;
     z2 = csqrt(zz);
     //printComplex("\tz2", z2);

     //printComplex("\tct", ct);
     zz = ct + z2;
     //printComplex("\tzz", zz);
     zz = clog(zz);
     //printComplex("\tzz", zz);
     /* multiply by 1/i = -i */
     w = zz * (-1.0 * I);
     return w;
}

double complex newlib_casinh(double complex z) {
     double complex w;

     w = -1.0 * I * newlib_casin(z * I);
     return w;
}

double complex other_casinh(double complex z) {
     return clog(z + csqrt(z * z + 1));
}

int main(void) {

     double complex yc;
     int x;
     for (x = 15; x < 21; x++) {

         printComplex("x: ", x);

         yc = csinh(x);
         printComplex("yc = csinh(x)", yc);
         printf("\n");

         printComplex("glibc casinh(yc)", casinh(yc));
         printComplex("newlib casinh(yc)", newlib_casinh(yc));
         printComplex("other casinh(yc)", other_casinh(yc));
         printf("\n\n");
     }

     }


When I compile this on my linux machine I get:

x: : 15.0000000000000000 0.0000000000000000 * j
yc = csinh(x): 1634508.6862359023652971 0.0000000000000000 * j

glibc casinh(yc): 15.0000000000000000 0.0000000000000000 * j
newlib casinh(yc): 14.9998785788736946 -0.0000000000000000 * j
other casinh(yc): 15.0000000000000000 0.0000000000000000 * j

...

a few more lines edited away for brevity

...


x: : 19.0000000000000000 0.0000000000000000 * j
yc = csinh(x): 89241150.4815936386585236 0.0000000000000000 * j

glibc casinh(yc): 19.0000000000000000 0.0000000000000000 * j
newlib casinh(yc): 18.0218266945585768 -0.0000000000000000 * j
other casinh(yc): 19.0000000000000000 0.0000000000000000 * j


x: : 20.0000000000000000 0.0000000000000000 * j
yc = csinh(x): 242582597.7048951387405396 0.0000000000000000 * j

glibc casinh(yc): 20.0000000000000000 0.0000000000000000 * j
newlib casinh(yc): inf -nan * j
other casinh(yc): 20.0000000000000000 0.0000000000000000 * j


It seems with the input getting larger and larger casinh() looses 
precision until it stops working completely.

Glibc's version of casinh() doesn't misbehave like that, and also a 
simple function other_casinh() that I created with an expression for the 
inverse hyperbolic sine behaves as expected.


Am I doing something wrong with this? are there disadvantages to 
other_casinh() with respect to newlib's implementation of it that I 
don't know about? other_casinh() also seems to result in slightly 
smaller code size when I compile it for ARM cortex M0.


Thank you,

Leo