Re: asinh() gives NaN on Linux/x86/glibc with optimization on

[Jeffrey A Law <law at hurl dot cygnus dot com>]

  In message <75ph1i$phv@gap.cco.caltech.edu>you write:
  > H.J. Lu <hjl@lucon.org> writes:
  > [me saying asinh of negative number gives wrong answer]
  > >
  > ># gcc -v
  > >Reading specs from /usr/lib/gcc-lib/i586-pc-linux/egcs-2.91.60/specs
  > >gcc version egcs-2.91.60 19981201 (egcs-1.1.1 release)
  > ># gcc -O6 x.c
  > ># a.out
  > >Input argument for asinh: -100
  > >asinh(-100.000000)=-5.298342
  > ># ldd a.out
  > >        libc.so.6 => /lib/libc.so.6 (0x40007000)
  > >	/lib/ld-linux.so.2 => /lib/ld-linux.so.2 (0x2aaaa000)
  > ># ls -l /lib/libc.so.6
  > >lrwxrwxrwx   1 root     root           26 Dec 11 08:16 /lib/libc.so.6 ->
  > >glibc-2.0-981211/libc.so.6
  > >
  > >I am running glibc 2.0.7 981211.
  > 
  > Hmm. That's funny. That's exactly what I did and I get NaN's for
  > answer. I'm using the version of glibc that came with RedHat 5.2, so
  > maybe that's the problem. OTOH, what kind of a bug in the library
  > would cause different optimizations in the *calling* program to give
  > different results?
  > 
  > BTW, where do I get the new snapshots of glibc? The snapshot
  > repository in ftp.kernel.org doesn't seem to be updated anymore.
I presume you're running glibc-2.0.7?

  Here's an analysis of an earlier problem with asinh with glibc-2.0.7:
  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.)