From mboxrd@z Thu Jan  1 00:00:00 1970
Return-Path: <glibc-bugs-return-7817-listarch-glibc-bugs=sources.redhat.com@sourceware.org>
Received: (qmail 32724 invoked by alias); 26 Feb 2008 07:04:31 -0000
Received: (qmail 32532 invoked by uid 48); 26 Feb 2008 07:03:51 -0000
Date: Tue, 26 Feb 2008 07:04:00 -0000
From: "petrosyan at gmail dot com" <sourceware-bugzilla@sourceware.org>
To: glibc-bugs@sources.redhat.com
Message-ID: <20080226070350.5794.petrosyan@gmail.com>
Reply-To: sourceware-bugzilla@sourceware.org
Subject: [Bug math/5794] New: expm1l(x) overflows when 11356.1768 < x <= 11356.5234
X-Bugzilla-Reason: CC
Mailing-List: contact glibc-bugs-help@sourceware.org; run by ezmlm
Precedence: bulk
List-Id: <glibc-bugs.sourceware.org>
List-Subscribe: <mailto:glibc-bugs-subscribe@sourceware.org>
List-Post: <mailto:glibc-bugs@sourceware.org>
List-Help: <mailto:glibc-bugs-help@sourceware.org>, <http://sourceware.org/lists.html#faqs>
Sender: glibc-bugs-owner@sourceware.org
X-SW-Source: 2008-02/txt/msg00112.txt.bz2

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.