[Bug math/5794] New: expm1l(x) overflows when 11356.1768 < x <= 11356.5234

["petrosyan at gmail dot com" <sourceware-bugzilla at sourceware dot org>]

The value of (e^x - 1) is representable as a long double when x <=
11356.5234, but expm1l(x) overflows to infinity when x > 11356.1768
(approximately).

Steps to Reproduce:
Compile and run the following program:

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

main ()
{
  long double a, b;
  const long double x = 11356.35L;

  a = expl (x) - 1.0L;
  b = expm1l (x);
  if (a == b) {
    printf ("e^%.2Lf - 1 = %Lg\n", b);
    return 0;
  } else {
    printf ("e^%.2Lf - 1 = %Lg; but expm1l(%.2Lf) returned %Lg\n",
	    x, a, x, b);
    return 1;
  }
}


Actual Results:  e^11356.35 - 1 = 1.00032e+4932; but expm1l(11356.35)
returned inf


Expected Results:  e^11356.35 - 1 = 1.00032e+4932


Additional info:

This bug probably isn't even worth reporting, except for completeness. 
expm1l is designed to compute e^x-1 more accurately when x is negative, so
[at least the i386 version] uses a slightly different algorithm than expl. 
This algorithm overflows for large x where expl would not.  But then, expm1
really isn't intended for large x anyway.

In my own port of the math library, I circumvent the problem by calling
__ieee754_expl when x >= 64, because at that magnitude, -1 is too small to
affect the result.  This is probably overkill, but it works.

-- 
           Summary: expm1l(x) overflows when 11356.1768 < x <= 11356.5234
           Product: glibc
           Version: unspecified
            Status: NEW
          Severity: minor
          Priority: P2
         Component: math
        AssignedTo: aj at suse dot de
        ReportedBy: petrosyan at gmail dot com
                CC: glibc-bugs at sources dot redhat dot com


http://sourceware.org/bugzilla/show_bug.cgi?id=5794

------- You are receiving this mail because: -------
You are on the CC list for the bug, or are watching someone who is.