mirror of
https://sourceware.org/git/glibc.git
synced 2024-11-23 13:30:06 +00:00
4c327f2ad8
This patch fixes bug 16408, ldbl-128ibm expm1l returning NaN for some large arguments. The basic problem is that the approach of converting the exponent to the form n * log(2) + y, where -0.5 <= y <= 0.5, then computing 2^n * expm1(y) + (2^n - 1) falls over when 2^n overflows (starting slightly before the point where expm1 overflows, when y is negative and n is the least integer for which 2^n overflows). The ldbl-128 code, and the x86/x86_64 code, make expm1l fall back to expl for large positive arguments to avoid this issue. This patch makes the ldbl-128ibm code do the same. (The problem appears for the particular argument in the testsuite because the ldbl-128ibm code also uses an overflow threshold that's for ldbl-128 and is too big for ldbl-128ibm, but the problem described applies for large non-overflowing cases as well, although during the freeze is not a suitable time for making the expm1 tests cover cases close to overflow more thoroughly.) This leaves some code for large positive arguments in expm1l that is now dead. To keep the code for ldbl-128 and ldbl-128ibm similar, and to avoid unnecessary changes during the freeze, the patch doesn't remove it; instead I propose to file a bug in Bugzilla as a reminder that this code (for overflow, including errno setting, and for arguments of +Inf) is no longer needed and should be removed from both those expm1l implementations. Tested powerpc32. * sysdeps/ieee754/ldbl-128ibm/s_expm1l.c (__expm1l): Use __expl for large positive arguments.
167 lines
4.4 KiB
C
167 lines
4.4 KiB
C
/* expm1l.c
|
|
*
|
|
* Exponential function, minus 1
|
|
* 128-bit long double precision
|
|
*
|
|
*
|
|
*
|
|
* SYNOPSIS:
|
|
*
|
|
* long double x, y, expm1l();
|
|
*
|
|
* y = expm1l( x );
|
|
*
|
|
*
|
|
*
|
|
* DESCRIPTION:
|
|
*
|
|
* Returns e (2.71828...) raised to the x power, minus one.
|
|
*
|
|
* Range reduction is accomplished by separating the argument
|
|
* into an integer k and fraction f such that
|
|
*
|
|
* x k f
|
|
* e = 2 e.
|
|
*
|
|
* An expansion x + .5 x^2 + x^3 R(x) approximates exp(f) - 1
|
|
* in the basic range [-0.5 ln 2, 0.5 ln 2].
|
|
*
|
|
*
|
|
* ACCURACY:
|
|
*
|
|
* Relative error:
|
|
* arithmetic domain # trials peak rms
|
|
* IEEE -79,+MAXLOG 100,000 1.7e-34 4.5e-35
|
|
*
|
|
*/
|
|
|
|
/* Copyright 2001 by Stephen L. Moshier
|
|
|
|
This library is free software; you can redistribute it and/or
|
|
modify it under the terms of the GNU Lesser General Public
|
|
License as published by the Free Software Foundation; either
|
|
version 2.1 of the License, or (at your option) any later version.
|
|
|
|
This library is distributed in the hope that it will be useful,
|
|
but WITHOUT ANY WARRANTY; without even the implied warranty of
|
|
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
|
|
Lesser General Public License for more details.
|
|
|
|
You should have received a copy of the GNU Lesser General Public
|
|
License along with this library; if not, see
|
|
<http://www.gnu.org/licenses/>. */
|
|
|
|
#include <errno.h>
|
|
#include <math.h>
|
|
#include <math_private.h>
|
|
#include <math_ldbl_opt.h>
|
|
|
|
/* exp(x) - 1 = x + 0.5 x^2 + x^3 P(x)/Q(x)
|
|
-.5 ln 2 < x < .5 ln 2
|
|
Theoretical peak relative error = 8.1e-36 */
|
|
|
|
static const long double
|
|
P0 = 2.943520915569954073888921213330863757240E8L,
|
|
P1 = -5.722847283900608941516165725053359168840E7L,
|
|
P2 = 8.944630806357575461578107295909719817253E6L,
|
|
P3 = -7.212432713558031519943281748462837065308E5L,
|
|
P4 = 4.578962475841642634225390068461943438441E4L,
|
|
P5 = -1.716772506388927649032068540558788106762E3L,
|
|
P6 = 4.401308817383362136048032038528753151144E1L,
|
|
P7 = -4.888737542888633647784737721812546636240E-1L,
|
|
Q0 = 1.766112549341972444333352727998584753865E9L,
|
|
Q1 = -7.848989743695296475743081255027098295771E8L,
|
|
Q2 = 1.615869009634292424463780387327037251069E8L,
|
|
Q3 = -2.019684072836541751428967854947019415698E7L,
|
|
Q4 = 1.682912729190313538934190635536631941751E6L,
|
|
Q5 = -9.615511549171441430850103489315371768998E4L,
|
|
Q6 = 3.697714952261803935521187272204485251835E3L,
|
|
Q7 = -8.802340681794263968892934703309274564037E1L,
|
|
/* Q8 = 1.000000000000000000000000000000000000000E0 */
|
|
/* C1 + C2 = ln 2 */
|
|
|
|
C1 = 6.93145751953125E-1L,
|
|
C2 = 1.428606820309417232121458176568075500134E-6L,
|
|
/* ln (2^16384 * (1 - 2^-113)) */
|
|
maxlog = 1.1356523406294143949491931077970764891253E4L,
|
|
/* ln 2^-114 */
|
|
minarg = -7.9018778583833765273564461846232128760607E1L, big = 1e290L;
|
|
|
|
|
|
long double
|
|
__expm1l (long double x)
|
|
{
|
|
long double px, qx, xx;
|
|
int32_t ix, lx, sign;
|
|
int k;
|
|
double xhi;
|
|
|
|
/* Detect infinity and NaN. */
|
|
xhi = ldbl_high (x);
|
|
EXTRACT_WORDS (ix, lx, xhi);
|
|
sign = ix & 0x80000000;
|
|
ix &= 0x7fffffff;
|
|
if (!sign && ix >= 0x40600000)
|
|
return __expl (x);
|
|
if (ix >= 0x7ff00000)
|
|
{
|
|
/* Infinity. */
|
|
if (((ix - 0x7ff00000) | lx) == 0)
|
|
{
|
|
if (sign)
|
|
return -1.0L;
|
|
else
|
|
return x;
|
|
}
|
|
/* NaN. No invalid exception. */
|
|
return x;
|
|
}
|
|
|
|
/* expm1(+- 0) = +- 0. */
|
|
if ((ix | lx) == 0)
|
|
return x;
|
|
|
|
/* Overflow. */
|
|
if (x > maxlog)
|
|
{
|
|
__set_errno (ERANGE);
|
|
return (big * big);
|
|
}
|
|
|
|
/* Minimum value. */
|
|
if (x < minarg)
|
|
return (4.0/big - 1.0L);
|
|
|
|
/* Express x = ln 2 (k + remainder), remainder not exceeding 1/2. */
|
|
xx = C1 + C2; /* ln 2. */
|
|
px = __floorl (0.5 + x / xx);
|
|
k = px;
|
|
/* remainder times ln 2 */
|
|
x -= px * C1;
|
|
x -= px * C2;
|
|
|
|
/* Approximate exp(remainder ln 2). */
|
|
px = (((((((P7 * x
|
|
+ P6) * x
|
|
+ P5) * x + P4) * x + P3) * x + P2) * x + P1) * x + P0) * x;
|
|
|
|
qx = (((((((x
|
|
+ Q7) * x
|
|
+ Q6) * x + Q5) * x + Q4) * x + Q3) * x + Q2) * x + Q1) * x + Q0;
|
|
|
|
xx = x * x;
|
|
qx = x + (0.5 * xx + xx * px / qx);
|
|
|
|
/* exp(x) = exp(k ln 2) exp(remainder ln 2) = 2^k exp(remainder ln 2).
|
|
|
|
We have qx = exp(remainder ln 2) - 1, so
|
|
exp(x) - 1 = 2^k (qx + 1) - 1
|
|
= 2^k qx + 2^k - 1. */
|
|
|
|
px = __ldexpl (1.0L, k);
|
|
x = px * qx + (px - 1.0);
|
|
return x;
|
|
}
|
|
libm_hidden_def (__expm1l)
|
|
long_double_symbol (libm, __expm1l, expm1l);
|