C99 and C11 allow but do not require ceil, floor, round and trunc to
raise the "inexact" exception for noninteger arguments. TS 18661-1
requires that this exception not be raised by these functions. This
aligns them with general IEEE semantics, where "inexact" is only
raised if the final step of rounding the infinite-precision result to
the result type is inexact; for these functions, the
infinite-precision integer result is always representable in the
result type, so "inexact" should never be raised.
The generic implementations of ceil, floor and round functions contain
code to force "inexact" to be raised. This patch removes it for floor
functions to align them with TS 18661-1 in this regard. Note that
some architecture-specific versions may still raise "inexact", so the
tests are not updated and the bug is not yet fixed.
Tested for x86_64, x86 and mips64.
[BZ #15479]
* sysdeps/ieee754/dbl-64/s_floor.c: Do not mention "inexact"
exception in comment.
(huge): Remove variable.
(__floor): Do not force "inexact" exception.
* sysdeps/ieee754/dbl-64/wordsize-64/s_floor.c: Do not mention
"inexact" exception in comment.
(huge): Remove variable.
(__floor): Do not force "inexact" exception.
* sysdeps/ieee754/flt-32/s_floorf.c: Do not mention "inexact"
exception in comment.
(huge): Remove variable.
(__floorf): Do not force "inexact" exception.
* sysdeps/ieee754/ldbl-128/s_floorl.c: Do not mention "inexact"
exception in comment.
(huge): Remove variable.
(__floorl): Do not force "inexact" exception.
C99 and C11 allow but do not require ceil, floor, round and trunc to
raise the "inexact" exception for noninteger arguments. TS 18661-1
requires that this exception not be raised by these functions. This
aligns them with general IEEE semantics, where "inexact" is only
raised if the final step of rounding the infinite-precision result to
the result type is inexact; for these functions, the
infinite-precision integer result is always representable in the
result type, so "inexact" should never be raised.
The generic implementations of ceil, floor and round functions contain
code to force "inexact" to be raised. This patch removes it for ceil
functions to align them with TS 18661-1 in this regard. Note that
some architecture-specific versions may still raise "inexact", so the
tests are not updated and the bug is not yet fixed.
Tested for x86_64, x86 and mips64.
[BZ #15479]
* sysdeps/ieee754/dbl-64/s_ceil.c: Do not mention "inexact"
exception in comment.
(huge): Remove variable.
(__ceil): Do not force "inexact" exception.
* sysdeps/ieee754/dbl-64/wordsize-64/s_ceil.c: Do not mention
"inexact" exception in comment.
(huge): Remove variable.
(__ceil): Do not force "inexact" exception.
* sysdeps/ieee754/flt-32/s_ceilf.c (huge): Remove variable.
(__ceilf): Do not force "inexact" exception.
* sysdeps/ieee754/ldbl-128/s_ceill.c: Do not mention "inexact"
exception in comment.
(huge): Remove variable.
(__ceill): Do not force "inexact" exception.
The changes to restrict implementation-namespace symbol aliases such
as __finitel to compat symbols used code for __finitel in libm
analogous to that for __finitel in libc. However, the versions for
the two symbols are actually different, GLIBC_2.0 in libc and
GLIBC_2.1 in libm. This patch fixes the handling of the libm compat
symbol.
Tested for mips (o32), where it fixes an ABI test failure.
* sysdeps/ieee754/dbl-64/s_finite.c
[NO_LONG_DOUBLE && LDBL_CLASSIFY_COMPAT] (__finitel): Define
compat symbol at version GLIBC_2_1 and use GLIBC_2_1 in
SHLIB_COMPAT condition for libm, not GLIBC_2_0.
* sysdeps/ieee754/dbl-64/wordsize-64/s_finite.c
[NO_LONG_DOUBLE && LDBL_CLASSIFY_COMPAT] (__finitel): Likewise.
When looking at the code generated for pow() on ppc64 I noticed quite
a few sign extensions. Making the array indices unsigned reduces the
number of sign extensions from 24 to 7.
Tested for powerpc64le and x86_64.
Like the previous change, exploit the fact that computation for sin
and cos is identical except that it is apart by a quadrant. Also
remove csloww, csloww1 and csloww2 since they can easily be expressed
in terms of sloww, sloww1 and sloww2.
The sin and cos computation for this range of input is identical
except for a difference in quadrants by 1. Exploit that fact and the
common argument reduction to reduce computations for sincos.
Range reduction needs to be done only once for sin and cos, so copy
over all of the relevant functions (__sin, __cos, reduce_and_compute)
and consolidate common code.
If a platform does not define "long-double-fcts = yes" in its
Makefiles and it does define __NO_LONG_DOUBLE_MATH in its installed
headers, it will currently create exported symbols for __finitel,
__isinfl, and __isnanl that can't be reached from userspace by
correct use of the finite(), isinf(), or isnan() macros in <math.h>.
To avoid this situation, by default for such platforms we now no
longer export these symbols, thus causing appropriate link-time
errors. However, for platforms that previously exported these
symbols, we continue to do so as compat symbols; this is enabled
by adding LDBL_CLASSIFY_COMPAT to math_private.h for the platform.
For tile, remove the now-unnecessary exports of those functions from
libc and libm.
Various sysdeps/ieee754/dbl-64 functions use double constants defined
using a union between a double and two ints, with separate big-endian
and little-endian definitions of the constants.
With modern C, this is unnecessary complication; hex float constants
(or __builtin_inf etc.) suffice to specify the exact value desired,
and so can avoid separate versions for each endianness. Having this
complication also complicates cleanups such as removing slow paths
from these library functions, as they need to make sure to remove both
copies of variables that are no longer used after such a cleanup (and
in at least one case, proper removal of a slow path will also involve
removing slow-path-only values from the middle of an array - an array
with both big-endian and little-endian copies - and adjusting other
references to that array).
So it makes sense to clean up the code to define these constants using
hex floats and so eliminate the endianness conditional. This patch
does so in the case of sqrt, where the two constants are such that it
makes sense just to put them directly in the code using them and
eliminate the names for them altogether.
Tested for arm (the code generated for sqrt does change, though not in
any significant way).
* sysdeps/ieee754/dbl-64/e_sqrt.c: Do not include uroot.h.
(__ieee754_sqrt): Use hex float constants instead of tm256.x and
t512.x.
* sysdeps/ieee754/dbl-64/uroot.h: Remove file.
Include the __sin and __cos functions as local static copies to allow
deper optimization of the functions. This change shows an improvement
of about 17% in the min case and 12.5% in the mean case for the sincos
microbenchmark on x86_64.
* sysdeps/ieee754/dbl-64/s_sin.c (__sin)[IN_SINCOS]: Mark function
static and don't set or restore rounding.
(__cos)[IN_SINCOS]: Likewise.
* sysdeps/ieee754/dbl-64/s_sincos.c: Include s_sin.c.
(__sincos): Set and restore rounding mode. Remove check for infinite
or NaN input.
For some large arguments, the dbl-64 implementation of remainder gives
zero results with the wrong sign, resulting from a subtraction that is
mathematically correct but does not guarantee that a zero result has
the sign of the first argument to remainder. This patch adds an
appropriate check for this case, similar to other implementations of
remainder in the case of equality, and adds tests of remainder on
inputs already used to test remquo.
Tested for x86_64 and x86.
[BZ #19201]
* sysdeps/ieee754/dbl-64/e_remainder.c (__ieee754_remainder):
Check for zero remainder in case of large exponents and ensure
correct sign of result in that case.
* math/libm-test.inc (remainder_test_data): Add more tests.
C11 defines standard <float.h> macros *_TRUE_MIN for the least
positive subnormal value of a type. Now that we build with
-std=gnu11, we can use these macros in glibc. This patch replaces
previous uses of the GCC predefines __*_DENORM_MIN__ (used in
<float.h> to define *_TRUE_MIN), as well as *_DENORM_MIN references in
comments.
Tested for x86_64 and x86 (testsuite, and that installed shared
libraries are unchanged by the patch). Also tested for powerpc that
installed stripped shared libraries are unchanged by the patch.
* math/libm-test.inc (min_subnorm_value): Use LDBL_TRUE_MIN,
DBL_TRUE_MIN and FLT_TRUE_MIN instead of __LDBL_DENORM_MIN__,
__DBL_DENORM_MIN__ and __FLT_DENORM_MIN__.
* sysdeps/ieee754/dbl-64/s_fma.c (__fma): Refer to DBL_TRUE_MIN
instead of DBL_DENORM_MIN in comment.
* sysdeps/ieee754/ldbl-128/s_fmal.c (__fmal): Refer to
LDBL_TRUE_MIN instead of LDBL_DENORM_MIN in comment.
* sysdeps/ieee754/ldbl-128ibm/s_nextafterl.c: Include <float.h>.
(__nextafterl): Use LDBL_TRUE_MIN instead of __LDBL_DENORM_MIN__.
* sysdeps/ieee754/ldbl-96/s_fmal.c (__fmal): Refer to
LDBL_TRUE_MIN instead of LDBL_DENORM_MIN in comment.
One common case of __GNUC_PREREQ (4, 7) conditionals is use of
diagnostic control pragmas for -Wmaybe-uninitialized, an option
introduced in GCC 4.7 where older GCC needed -Wuninitialized to be
controlled instead if the warning appeared with older GCC. This patch
removes such conditionals.
(There remain several older uses of -Wno-uninitialized in makefiles
that still need to be converted to diagnostic control pragmas if the
issue is still present with current sources and supported GCC
versions, and it's likely that in most cases those pragmas also will
end up controlling -Wmaybe-uninitialized.)
Tested for x86_64 and x86 (testsuite, and that installed stripped
shared libraries are unchanged by the patch, except for libresolv
since res_send.c contains assertions whose line numbers are changed by
the patch).
* resolv/res_send.c (send_vc) [__GNUC_PREREQ (4, 7)]: Make code
unconditional.
* soft-fp/fmadf4.c [__GNUC_PREREQ (4, 7)]: Likewise.
[!__GNUC_PREREQ (4, 7)]: Remove conditional code.
* soft-fp/fmasf4.c [__GNUC_PREREQ (4, 7)]: Make code
unconditional.
[!__GNUC_PREREQ (4, 7)]: Remove conditional code.
* soft-fp/fmatf4.c [__GNUC_PREREQ (4, 7)]: Make code
unconditional.
[!__GNUC_PREREQ (4, 7)]: Remove conditional code.
* stdlib/setenv.c
[((__GNUC__ << 16) + __GNUC_MINOR__) >= ((4 << 16) + 7)]: Make
code unconditional.
[!(((__GNUC__ << 16) + __GNUC_MINOR__) >= ((4 << 16) + 7))]:
Remove conditional code.
* sysdeps/ieee754/dbl-64/e_lgamma_r.c
(__ieee754_lgamma_r) [__GNUC_PREREQ (4, 7)]: Make code
unconditional.
(__ieee754_lgamma_r) [!__GNUC_PREREQ (4, 7)]: Remove conditional
code.
* sysdeps/ieee754/flt-32/e_lgammaf_r.c
(__ieee754_lgammaf_r) [__GNUC_PREREQ (4, 7)]: Make code
unconditional.
(__ieee754_lgammaf_r) [!__GNUC_PREREQ (4, 7)]: Remove conditional
code.
* sysdeps/ieee754/ldbl-128/k_tanl.c
(__kernel_tanl) [__GNUC_PREREQ (4, 7)]: Make code unconditional.
(__kernel_tanl) [!__GNUC_PREREQ (4, 7)]: Remove conditional code.
* sysdeps/ieee754/ldbl-128ibm/k_tanl.c
(__kernel_tanl) [__GNUC_PREREQ (4, 7)]: Make code unconditional.
(__kernel_tanl) [!__GNUC_PREREQ (4, 7)]: Remove conditional code.
* sysdeps/ieee754/ldbl-96/e_lgammal_r.c
(__ieee754_lgammal_r) [__GNUC_PREREQ (4, 7)]: Make code
unconditional.
(__ieee754_lgammal_r) [!__GNUC_PREREQ (4, 7)]: Remove conditional
code.
* sysdeps/ieee754/ldbl-96/k_tanl.c
(__kernel_tanl) [__GNUC_PREREQ (4, 7)]: Make code unconditional.
(__kernel_tanl) [!__GNUC_PREREQ (4, 7)]: Remove conditional code.
j1 and jn can underflow for small arguments, but fail to set errno
when underflowing to 0. This patch fixes them to set errno in that
case.
Tested for x86_64, x86, mips64 and powerpc.
[BZ #18611]
* sysdeps/ieee754/dbl-64/e_j1.c (__ieee754_j1): Set errno and
avoid excess range and precision on underflow.
* sysdeps/ieee754/dbl-64/e_jn.c (__ieee754_jn): Likewise.
* sysdeps/ieee754/flt-32/e_j1f.c (__ieee754_j1f): Likewise.
* sysdeps/ieee754/flt-32/e_jnf.c (__ieee754_jnf): Likewise.
* sysdeps/ieee754/ldbl-128/e_j1l.c (__ieee754_j1l): Set errno on
underflow.
* sysdeps/ieee754/ldbl-128/e_jnl.c (__ieee754_jnl): Likewise.
* sysdeps/ieee754/ldbl-128ibm/e_jnl.c (__ieee754_jnl): Likewise.
* sysdeps/ieee754/ldbl-96/e_j1l.c (__ieee754_j1l): Likewise.
* sysdeps/ieee754/ldbl-96/e_jnl.c (__ieee754_jnl): Likewise.
* math/auto-libm-test-in: Do not allow missing errno setting for
tests of j1 and jn.
* math/auto-libm-test-out: Regenerated.
For 32-bit MIPS and some other systems, various of the lrint, llrint,
lround, llround functions can be missing exceptions on overflow
because casts do not (in current GCC) result in the proper
exceptions. In the MIPS case there are two problems here: MIPS I code
generation uses an assembler macro that doesn't raise exceptions,
while the libgcc conversions of floating-point values to long long
also do not raise "invalid" on all overflow cases (and can raise
spurious "inexact").
This patch adds support in the generic code (only the functions for
which this problem has actually been seen) for forcing the "invalid"
exception in the problem cases, and enables that support for the
affected MIPS cases.
Tested for MIPS; also tested for x86_64 and x86 that installed
stripped shared libraries are unchanged by this patch.
[BZ #16399]
* sysdeps/generic/fix-fp-int-convert-overflow.h: New file.
* sysdeps/ieee754/dbl-64/s_llrint.c: Include <fenv.h>, <limits.h>
and <fix-fp-int-convert-overflow.h>.
(__llrint) [FE_INVALID]: Force FE_INVALID exception as needed if
FIX_DBL_LLONG_CONVERT_OVERFLOW.
* sysdeps/ieee754/dbl-64/s_llround.c: Include <fenv.h>, <limits.h>
and <fix-fp-int-convert-overflow.h>.
(__llround) [FE_INVALID]: Force FE_INVALID exception as needed if
FIX_DBL_LLONG_CONVERT_OVERFLOW.
* sysdeps/ieee754/dbl-64/s_lrint.c: Include
<fix-fp-int-convert-overflow.h>.
(__lrint) [FE_INVALID]: Force FE_INVALID exception as needed if
FIX_DBL_LLONG_CONVERT_OVERFLOW.
* sysdeps/ieee754/dbl-64/s_lround.c: Include
<fix-fp-int-convert-overflow.h>.
(__lround) [FE_INVALID]: Force FE_INVALID exception as needed if
FIX_DBL_LLONG_CONVERT_OVERFLOW.
* sysdeps/ieee754/flt-32/s_llrintf.c: Include <fenv.h>, <limits.h>
and <fix-fp-int-convert-overflow.h>.
(__llrintf) [FE_INVALID]: Force FE_INVALID exception as needed if
FIX_DBL_LLONG_CONVERT_OVERFLOW.
* sysdeps/ieee754/flt-32/s_llroundf.c: Include <fenv.h>,
<limits.h> and <fix-fp-int-convert-overflow.h>.
(__llroundf) [FE_INVALID]: Force FE_INVALID exception as needed if
FIX_DBL_LLONG_CONVERT_OVERFLOW.
* sysdeps/ieee754/flt-32/s_lrintf.c: Include <fenv.h>, <limits.h>
and <fix-fp-int-convert-overflow.h>.
(__lrintf) [FE_INVALID]: Force FE_INVALID exception as needed if
FIX_DBL_LLONG_CONVERT_OVERFLOW.
* sysdeps/ieee754/flt-32/s_lroundf.c: Include <fenv.h>, <limits.h>
and <fix-fp-int-convert-overflow.h>.
(__lroundf) [FE_INVALID]: Force FE_INVALID exception as needed if
FIX_DBL_LLONG_CONVERT_OVERFLOW.
* sysdeps/mips/mips32/fpu/fix-fp-int-convert-overflow.h: New file.
The dbl-64 implementation of lrint produces incorrect results for some
arguments with 64-bit long because a 32-bit (unsigned) low part of the
mantissa is shifted left, losing high bits in the process. This patch
fixes this by casting to long int before shifting, as in lround (as
this case only applies for 64-bit long, there are no issues with
sign-extension).
Tested for mips64 (n64).
[BZ #19095]
* sysdeps/ieee754/dbl-64/s_lrint.c (__lrint): Cast low part of
mantissa to long int before shifting left.
The dbl-64, ldbl-96 and ldbl-128 implementations of lrint and llrint
fail to produce "invalid" exceptions in cases where the rounded result
overflows the target type, but truncating the floating-point argument
to the next integer towards zero does not overflow it (so in
particular casts do not produce such exceptions). (This issue cannot
arise for float, or for double with 64-bit target type, or for ldbl-96
with 64-bit target type and negative arguments, because of
insufficient precision in the floating-point type for arguments with
the relevant property to exist. It also obviously cannot arise in
FE_TOWARDZERO mode.)
This patch fixes these problems by inserting checks for the special
cases that can occur in each implementation, and explicitly raising
FE_INVALID (and avoiding the cast if it might raise spurious
FE_INEXACT, while raising FE_INEXACT explicitly in the cases where it
is needed; unlike lround and llround, FE_INEXACT is required, not
optional, for these functions for a within-range inexact result).
The fixes are conditional on FE_INVALID or FE_INEXACT being defined.
If any future architecture supports one but not both of those
exceptions, the code will fail to compile and need fixing to handle
that case (this seemed better than conditioning on both macros being
defined, resulting in code that would compile but quietly miss
exceptions on such a system).
Tested for x86_64, x86 and mips64. Tested the ldbl-96 changes (only
relevant for ia64, it appears) on x86_64 by removing the x86_64
versions of lrintl / llrintl.
[BZ #19094]
* sysdeps/ieee754/dbl-64/s_lrint.c: Include <fenv.h> and
<limits.h>.
(__lrint) [FE_INVALID || FE_INEXACT]: Force FE_INVALID exception
when result overflows but exception would not result from cast.
* sysdeps/ieee754/ldbl-128/s_llrintl.c: Include <fenv.h> and
<limits.h>.
(__llrintl) [FE_INVALID || FE_INEXACT]: Force FE_INVALID exception
when result overflows but exception would not result from cast.
* sysdeps/ieee754/ldbl-128/s_lrintl.c: Include <fenv.h> and
<limits.h>.
(__lrintl) [FE_INVALID || FE_INEXACT]: Force FE_INVALID exception
when result overflows but exception would not result from cast.
* sysdeps/ieee754/ldbl-96/s_llrintl.c: Include <fenv.h> and
<limits.h>.
(__llrintl) [FE_INVALID || FE_INEXACT]: Force FE_INVALID exception
when result overflows but exception would not result from cast.
* sysdeps/ieee754/ldbl-96/s_lrintl.c: Include <fenv.h> and
<limits.h>.
(__lrintl) [FE_INVALID || FE_INEXACT]: Force FE_INVALID exception
when result overflows but exception would not result from cast.
* math/libm-test.inc (lrint_test_data): Add more tests.
(llrint_test_data): Likewise.
The dbl-64, ldbl-96 and ldbl-128 implementations of lround and llround
fail to produce "invalid" exceptions in cases where the rounded result
overflows the target type, but truncating the floating-point argument
to the next integer towards zero does not overflow it (so in
particular casts do not produce such exceptions). (This issue cannot
arise for float, or for double with 64-bit target type, or for ldbl-96
with 64-bit target type and negative arguments, because of
insufficient precision in the floating-point type for arguments with
the relevant property to exist.)
This patch fixes these problems by inserting checks for the special
cases that can occur in each implementation, and explicitly raising
FE_INVALID (and avoiding the cast if it might raise spurious
FE_INEXACT).
Tested for x86_64, x86 and mips64.
[BZ #19088]
* sysdeps/ieee754/dbl-64/s_lround.c: Include <fenv.h> and
<limits.h>.
(__lround) [FE_INVALID]: Force FE_INVALID exception when result
overflows but exception would not result from cast.
* sysdeps/ieee754/dbl-64/wordsize-64/s_lround.c: Include <fenv.h>
and <limits.h>.
(__lround) [FE_INVALID]: Force FE_INVALID exception when result
overflows but exception would not result from cast.
* sysdeps/ieee754/ldbl-128/s_llroundl.c: Include <fenv.h> and
<limits.h>.
(__llroundl) [FE_INVALID]: Force FE_INVALID exception when result
overflows but exception would not result from cast.
* sysdeps/ieee754/ldbl-128/s_lroundl.c: Include <fenv.h> and
<limits.h>.
(__lroundl) [FE_INVALID]: Force FE_INVALID exception when result
overflows but exception would not result from cast.
* sysdeps/ieee754/ldbl-96/s_llroundl.c: Include <fenv.h> and
<limits.h>.
(__llroundl) [FE_INVALID]: Force FE_INVALID exception when result
overflows but exception would not result from cast.
* sysdeps/ieee754/ldbl-96/s_lroundl.c: Include <fenv.h> and
<limits.h>.
(__lroundl) [FE_INVALID]: Force FE_INVALID exception when result
overflows but exception would not result from cast.
* math/libm-test.inc (lround_test_data): Add more tests.
(llround_test_data): Likewise.
This patch enables use of sysdeps/ieee754/dbl-64/wordsize-64 for
MIPS64 (both n64 and n32), removing a #error in one case now that case
has been tested and found to work.
Tested for mips64 (n64 and n32).
* sysdeps/mips/mips64/Implies: Use ieee754/dbl-64/wordsize-64.
* sysdeps/ieee754/dbl-64/wordsize-64/s_issignaling.c
(__issignaling) [HIGH_ORDER_BIT_IS_SET_FOR_SNAN]: Remove #error.
The implementation of lround in dbl-64/wordsize-64 as an alias or
wrapper for llround is always incorrect when long is not 64-bit,
because it misses required exceptions in overflow cases, as shown by
my recently added tests. This patch removes that alias / wrapper in
the non-LP64 case, together with the REGISTER_CAST_INT32_TO_INT64
macro, restoring the previous version of lround for dbl-64/wordsize-64
(newly conditioned on !_LP64).
Tested for x86_64, and for mips64 with use of dbl-64/wordsize-64
enabled.
[BZ #19079]
* sysdeps/ieee754/dbl-64/wordsize-64/s_lround.c: Restore previous
file, conditioned on [!_LP64].
* sysdeps/ieee754/dbl-64/wordsize-64/s_llround.c
[!_LP64] (__lround): Do not define as function or alias.
[!_LP64] (lround): Likewise.
[!_LP64] (__lroundl): Likewise.
[!_LP64] (lroundl): Likewise.
* sysdeps/tile/sysdep.h (REGISTER_CAST_INT32_TO_INT64): Remove
macro.
* sysdeps/x86_64/x32/sysdep.h (REGISTER_CAST_INT32_TO_INT64):
Likewise.
On powerpc32 hard-float, older processors (ones where fcfid is not
available for 32-bit code), GCC generates conversions from integers to
floating point that wrongly convert integer 0 to -0 instead of +0 in
FE_DOWNWARD mode. This in turn results in logb and a few other
functions wrongly returning -0 when they should return +0.
This patch works around this issue in glibc as I proposed in
<https://sourceware.org/ml/libc-alpha/2015-09/msg00728.html>, so that
the affected functions can be correct and the affected tests pass in
the absence of a GCC fix for this longstanding issue (GCC bug 67771 -
if fixed, of course we can put in GCC version conditionals, and
eventually phase out the workarounds). A new macro
FIX_INT_FP_CONVERT_ZERO is added in a new sysdeps header
fix-int-fp-convert-zero.h, and the powerpc32/fpu version of that
header defines the macro based on the results of a configure test for
whether such conversions use the fcfid instruction.
Tested for x86_64 (that installed stripped shared libraries are
unchanged by the patch) and powerpc (that HAVE_PPC_FCFID comes out to
0 as expected and that the relevant tests are fixed). Also tested a
build with GCC configured for -mcpu=power4 and verified that
HAVE_PPC_FCFID comes out to 1 in that case.
There are still some other issues to fix to get test-float and
test-double passing cleanly for older powerpc32 processors (apart from
the need for an ulps regeneration for powerpc). (test-ldouble will be
harder to get passing cleanly, but with a combination of selected
fixes to ldbl-128ibm code that don't involve significant performance
issues, allowing spurious underflow and inexact exceptions for that
format, and lots of XFAILing for the default case of unpatched libgcc,
it should be doable.)
[BZ #887]
[BZ #19049]
[BZ #19050]
* sysdeps/generic/fix-int-fp-convert-zero.h: New file.
* sysdeps/ieee754/dbl-64/e_log10.c: Include
<fix-int-fp-convert-zero.h>.
(__ieee754_log10): Adjust signs as needed if FIX_INT_FP_CONVERT_ZERO.
* sysdeps/ieee754/dbl-64/e_log2.c: Include
<fix-int-fp-convert-zero.h>.
(__ieee754_log2): Adjust signs as needed if FIX_INT_FP_CONVERT_ZERO.
* sysdeps/ieee754/dbl-64/s_erf.c: Include
<fix-int-fp-convert-zero.h>.
(__erfc): Adjust signs as needed if FIX_INT_FP_CONVERT_ZERO.
* sysdeps/ieee754/dbl-64/s_logb.c: Include
<fix-int-fp-convert-zero.h>.
(__logb): Adjust signs as needed if FIX_INT_FP_CONVERT_ZERO.
* sysdeps/ieee754/flt-32/e_log10f.c: Include
<fix-int-fp-convert-zero.h>.
(__ieee754_log10f): Adjust signs as needed if FIX_INT_FP_CONVERT_ZERO.
* sysdeps/ieee754/flt-32/e_log2f.c: Include
<fix-int-fp-convert-zero.h>.
(__ieee754_log2f): Adjust signs as needed if FIX_INT_FP_CONVERT_ZERO.
* sysdeps/ieee754/flt-32/s_erff.c: Include
<fix-int-fp-convert-zero.h>.
(__erfcf): Adjust signs as needed if FIX_INT_FP_CONVERT_ZERO.
* sysdeps/ieee754/flt-32/s_logbf.c: Include
<fix-int-fp-convert-zero.h>.
(__logbf): Adjust signs as needed if FIX_INT_FP_CONVERT_ZERO.
* sysdeps/ieee754/ldbl-128ibm/s_erfl.c: Include
<fix-int-fp-convert-zero.h>.
(__erfcl): Adjust signs as needed if FIX_INT_FP_CONVERT_ZERO.
* sysdeps/ieee754/ldbl-128ibm/s_logbl.c: Include
<fix-int-fp-convert-zero.h>.
(__logbl): Adjust signs as needed if FIX_INT_FP_CONVERT_ZERO.
* sysdeps/powerpc/powerpc32/fpu/configure.ac: New file.
* sysdeps/powerpc/powerpc32/fpu/configure: New generated file.
* sysdeps/powerpc/powerpc32/fpu/fix-int-fp-convert-zero.h: New
file.
* config.h.in [_LIBC] (HAVE_PPC_FCFID): New macro.
For arguments with X^2 + Y^2 close to 1, clog and clog10 avoid large
errors from log(hypot) by computing X^2 + Y^2 - 1 in a way that avoids
cancellation error and then using log1p.
However, the thresholds for using that approach still result in log
being used on argument as large as sqrt(13/16) > 0.9, leading to
significant errors, in some cases above the 9ulp maximum allowed in
glibc libm. This patch arranges for the approach using log1p to be
used in any cases where |X|, |Y| < 1 and X^2 + Y^2 >= 0.5 (with the
existing allowance for cases where one of X and Y is very small),
adjusting the __x2y2m1 functions to work with the wider range of
inputs. This way, log only gets used on arguments below sqrt(1/2) (or
substantially above 1), where the error involved is much less.
Tested for x86_64, x86, mips64 and powerpc. For the ulps regeneration
I removed the existing clog and clog10 ulps before regenerating to
allow any reduced ulps to appear. Tests added include those found by
random test generation to produce large ulps either before or after
the patch, and some found by trying inputs close to the (0.75, 0.5)
threshold where the potential errors from using log are largest.
[BZ #19016]
* sysdeps/generic/math_private.h (__x2y2m1f): Update comment to
allow more cases with X^2 + Y^2 >= 0.5.
* sysdeps/ieee754/dbl-64/x2y2m1.c (__x2y2m1): Likewise. Add -1 as
normal element in sum instead of special-casing based on values of
arguments.
* sysdeps/ieee754/dbl-64/x2y2m1f.c (__x2y2m1f): Update comment.
* sysdeps/ieee754/ldbl-128/x2y2m1l.c (__x2y2m1l): Likewise. Add
-1 as normal element in sum instead of special-casing based on
values of arguments.
* sysdeps/ieee754/ldbl-128ibm/x2y2m1l.c (__x2y2m1l): Likewise.
* sysdeps/ieee754/ldbl-96/x2y2m1.c [FLT_EVAL_METHOD != 0]
(__x2y2m1): Update comment.
* sysdeps/ieee754/ldbl-96/x2y2m1l.c (__x2y2m1l): Likewise. Add -1
as normal element in sum instead of special-casing based on values
of arguments.
* math/s_clog.c (__clog): Handle more cases using log1p without
hypot.
* math/s_clog10.c (__clog10): Likewise.
* math/s_clog10f.c (__clog10f): Likewise.
* math/s_clog10l.c (__clog10l): Likewise.
* math/s_clogf.c (__clogf): Likewise.
* math/s_clogl.c (__clogl): Likewise.
* math/auto-libm-test-in: Add more tests of clog and clog10.
* math/auto-libm-test-out: Regenerated.
* sysdeps/i386/fpu/libm-test-ulps: Update.
* sysdeps/x86_64/fpu/libm-test-ulps: Likewise.
Similar to various other bugs in this area, pow functions can fail to
raise the underflow exception when the result is tiny and inexact but
one or more low bits of the intermediate result that is scaled down
(or, in the i386 case, converted from a wider evaluation format) are
zero. This patch forces the exception in a similar way to previous
fixes, thereby concluding the fixes for known bugs with missing
underflow exceptions currently filed in Bugzilla.
Tested for x86_64, x86, mips64 and powerpc.
[BZ #18825]
* sysdeps/i386/fpu/i386-math-asm.h (FLT_NARROW_EVAL_UFLOW_NONNAN):
New macro.
(DBL_NARROW_EVAL_UFLOW_NONNAN): Likewise.
(LDBL_CHECK_FORCE_UFLOW_NONNAN): Likewise.
* sysdeps/i386/fpu/e_pow.S: Use DEFINE_DBL_MIN.
(__ieee754_pow): Use DBL_NARROW_EVAL_UFLOW_NONNAN instead of
DBL_NARROW_EVAL, reloading the PIC register as needed.
* sysdeps/i386/fpu/e_powf.S: Use DEFINE_FLT_MIN.
(__ieee754_powf): Use FLT_NARROW_EVAL_UFLOW_NONNAN instead of
FLT_NARROW_EVAL. Use separate return path for case when first
argument is NaN.
* sysdeps/i386/fpu/e_powl.S: Include <i386-math-asm.h>. Use
DEFINE_LDBL_MIN.
(__ieee754_powl): Use LDBL_CHECK_FORCE_UFLOW_NONNAN, reloading the
PIC register.
* sysdeps/ieee754/dbl-64/e_pow.c (__ieee754_pow): Use
math_check_force_underflow_nonneg.
* sysdeps/ieee754/flt-32/e_powf.c (__ieee754_powf): Force
underflow for subnormal result.
* sysdeps/ieee754/ldbl-128/e_powl.c (__ieee754_powl): Likewise.
* sysdeps/ieee754/ldbl-128ibm/e_powl.c (__ieee754_powl): Use
math_check_force_underflow_nonneg.
* sysdeps/x86/fpu/powl_helper.c (__powl_helper): Use
math_check_force_underflow.
* sysdeps/x86_64/fpu/x86_64-math-asm.h
(LDBL_CHECK_FORCE_UFLOW_NONNAN): New macro.
* sysdeps/x86_64/fpu/e_powl.S: Include <x86_64-math-asm.h>. Use
DEFINE_LDBL_MIN.
(__ieee754_powl): Use LDBL_CHECK_FORCE_UFLOW_NONNAN.
* math/auto-libm-test-in: Add more tests of pow.
* math/auto-libm-test-out: Regenerated.
Similar to various other bugs in this area, hypot functions can fail
to raise the underflow exception when the result is tiny and inexact
but one or more low bits of the intermediate result that is scaled
down (or, in the i386 case, converted from a wider evaluation format)
are zero. This patch forces the exception in a similar way to
previous fixes.
Note that this issue cannot arise for implementations of hypotf using
double (or wider) for intermediate evaluation (if hypotf should
underflow, that means the double square root is being computed of some
number of the form N*2^-298, for 0 < N < 2^46, which is exactly
represented as a double, and whatever the rounding mode such a square
root cannot have a mantissa with all zeroes after the initial 23
bits). Thus no changes are made to hypotf implementations in this
patch, only to hypot and hypotl.
Tested for x86_64, x86, mips64 and powerpc.
[BZ #18803]
* sysdeps/i386/fpu/e_hypot.S: Use DEFINE_DBL_MIN.
(MO): New macro.
(__ieee754_hypot) [PIC]: Load PIC register.
(__ieee754_hypot): Use DBL_NARROW_EVAL_UFLOW_NONNEG instead of
DBL_NARROW_EVAL.
* sysdeps/ieee754/dbl-64/e_hypot.c (__ieee754_hypot): Use
math_check_force_underflow_nonneg in case where result might be
tiny.
* sysdeps/ieee754/ldbl-128/e_hypotl.c (__ieee754_hypotl):
Likewise.
* sysdeps/ieee754/ldbl-128ibm/e_hypotl.c (__ieee754_hypotl):
Likewise.
* sysdeps/ieee754/ldbl-96/e_hypotl.c (__ieee754_hypotl): Likewise.
* sysdeps/powerpc/fpu/e_hypot.c (__ieee754_hypot): Likewise.
* math/auto-libm-test-in: Add more tests of hypot.
* math/auto-libm-test-out: Regenerated.
Where glibc code needs to avoid excess range and precision in
floating-point arithmetic, code variously uses either asms or volatile
to force the results of that arithmetic to memory; mostly this is
conditional on FLT_EVAL_METHOD, but in the case of lrint / llrint
functions some use of volatile is unconditional (and is present
unnecessarily in versions for long double). This patch make such code
use the recently-added math_narrow_eval macro consistently, removing
the unnecessary uses of volatile in long double lrint / llrint
implementations completely.
Tested for x86_64, x86, mips64 and powerpc.
* math/s_nexttowardf.c (__nexttowardf): Use math_narrow_eval.
* stdlib/strtod_l.c: Include <math_private.h>.
(overflow_value): Use math_narrow_eval.
(underflow_value): Likewise.
* sysdeps/i386/fpu/s_nexttoward.c (__nexttoward): Likewise.
* sysdeps/i386/fpu/s_nexttowardf.c (__nexttowardf): Likewise.
* sysdeps/ieee754/dbl-64/e_gamma_r.c (gamma_positive): Likewise.
(__ieee754_gamma_r): Likewise.
* sysdeps/ieee754/dbl-64/gamma_productf.c (__gamma_productf):
Likewise.
* sysdeps/ieee754/dbl-64/k_rem_pio2.c (__kernel_rem_pio2):
Likewise.
* sysdeps/ieee754/dbl-64/lgamma_neg.c (__lgamma_neg): Likewise.
* sysdeps/ieee754/dbl-64/s_erf.c (__erfc): Likewise.
* sysdeps/ieee754/dbl-64/s_llrint.c (__llrint): Likewise.
* sysdeps/ieee754/dbl-64/s_lrint.c (__lrint): Likewise.
* sysdeps/ieee754/flt-32/e_gammaf_r.c (gammaf_positive): Likewise.
(__ieee754_gammaf_r): Likewise.
* sysdeps/ieee754/flt-32/k_rem_pio2f.c (__kernel_rem_pio2f):
Likewise.
* sysdeps/ieee754/flt-32/lgamma_negf.c (__lgamma_negf): Likewise.
* sysdeps/ieee754/flt-32/s_erff.c (__erfcf): Likewise.
* sysdeps/ieee754/flt-32/s_llrintf.c (__llrintf): Likewise.
* sysdeps/ieee754/flt-32/s_lrintf.c (__lrintf): Likewise.
* sysdeps/ieee754/ldbl-128/s_llrintl.c (__llrintl): Do not use
volatile.
* sysdeps/ieee754/ldbl-128/s_lrintl.c (__lrintl): Likewise.
* sysdeps/ieee754/ldbl-128/s_nexttoward.c (__nexttoward): Use
math_narrow_eval.
* sysdeps/ieee754/ldbl-128ibm/s_nexttoward.c (__nexttoward):
Likewise.
* sysdeps/ieee754/ldbl-128ibm/s_nexttowardf.c (__nexttowardf):
Likewise.
* sysdeps/ieee754/ldbl-96/gamma_product.c (__gamma_product):
Likewise.
* sysdeps/ieee754/ldbl-96/s_llrintl.c (__llrintl): Do not use
volatile.
* sysdeps/ieee754/ldbl-96/s_lrintl.c (__lrintl): Likewise.
* sysdeps/ieee754/ldbl-96/s_nexttoward.c (__nexttoward): Use
math_narrow_eval.
* sysdeps/ieee754/ldbl-96/s_nexttowardf.c (__nexttowardf):
Likewise.
* sysdeps/ieee754/ldbl-opt/s_nexttowardfd.c (__nldbl_nexttowardf):
Likewise.
Various i386 libm functions return values with excess range and
precision; Wilco Dijkstra's patches to make isfinite etc. expand
inline cause this pre-existing issue to result in test failures (when
e.g. a result that overflows float but not long double gets counted as
overflowing for some purposes but not others).
This patch addresses those cases arising from functions defined in C,
adding a math_narrow_eval macro that forces values to memory to
eliminate excess precision if FLT_EVAL_METHOD indicates this is
needed, and is a no-op otherwise. I'll convert existing uses of
volatile and asm for this purpose to use the new macro later, once
i386 has clean test results again (which requires fixes for .S files
as well).
Tested for x86_64 and x86. Committed.
[BZ #18980]
* sysdeps/generic/math_private.h: Include <float.h>.
(math_narrow_eval): New macro.
[FLT_EVAL_METHOD != 0] (excess_precision): Likewise.
* sysdeps/ieee754/dbl-64/e_cosh.c (__ieee754_cosh): Use
math_narrow_eval on overflowing return value.
* sysdeps/ieee754/dbl-64/e_lgamma_r.c (__ieee754_lgamma_r):
Likewise.
* sysdeps/ieee754/dbl-64/e_sinh.c (__ieee754_sinh): Likewise.
* sysdeps/ieee754/flt-32/e_coshf.c (__ieee754_coshf): Likewise.
* sysdeps/ieee754/flt-32/e_lgammaf_r.c (__ieee754_lgammaf_r):
Likewise.
* sysdeps/ieee754/flt-32/e_sinhf.c (__ieee754_sinhf): Likewise.
Bug 15384 notes that in __finite, two different constants are used
that could be the same constant (the result only depends on the
exponent of the floating-point representation), and that using the
same constant is better for architectures where constants need loading
from a constant pool. This patch implements that change.
Tested for x86_64, mips64 and powerpc.
[BZ #15384]
* sysdeps/ieee754/dbl-64/s_finite.c (FINITE): Use same constant as
bit-mask as in subtraction.
* sysdeps/ieee754/dbl-64/wordsize-64/s_finite.c (__finite):
Likewise.
* sysdeps/ieee754/flt-32/s_finitef.c (FINITEF): Likewise.
* sysdeps/ieee754/ldbl-128/s_finitel.c (__finitel): Likewise.
* sysdeps/ieee754/ldbl-128ibm/s_finitel.c (__finitel): Likewise.
Similar to various other bugs in this area, tgamma functions can fail
to raise the underflow exception when the result is tiny and inexact
but one or more low bits of the intermediate result that is scaled
down are zero. This patch forces the exception in a similar way to
previous fixes.
Tested for x86_64, x86, mips64 and powerpc.
[BZ #18951]
* sysdeps/ieee754/dbl-64/e_gamma_r.c (__ieee754_gamma_r): Force
underflow exception for small results.
* sysdeps/ieee754/flt-32/e_gammaf_r.c (__ieee754_gammaf_r):
Likewise.
* sysdeps/ieee754/ldbl-128/e_gammal_r.c (__ieee754_gammal_r):
Likewise.
* sysdeps/ieee754/ldbl-128ibm/e_gammal_r.c (__ieee754_gammal_r):
Likewise.
* sysdeps/ieee754/ldbl-96/e_gammal_r.c (__ieee754_gammal_r):
Likewise.
* math/auto-libm-test-in: Add more tests of tgamma.
* math/auto-libm-test-out: Regenerated.
As noted in bug 6803, scalbn fails to set errno on overflow and
underflow. This patch fixes this by making scalbn an alias of ldexp,
which has exactly the same semantics (for floating-point types with
radix 2) and already has wrappers that deal with setting errno,
instead of an alias of the internal __scalbn (which ldexp calls).
Notes:
* Where compat symbols were defined for scalbn functions, I didn't
change what they point to (to keep the patch minimal), so such
compat symbols continue to go directly to the non-errno-setting
functions.
* Mike, I didn't do anything with the IA64 versions of these
functions, where I think both the ldexp and scalbn functions already
deal with setting errno. As a cleanup (not needed to fix this bug)
however you might want to make those functions into aliases for
IA64; there is no need for them to be separate function
implementations at all.
* This concludes the fix for bug 6803 since the scalb and scalbln
cases of that bug were fixed some time ago.
Tested for x86_64, x86, mips64 and powerpc.
[BZ #6803]
* math/s_ldexp.c (scalbn): Define as weak alias of __ldexp.
[NO_LONG_DOUBLE] (scalbnl): Define as weak alias of __ldexp.
* math/s_ldexpf.c (scalbnf): Define as weak alias of __ldexpf.
* math/s_ldexpl.c (scalbnl): Define as weak alias of __ldexpl.
* sysdeps/i386/fpu/s_scalbn.S (scalbn): Remove alias.
* sysdeps/i386/fpu/s_scalbnf.S (scalbnf): Likewise.
* sysdeps/i386/fpu/s_scalbnl.S (scalbnl): Likewise.
* sysdeps/ieee754/dbl-64/s_scalbn.c (scalbn): Likewise.
[NO_LONG_DOUBLE] (scalbnl): Likewise.
* sysdeps/ieee754/dbl-64/wordsize-64/s_scalbn.c (scalbn):
Likewise.
[NO_LONG_DOUBLE] (scalbnl): Likewise.
* sysdeps/ieee754/flt-32/s_scalbnf.c (scalbnf): Likewise.
* sysdeps/ieee754/ldbl-128/s_scalbnl.c (scalbnl): Likewise.
* sysdeps/ieee754/ldbl-128ibm/s_scalbnl.c (scalbnl): Remove
long_double_symbol calls.
* sysdeps/ieee754/ldbl-64-128/s_scalbnl.c (scalbnl): Likewise.
* sysdeps/ieee754/ldbl-opt/s_ldexpl.c (__ldexpl_2): Define as
strong alias of __ldexpl.
(scalbnl): Define using long_double_symbol.
* sysdeps/m68k/m680x0/fpu/s_scalbn.c (__CONCATX(scalbn,suffix)):
Remove alias.
* sysdeps/sparc/sparc64/soft-fp/s_scalbnl.c (scalbnl): Likewise.
* sysdeps/x86_64/fpu/s_scalbnl.S (scalbnl): Likewise.
* math/libm-test.inc (scalbn_test_data): Add errno expectations.
(scalbln_test_data): Add more errno expectations.
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.
The existing implementations of lgamma functions (except for the ia64
versions) use the reflection formula for negative arguments. This
suffers large inaccuracy from cancellation near zeros of lgamma (near
where the gamma function is +/- 1).
This patch fixes this inaccuracy. For arguments above -2, there are
no zeros and no large cancellation, while for sufficiently large
negative arguments the zeros are so close to integers that even for
integers +/- 1ulp the log(gamma(1-x)) term dominates and cancellation
is not significant. Thus, it is only necessary to take special care
about cancellation for arguments around a limited number of zeros.
Accordingly, this patch uses precomputed tables of relevant zeros,
expressed as the sum of two floating-point values. The log of the
ratio of two sines can be computed accurately using log1p in cases
where log would lose accuracy. The log of the ratio of two gamma(1-x)
values can be computed using Stirling's approximation (the difference
between two values of that approximation to lgamma being computable
without computing the two values and then subtracting), with
appropriate adjustments (which don't reduce accuracy too much) in
cases where 1-x is too small to use Stirling's approximation directly.
In the interval from -3 to -2, using the ratios of sines and of
gamma(1-x) can still produce too much cancellation between those two
parts of the computation (and that interval is also the worst interval
for computing the ratio between gamma(1-x) values, which computation
becomes more accurate, while being less critical for the final result,
for larger 1-x). Because this can result in errors slightly above
those accepted in glibc, this interval is instead dealt with by
polynomial approximations. Separate polynomial approximations to
(|gamma(x)|-1)(x-n)/(x-x0) are used for each interval of length 1/8
from -3 to -2, where n (-3 or -2) is the nearest integer to the
1/8-interval and x0 is the zero of lgamma in the relevant half-integer
interval (-3 to -2.5 or -2.5 to -2).
Together, the two approaches are intended to give sufficient accuracy
for all negative arguments in the problem range. Outside that range,
the previous implementation continues to be used.
Tested for x86_64, x86, mips64 and powerpc. The mips64 and powerpc
testing shows up pre-existing problems for ldbl-128 and ldbl-128ibm
with large negative arguments giving spurious "invalid" exceptions
(exposed by newly added tests for cases this patch doesn't affect the
logic for); I'll address those problems separately.
[BZ #2542]
[BZ #2543]
[BZ #2558]
* sysdeps/ieee754/dbl-64/e_lgamma_r.c (__ieee754_lgamma_r): Call
__lgamma_neg for arguments from -28.0 to -2.0.
* sysdeps/ieee754/flt-32/e_lgammaf_r.c (__ieee754_lgammaf_r): Call
__lgamma_negf for arguments from -15.0 to -2.0.
* sysdeps/ieee754/ldbl-128/e_lgammal_r.c (__ieee754_lgammal_r):
Call __lgamma_negl for arguments from -48.0 or -50.0 to -2.0.
* sysdeps/ieee754/ldbl-96/e_lgammal_r.c (__ieee754_lgammal_r):
Call __lgamma_negl for arguments from -33.0 to -2.0.
* sysdeps/ieee754/dbl-64/lgamma_neg.c: New file.
* sysdeps/ieee754/dbl-64/lgamma_product.c: Likewise.
* sysdeps/ieee754/flt-32/lgamma_negf.c: Likewise.
* sysdeps/ieee754/flt-32/lgamma_productf.c: Likewise.
* sysdeps/ieee754/ldbl-128/lgamma_negl.c: Likewise.
* sysdeps/ieee754/ldbl-128/lgamma_productl.c: Likewise.
* sysdeps/ieee754/ldbl-128ibm/lgamma_negl.c: Likewise.
* sysdeps/ieee754/ldbl-128ibm/lgamma_productl.c: Likewise.
* sysdeps/ieee754/ldbl-96/lgamma_negl.c: Likewise.
* sysdeps/ieee754/ldbl-96/lgamma_product.c: Likewise.
* sysdeps/ieee754/ldbl-96/lgamma_productl.c: Likewise.
* sysdeps/generic/math_private.h (__lgamma_negf): New prototype.
(__lgamma_neg): Likewise.
(__lgamma_negl): Likewise.
(__lgamma_product): Likewise.
(__lgamma_productl): Likewise.
* math/Makefile (libm-calls): Add lgamma_neg and lgamma_product.
* math/auto-libm-test-in: Add more tests of lgamma.
* math/auto-libm-test-out: Regenerated.
* sysdeps/i386/fpu/libm-test-ulps: Update.
* sysdeps/x86_64/fpu/libm-test-ulps: Likewise.
Various fma implementations have logic that, when computing fma (x, y,
z) where z is large (so care needs taking to avoid internal overflow)
but x * y is small, scale x * y up instead of down to avoid internal
underflows resulting from scaling down. (In these cases, x * y is
small enough that only its sign actually matters rather than the exact
value.)
The threshold for scaling up instead of down was correct for "if the
unscaled values were multiplied, the low part of the multiplication
could underflow", and the scaling was sufficient to ensure that the
low part of the multiplication did not underflow (given that cases of
very small x * y - less than half the least subnormal - were
previously dealt with). However, the choice in the functions wasn't
between scaling up or no scaling, but between scaling up and scaling
down (scaling down actually being needed when x * y isn't so small
compared to z and so the exact value does matter). Thus a larger
threshold is needed to ensure that scaling down doesn't produce values
the multiplication of whose low parts underflows. This patch
increases the thresholds accordingly.
Tested for x86_64, x86 and mips64 (with the MIPS version of s_fmal.c
removed so that the ldbl-128 version gets tested instead of the
soft-fp one).
[BZ #18824]
* sysdeps/ieee754/dbl-64/s_fma.c (__fma): Increase threshold for
scaling x * y up instead of down.
* sysdeps/ieee754/ldbl-128/s_fmal.c (__fmal): Likewise.
* sysdeps/ieee754/ldbl-96/s_fmal.c (__fmal): Likewise.
* math/auto-libm-test-in: Add more tests of fma.
* math/auto-libm-test-out: Regenerated.
Similar to various other bugs in this area, some tanh implementations
do not raise the underflow exception for subnormal arguments, when the
result is tiny and inexact. This patch forces the exception in a
similar way to previous fixes.
Tested for x86_64, x86, mips64 and powerpc.
[BZ #16520]
* sysdeps/ieee754/dbl-64/s_tanh.c: Include <float.h>.
(__tanh): Force underflow exception for arguments with small
absolute value.
* sysdeps/ieee754/flt-32/s_tanhf.c: Include <float.h>.
(__tanhf): Force underflow exception for arguments with small
absolute value.
* sysdeps/ieee754/ldbl-128/s_tanhl.c: Include <float.h>.
(__tanhl): Force underflow exception for arguments with small
absolute value.
* sysdeps/ieee754/ldbl-128ibm/s_tanhl.c: Include <float.h>.
(__tanhl): Force underflow exception for arguments with small
absolute value.
* sysdeps/ieee754/ldbl-96/s_tanhl.c: Include <float.h>.
(__tanhl): Force underflow exception for arguments with small
absolute value.
* math/auto-libm-test-in: Add more tests of tanh.
* math/auto-libm-test-out: Regenerated.
* sysdeps/i386/fpu/libm-test-ulps: Update.
Similar to various other bugs in this area, some tan implementations
do not raise the underflow exception for subnormal arguments, when the
result is tiny and inexact. This patch forces the exception in a
similar way to previous fixes.
Tested for x86_64, x86, mips64 and powerpc.
[BZ #16517]
* sysdeps/ieee754/dbl-64/s_tan.c: Include <float.h>.
(tan): Force underflow exception for arguments with small absolute
value.
* sysdeps/ieee754/flt-32/k_tanf.c: Include <float.h>.
(__kernel_tanf): Force underflow exception for arguments with
small absolute value.
* sysdeps/ieee754/ldbl-128/k_tanl.c: Include <float.h>.
(__kernel_tanl): Force underflow exception for arguments with
small absolute value.
* sysdeps/ieee754/ldbl-128ibm/k_tanl.c: Include <float.h>.
(__kernel_tanl): Force underflow exception for arguments with
small absolute value.
* sysdeps/ieee754/ldbl-96/k_tanl.c: Include <float.h>.
(__kernel_tanl): Force underflow exception for arguments with
small absolute value.
* math/auto-libm-test-in: Add more tests of tan.
* math/auto-libm-test-out: Regenerated.
Similar to various other bugs in this area, some sinh implementations
do not raise the underflow exception for subnormal arguments, when the
result is tiny and inexact. This patch forces the exception in a
similar way to previous fixes.
Tested for x86_64, x86, mips64 and powerpc.
[BZ #16519]
* sysdeps/ieee754/dbl-64/e_sinh.c: Include <float.h>.
(__ieee754_sinh): Force underflow exception for arguments with
small absolute value.
* sysdeps/ieee754/flt-32/e_sinhf.c: Include <float.h>.
(__ieee754_sinhf): Force underflow exception for arguments with
small absolute value.
* sysdeps/ieee754/ldbl-128/e_sinhl.c: Include <float.h>.
(__ieee754_sinhl): Force underflow exception for arguments with
small absolute value.
* sysdeps/ieee754/ldbl-128ibm/e_sinhl.c: Include <float.h>.
(__ieee754_sinhl): Force underflow exception for arguments with
small absolute value.
* sysdeps/ieee754/ldbl-96/e_sinhl.c: Include <float.h>.
(__ieee754_sinhl): Force underflow exception for arguments with
small absolute value.
* math/auto-libm-test-in: Add more tests of sinh.
* math/auto-libm-test-out: Regenerated.
* sysdeps/i386/fpu/libm-test-ulps: Update.
In non-default rounding modes, tgamma can be slightly less accurate
than permitted by glibc's accuracy goals.
Part of the problem is error accumulation, addressed in this patch by
setting round-to-nearest for internal computations. However, there
was also a bug in the code dealing with computing pow (x + n, x + n)
where x + n is not exactly representable, providing another source of
error even in round-to-nearest mode; it was necessary to address both
bugs to get errors for all testcases within glibc's accuracy goals.
Given this second fix, accuracy in round-to-nearest mode is also
improved (hence regeneration of ulps for tgamma should be from scratch
- truncate libm-test-ulps or at least remove existing tgamma entries -
so that the expected ulps can be reduced).
Some additional complications also arose. Certain tgamma tests should
strictly, according to IEEE semantics, overflow or not depending on
the rounding mode; this is beyond the scope of glibc's accuracy goals
for any function without exactly-determined results, but
gen-auto-libm-tests doesn't handle being lax there as it does for
underflow. (libm-test.inc also doesn't handle being lax about whether
the result in cases very close to the overflow threshold is infinity
or a finite value close to overflow, but that doesn't cause problems
in this case though I've seen it cause problems with random test
generation for some functions.) Thus, spurious-overflow markings,
with a comment, are added to auto-libm-test-in (no bug in Bugzilla
because the issue is with the testsuite, not a user-visible bug in
glibc). And on x86, after the patch I saw ERANGE issues as previously
reported by Carlos (see my commentary in
<https://sourceware.org/ml/libc-alpha/2015-01/msg00485.html>), which
needed addressing by ensuring excess range and precision were
eliminated at various points if FLT_EVAL_METHOD != 0.
I also noticed and fixed a cosmetic issue where 1.0f was used in long
double functions and should have been 1.0L.
This completes the move of all functions to testing in all rounding
modes with ALL_RM_TEST, so gen-libm-have-vector-test.sh is updated to
remove the workaround for some functions not using ALL_RM_TEST.
Tested for x86_64, x86, mips64 and powerpc.
[BZ #18613]
* sysdeps/ieee754/dbl-64/e_gamma_r.c (gamma_positive): Take log of
X_ADJ not X when adjusting exponent.
(__ieee754_gamma_r): Do intermediate computations in
round-to-nearest then adjust overflowing and underflowing results
as needed.
* sysdeps/ieee754/flt-32/e_gammaf_r.c (gammaf_positive): Take log
of X_ADJ not X when adjusting exponent.
(__ieee754_gammaf_r): Do intermediate computations in
round-to-nearest then adjust overflowing and underflowing results
as needed.
* sysdeps/ieee754/ldbl-128/e_gammal_r.c (gammal_positive): Take
log of X_ADJ not X when adjusting exponent.
(__ieee754_gammal_r): Do intermediate computations in
round-to-nearest then adjust overflowing and underflowing results
as needed. Use 1.0L not 1.0f as numerator of division.
* sysdeps/ieee754/ldbl-128ibm/e_gammal_r.c (gammal_positive): Take
log of X_ADJ not X when adjusting exponent.
(__ieee754_gammal_r): Do intermediate computations in
round-to-nearest then adjust overflowing and underflowing results
as needed. Use 1.0L not 1.0f as numerator of division.
* sysdeps/ieee754/ldbl-96/e_gammal_r.c (gammal_positive): Take log
of X_ADJ not X when adjusting exponent.
(__ieee754_gammal_r): Do intermediate computations in
round-to-nearest then adjust overflowing and underflowing results
as needed. Use 1.0L not 1.0f as numerator of division.
* math/libm-test.inc (tgamma_test_data): Remove one test. Moved
to auto-libm-test-in.
(tgamma_test): Use ALL_RM_TEST.
* math/auto-libm-test-in: Add one test of tgamma. Mark some other
tests of tgamma with spurious-overflow.
* math/auto-libm-test-out: Regenerated.
* math/gen-libm-have-vector-test.sh: Do not check for START.
* sysdeps/i386/fpu/libm-test-ulps: Update.
* sysdeps/x86_64/fpu/libm-test-ulps: Likewise.
Similar to various other bugs in this area, j1 and jn implementations
can fail to raise the underflow exception when the internal
computation is exact although the actual function is inexact. This
patch forces the exception in a similar way to other such fixes. (The
ldbl-128 / ldbl-128ibm j1l implementation is different and doesn't
need a change for this until spurious underflows in it are fixed.)
Tested for x86_64, x86, mips64 and powerpc.
[BZ #16559]
* sysdeps/ieee754/dbl-64/e_j1.c: Include <float.h>.
(__ieee754_j1): Force underflow exception for small results.
* sysdeps/ieee754/dbl-64/e_jn.c (__ieee754_jn): Likewise.
* sysdeps/ieee754/flt-32/e_j1f.c: Include <float.h>.
(__ieee754_j1f): Force underflow exception for small results.
* sysdeps/ieee754/flt-32/e_jnf.c (__ieee754_jnf): Likewise.
* sysdeps/ieee754/ldbl-128/e_jnl.c (__ieee754_jnl): Likewise.
* sysdeps/ieee754/ldbl-128ibm/e_jnl.c (__ieee754_jnl): Likewise.
* sysdeps/ieee754/ldbl-96/e_j1l.c: Include <float.h>.
(__ieee754_j1l): Force underflow exception for small results.
* sysdeps/ieee754/ldbl-96/e_jnl.c (__ieee754_jnl): Likewise.
* math/auto-libm-test-in: Add more tests of j1 and jn.
* math/auto-libm-test-out: Regenerated.
Some existing jn tests, if run in non-default rounding modes, produce
errors above those accepted in glibc, which causes problems for moving
tests of jn to use ALL_RM_TEST. This patch makes jn set rounding
to-nearest internally, as was done for yn some time ago, then computes
the appropriate underflowing value for results that underflowed to
zero in to-nearest, and moves the tests to ALL_RM_TEST. It does
nothing about the general inaccuracy of Bessel function
implementations in glibc, though it should make jn more accurate on
average in non-default rounding modes through reduced error
accumulation. The recomputation of results that underflowed to zero
should as a side-effect fix some cases of bug 16559, where jn just
used an exact zero, but that is *not* the goal of this patch and other
cases of that bug remain unfixed.
(Most of the changes in the patch are reindentation to add new scopes
for SET_RESTORE_ROUND*.)
Tested for x86_64, x86, powerpc and mips64.
[BZ #16559]
[BZ #18602]
* sysdeps/ieee754/dbl-64/e_jn.c (__ieee754_jn): Set
round-to-nearest internally then recompute results that
underflowed to zero in the original rounding mode.
* sysdeps/ieee754/flt-32/e_jnf.c (__ieee754_jnf): Likewise.
* sysdeps/ieee754/ldbl-128/e_jnl.c (__ieee754_jnl): Likewise.
* sysdeps/ieee754/ldbl-128ibm/e_jnl.c (__ieee754_jnl): Likewise.
* sysdeps/ieee754/ldbl-96/e_jnl.c (__ieee754_jnl): Likewise
* math/libm-test.inc (jn_test): Use ALL_RM_TEST.
* sysdeps/i386/fpu/libm-test-ulps: Update.
* sysdeps/x86_64/fpu/libm-test-ulps: Likewise.
Similar to various other bugs in this area, some sin and sincos
implementations do not raise the underflow exception for subnormal
arguments, when the result is tiny and inexact. This patch forces the
exception in a similar way to previous fixes.
Tested for x86_64, x86, mips64 and powerpc.
[BZ #16526]
[BZ #16538]
* sysdeps/ieee754/dbl-64/s_sin.c: Include <float.h>.
(__sin): Force underflow exception for arguments with small
absolute value.
* sysdeps/ieee754/flt-32/k_sinf.c: Include <float.h>.
(__kernel_sinf): Force underflow exception for arguments with
small absolute value.
* sysdeps/ieee754/ldbl-128/k_sincosl.c: Include <float.h>.
(__kernel_sincosl): Force underflow exception for arguments with
small absolute value.
* sysdeps/ieee754/ldbl-128/k_sinl.c: Include <float.h>.
(__kernel_sinl): Force underflow exception for arguments with
small absolute value.
* sysdeps/ieee754/ldbl-128ibm/k_sincosl.c: Include <float.h>.
(__kernel_sincosl): Force underflow exception for arguments with
small absolute value.
* sysdeps/ieee754/ldbl-128ibm/k_sinl.c: Include <float.h>.
(__kernel_sinl): Force underflow exception for arguments with
small absolute value.
* sysdeps/ieee754/ldbl-96/k_sinl.c: Include <float.h>.
(__kernel_sinl): Force underflow exception for arguments with
small absolute value.
* sysdeps/powerpc/fpu/k_sinf.c: Include <float.h>.
(__kernel_sinf): Force underflow exception for arguments with
small absolute value.
* math/auto-libm-test-in: Add more tests of sin and sincos.
* math/auto-libm-test-out: Regenerated.
The dbl-64 and flt-32 implementations of exp2 functions produce
spurious underflow exceptions. The underlying reason is the same in
both cases: the computation works as (2^a - 1)*2^b + 2^b for suitably
chosen a and b, where a has small magnitude so 2^a - 1 can be computed
with a low-degree polynomial approximation, and (2^a - 1)*2^b can
underflow even when the final result does not. This patch fixes this
by adjusting the threshold for when scaling is used to avoid
intermediate underflow so it works for any possible value of a where
the final result would not underflow.
Tested for x86_64 and x86.
[BZ #18219]
* sysdeps/ieee754/dbl-64/e_exp2.c (__ieee754_exp2): Reduce
threshold on absolute value of exponent for which scaling is used.
* sysdeps/ieee754/flt-32/e_exp2f.c (__ieee754_exp2f): Likewise.
* math/auto-libm-test-in: Add more tests of exp2.
* math/auto-libm-test-out: Regenerated.
Similar to various other bugs in this area, some expm1 implementations
do not raise the underflow exception for subnormal arguments, when the
result is tiny and inexact. This patch forces the exception in a
similar way to previous fixes.
(The issue does not apply to the ldbl-* implementations or to those
for x86 / x86_64 long double. The change to
sysdeps/ieee754/dbl-64/wordsize-64/e_cosh.c is one I missed when
previously fixing bug 16354; the bug in that implementation was
previously latent, but the expm1 fixes stopped it being latent and so
required it to be fixed to avoid spurious underflows from cosh.)
Tested for x86_64 and x86.
[BZ #16353]
* sysdeps/i386/fpu/s_expm1.S (dbl_min): New object.
(__expm1): Force underflow exception for arguments with small
absolute value.
* sysdeps/i386/fpu/s_expm1f.S (flt_min): New object.
(__expm1f): Force underflow exception for arguments with small
absolute value.
* sysdeps/ieee754/dbl-64/s_expm1.c: Include <float.h>.
(__expm1): Force underflow exception for arguments with small
absolute value.
* sysdeps/ieee754/flt-32/s_expm1f.c: Include <float.h>.
(__expm1f): Force underflow exception for arguments with small
absolute value.
* sysdeps/ieee754/dbl-64/wordsize-64/e_cosh.c (__ieee754_cosh):
Check for small arguments before calling __expm1.
* math/auto-libm-test-in: Do not mark underflow exceptions as
possibly missing for bug 16353.
* math/auto-libm-test-out: Regenerated.
Similar to various other bugs in this area, some asinh implementations
do not raise the underflow exception for subnormal arguments, when the
result is tiny and inexact. This patch forces the exception in a
similar way to previous fixes.
Tested for x86_64, x86 and mips64.
[BZ #16350]
* sysdeps/i386/fpu/s_asinh.S (__asinh): Force underflow exception
for arguments with small absolute value.
* sysdeps/i386/fpu/s_asinhf.S (__asinhf): Likewise.
* sysdeps/i386/fpu/s_asinhl.S (__asinhl): Likewise.
* sysdeps/ieee754/dbl-64/s_asinh.c: Include <float.h>.
(__asinh): Force underflow exception for arguments with small
absolute value.
* sysdeps/ieee754/flt-32/s_asinhf.c: Include <float.h>.
(__asinhf): Force underflow exception for arguments with small
absolute value.
* sysdeps/ieee754/ldbl-128/s_asinhl.c: Include <float.h>.
(__asinhl): Force underflow exception for arguments with small
absolute value.
* sysdeps/ieee754/ldbl-128ibm/s_asinhl.c: Include <float.h>.
(__asinhl): Force underflow exception for arguments with small
absolute value.
* sysdeps/ieee754/ldbl-96/s_asinhl.c: Include <float.h>.
(__asinhl): Force underflow exception for arguments with small
absolute value.
* math/auto-libm-test-in: Do not mark underflow exceptions as
possibly missing for bug 16350.
* math/auto-libm-test-out: Regenerated.
