weird behavior with casinh()

[Paul Zimmermann]

       Hi,

indeed, with the test program we used for our comparison [1],
I obtain with Newlib 4.4.0:

casinh 0 -1 0x0p+0,0x1.0000002d413cdp+26 [inf] inf inf
libm gives (inf,-nan)
mpc gives (0x1.2b7088750621fp+4,0x1.921fb54442d18p+0)

i.e., for z=(0x0p+0,0x1.0000002d413cdp+26), the casinh() function from
Newlib gives (Inf,NaN) instead of (0x1.2b7088750621fp+4,0x1.921fb54442d18p+0).

Paul

[1] https://inria.hal.science/hal-04714173

> Date: Mon,  7 Oct 2024 00:29:24 +0000
> From: Leo Filippini <leo.filippini@posteo.org>
> 
> 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
> 
> 
>