mirror of
https://sourceware.org/git/glibc.git
synced 2024-11-23 21:40:12 +00:00
03d95bd483
This patch fixes the remaining part of bug 16560, spurious underflows from exp2 of arguments close to 0 (when the result is close to 1, so should not underflow), by just using 1+x instead of a more complicated calculation when the argument is sufficiently small. Tested for x86_64, x86 and mips64. [BZ #16560] * math/e_exp2l.c [LDBL_MANT_DIG == 106] (LDBL_EPSILON): Undefine and redefine. (__ieee754_exp2l): Do not multiply small fractional parts by M_LN2l. * sysdeps/i386/fpu/e_exp2l.S (__ieee754_exp2l): Just add 1 to small argument. * sysdeps/ieee754/dbl-64/e_exp2.c (__ieee754_exp2): Likewise. * sysdeps/ieee754/flt-32/e_exp2f.c (__ieee754_exp2f): Likewise. * sysdeps/x86_64/fpu/e_exp2l.S (__ieee754_exp2l): Likewise. * math/auto-libm-test-in: Add more tests of exp2. * math/auto-libm-test-out: Regenerated.
128 lines
3.8 KiB
C
128 lines
3.8 KiB
C
/* Double-precision floating point 2^x.
|
|
Copyright (C) 1997-2015 Free Software Foundation, Inc.
|
|
This file is part of the GNU C Library.
|
|
Contributed by Geoffrey Keating <geoffk@ozemail.com.au>
|
|
|
|
The GNU C 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.
|
|
|
|
The GNU C 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 the GNU C Library; if not, see
|
|
<http://www.gnu.org/licenses/>. */
|
|
|
|
/* The basic design here is from
|
|
Shmuel Gal and Boris Bachelis, "An Accurate Elementary Mathematical
|
|
Library for the IEEE Floating Point Standard", ACM Trans. Math. Soft.,
|
|
17 (1), March 1991, pp. 26-45.
|
|
It has been slightly modified to compute 2^x instead of e^x.
|
|
*/
|
|
#include <stdlib.h>
|
|
#include <float.h>
|
|
#include <ieee754.h>
|
|
#include <math.h>
|
|
#include <fenv.h>
|
|
#include <inttypes.h>
|
|
#include <math_private.h>
|
|
|
|
#include "t_exp2.h"
|
|
|
|
static const double TWO1023 = 8.988465674311579539e+307;
|
|
static const double TWOM1000 = 9.3326361850321887899e-302;
|
|
|
|
double
|
|
__ieee754_exp2 (double x)
|
|
{
|
|
static const double himark = (double) DBL_MAX_EXP;
|
|
static const double lomark = (double) (DBL_MIN_EXP - DBL_MANT_DIG - 1);
|
|
|
|
/* Check for usual case. */
|
|
if (__glibc_likely (isless (x, himark)))
|
|
{
|
|
/* Exceptional cases: */
|
|
if (__glibc_unlikely (!isgreaterequal (x, lomark)))
|
|
{
|
|
if (__isinf (x))
|
|
/* e^-inf == 0, with no error. */
|
|
return 0;
|
|
else
|
|
/* Underflow */
|
|
return TWOM1000 * TWOM1000;
|
|
}
|
|
|
|
static const double THREEp42 = 13194139533312.0;
|
|
int tval, unsafe;
|
|
double rx, x22, result;
|
|
union ieee754_double ex2_u, scale_u;
|
|
|
|
if (fabs (x) < DBL_EPSILON / 4.0)
|
|
return 1.0 + x;
|
|
|
|
{
|
|
SET_RESTORE_ROUND_NOEX (FE_TONEAREST);
|
|
|
|
/* 1. Argument reduction.
|
|
Choose integers ex, -256 <= t < 256, and some real
|
|
-1/1024 <= x1 <= 1024 so that
|
|
x = ex + t/512 + x1.
|
|
|
|
First, calculate rx = ex + t/512. */
|
|
rx = x + THREEp42;
|
|
rx -= THREEp42;
|
|
x -= rx; /* Compute x=x1. */
|
|
/* Compute tval = (ex*512 + t)+256.
|
|
Now, t = (tval mod 512)-256 and ex=tval/512 [that's mod, NOT %;
|
|
and /-round-to-nearest not the usual c integer /]. */
|
|
tval = (int) (rx * 512.0 + 256.0);
|
|
|
|
/* 2. Adjust for accurate table entry.
|
|
Find e so that
|
|
x = ex + t/512 + e + x2
|
|
where -1e6 < e < 1e6, and
|
|
(double)(2^(t/512+e))
|
|
is accurate to one part in 2^-64. */
|
|
|
|
/* 'tval & 511' is the same as 'tval%512' except that it's always
|
|
positive.
|
|
Compute x = x2. */
|
|
x -= exp2_deltatable[tval & 511];
|
|
|
|
/* 3. Compute ex2 = 2^(t/512+e+ex). */
|
|
ex2_u.d = exp2_accuratetable[tval & 511];
|
|
tval >>= 9;
|
|
unsafe = abs (tval) >= -DBL_MIN_EXP - 1;
|
|
ex2_u.ieee.exponent += tval >> unsafe;
|
|
scale_u.d = 1.0;
|
|
scale_u.ieee.exponent += tval - (tval >> unsafe);
|
|
|
|
/* 4. Approximate 2^x2 - 1, using a fourth-degree polynomial,
|
|
with maximum error in [-2^-10-2^-30,2^-10+2^-30]
|
|
less than 10^-19. */
|
|
|
|
x22 = (((.0096181293647031180
|
|
* x + .055504110254308625)
|
|
* x + .240226506959100583)
|
|
* x + .69314718055994495) * ex2_u.d;
|
|
math_opt_barrier (x22);
|
|
}
|
|
|
|
/* 5. Return (2^x2-1) * 2^(t/512+e+ex) + 2^(t/512+e+ex). */
|
|
result = x22 * x + ex2_u.d;
|
|
|
|
if (!unsafe)
|
|
return result;
|
|
else
|
|
return result * scale_u.d;
|
|
}
|
|
else
|
|
/* Return x, if x is a NaN or Inf; or overflow, otherwise. */
|
|
return TWO1023 * x;
|
|
}
|
|
strong_alias (__ieee754_exp2, __exp2_finite)
|