mirror of
https://sourceware.org/git/glibc.git
synced 2024-11-22 21:10:07 +00:00
903af5af9a
Various exp2 implementations in glibc can miss underflow exceptions when the scaling down part of the calculation is exact (or, in the x86 case, when the conversion from extended precision to the target precision is exact). This patch forces the exception in a similar way to previous fixes. The x86 exp2f changes may in fact not be needed for this purpose - it's likely to be the case that no argument of type float has an exp2 result so close to an exact subnormal float value that it equals that value when rounded to 64 bits (even taking account of variation between different x86 implementations). However, they are included for consistency with the changes to exp2 and so as to fix the exp2f part of bug 18875 by ensuring that excess range and precision is removed from underflowing return values. Tested for x86_64, x86 and mips64. [BZ #16521] [BZ #18875] * math/e_exp2l.c (__ieee754_exp2l): Force underflow exception for small results. * sysdeps/i386/fpu/e_exp2.S (dbl_min): New object. (MO): New macro. (__ieee754_exp2): For small results, force underflow exception and remove excess range and precision from return value. * sysdeps/i386/fpu/e_exp2f.S (flt_min): New object. (MO): New macro. (__ieee754_exp2f): For small results, force underflow exception and remove excess range and precision from return value. * sysdeps/i386/fpu/e_exp2l.S (ldbl_min): New object. (MO): New macro. (__ieee754_exp2l): Force underflow exception for small results. * 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 (ldbl_min): New object. (MO): New macro. (__ieee754_exp2l): Force underflow exception for small results. * math/auto-libm-test-in: Add more tests or exp2. * math/auto-libm-test-out: Regenerated.
67 lines
1.9 KiB
C
67 lines
1.9 KiB
C
/* Compute 2^x.
|
|
Copyright (C) 2012-2015 Free Software Foundation, Inc.
|
|
This file is part of the GNU C Library.
|
|
|
|
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/>. */
|
|
|
|
#include <math.h>
|
|
#include <math_private.h>
|
|
#include <float.h>
|
|
|
|
/* To avoid spurious underflows, use this definition to treat IBM long
|
|
double as approximating an IEEE-style format. */
|
|
#if LDBL_MANT_DIG == 106
|
|
# undef LDBL_EPSILON
|
|
# define LDBL_EPSILON 0x1p-106L
|
|
#endif
|
|
|
|
long double
|
|
__ieee754_exp2l (long double x)
|
|
{
|
|
if (__glibc_likely (isless (x, (long double) LDBL_MAX_EXP)))
|
|
{
|
|
if (__builtin_expect (isgreaterequal (x, (long double) (LDBL_MIN_EXP
|
|
- LDBL_MANT_DIG
|
|
- 1)), 1))
|
|
{
|
|
int intx = (int) x;
|
|
long double fractx = x - intx;
|
|
long double result;
|
|
if (fabsl (fractx) < LDBL_EPSILON / 4.0L)
|
|
result = __scalbnl (1.0L + fractx, intx);
|
|
else
|
|
result = __scalbnl (__ieee754_expl (M_LN2l * fractx), intx);
|
|
if (result < LDBL_MIN)
|
|
{
|
|
long double force_underflow = result * result;
|
|
math_force_eval (force_underflow);
|
|
}
|
|
return result;
|
|
}
|
|
else
|
|
{
|
|
/* Underflow or exact zero. */
|
|
if (isinf (x))
|
|
return 0;
|
|
else
|
|
return LDBL_MIN * LDBL_MIN;
|
|
}
|
|
}
|
|
else
|
|
/* Infinity, NaN or overflow. */
|
|
return LDBL_MAX * x;
|
|
}
|
|
strong_alias (__ieee754_exp2l, __exp2l_finite)
|