1.1b: asinh breaks at -O1

[Colin Douglas Howell]

Jeffrey A Law writes:
> 
>   In message < 19981001222047.A19471@zensunni >you write:
>   > [This problem has been reported as a bug through the Debian bugtracking
>   > system ( http://www.debian.org/Bugs/ ); please Cc: 27146@bugs.debian.org on
>   > on-topic replies]
>   > 
>   > OS: Linux 2.0.35
>   > Architecture: i686
>   > egcs: 1.1b
>   > libc: GNU libc 2.0.7 (aka libc6)
>   > 
>   > Code test.c:
>   > #include <stdlib.h>
>   > #include <math.h>
>   > 
>   > int main()
>   > {
>   >   double z = 3;
>   > 
>   >   printf("glib result asinh(%f)=%f\n", -z, asinh(-z));
>   >   printf("glib result asinh(%f)=%f\n", z, asinh(z));
>   >   exit(0);
>   > }
>   > 
>   > Behaviour: 
>   > 
>   > egcc -O0 test.c -lm ; ./a.out # As expected
>   > glib result asinh(-3.000000)=-1.818446
>   > glib result asinh(3.000000)=1.818446
>   > 
>   > egcc -O1 test.c -lm ; ./a.out # bug: asinh misbehaves with negative argumen
>   > t
>   > glib result asinh(-3.000000)=NaN
>   > glib result asinh(3.000000)=1.818446
> Can someone who either a.  Knows the x86 or b. has a debugger which
> will print the FP variables under linux take a look at this?
> 
> I'd do it, but I don't know the x86 well enough to just look at the
> assembly code and figure out what's going on.  And for reasons I
> don't want to know, gdb won't access the floating point regs on a
> linux box, so I can't run it under the debugger to find out why
> we're losing.

My gdb acts the same way.  Bummer.

I think I've figured it out, but looking at the x86 assembly code
actually wasn't that much help, since the code is correct.  The
problem is glibc 2.0.7's inline definition of asinh() in
/usr/include/__math.h.  This definition has a misplaced paren which
breaks it for negative arguments.  The inline definition is not
included when compiling under gcc 2.7.x, but is included when running
newer versions of gcc or egcs.  Since the inline definition is only
used when optimization is turned on, and since the actual library
definition of asinh() in glibc is apparently correct, the bug only
shows up when optimizing.

Here's glibc 2.0.7's definition of asinh() from /usr/include/__math.h:

__MATH_INLINE double asinh (double __x);
__MATH_INLINE double
asinh (double __x)
{
  register double __y = fabs (__x);

  return log1p ((__y * __y / (sqrt (__y * __y + 1.0) + 1.0) + __y)
                * __sgn1 (__x));
}

Here's a correct definition:

__MATH_INLINE double asinh (double __x);
__MATH_INLINE double
asinh (double __x)
{
  register double __y = fabs (__x);

  return log1p ((__y * __y / (sqrt (__y * __y + 1.0) + 1.0) + __y))
                * __sgn1 (__x);
}

The __sgn1() function (which computes the sign of its argument) should
go outside the log1p(), not inside it.  The argument of log1p() is
supposed to be positive here.  The __sgn1() is meant to invert the
sign of the result for negative arguments of asinh(), since asinh(-x)
= -asinh(x).

(Actually this correction gives an extra pair of parens around the
argument of log1p(), which could be removed.)

-- 
Colin Douglas Howell                    Systems Administrator
e-mail:  howell@cs.stanford.edu         Computer Facilities Group
office:  (650) 723-2491                 Computer Science Department
                                        Stanford University