nextafter and nexttoward fail to set errno on overflow and underflow.
This patch makes them do so in cases that should include all the cases
where such errno setting is required by glibc's goals for when to set
errno (but not all cases of underflow where the result is nonzero and
so glibc's goals do not require errno setting).
Tested for x86_64, x86, mips64 and powerpc.
[BZ #6799]
* math/s_nextafter.c: Include <errno.h>.
(__nextafter): Set errno on overflow and underflow.
* math/s_nexttowardf.c: Include <errno.h>.
(__nexttowardf): Set errno on overflow and underflow.
* sysdeps/i386/fpu/s_nextafterl.c: Include <errno.h>.
(__nextafterl): Set errno on overflow and underflow.
* sysdeps/i386/fpu/s_nexttoward.c: Include <errno.h>.
(__nexttoward): Set errno on overflow and underflow.
* sysdeps/i386/fpu/s_nexttowardf.c: Include <errno.h>.
(__nexttowardf): Set errno on overflow and underflow.
* sysdeps/ieee754/flt-32/s_nextafterf.c: Include <errno.h>.
(__nextafterf): Set errno on overflow and underflow.
* sysdeps/ieee754/ldbl-128/s_nextafterl.c: Include <errno.h>.
(__nextafterl): Set errno on overflow and underflow.
* sysdeps/ieee754/ldbl-128/s_nexttoward.c: Include <errno.h>.
(__nexttoward): Set errno on overflow and underflow.
* sysdeps/ieee754/ldbl-128/s_nexttowardf.c: Include <errno.h>.
(__nexttowardf): Set errno on overflow and underflow.
* sysdeps/ieee754/ldbl-128ibm/s_nextafterl.c: Include <errno.h>.
(__nextafterl): Set errno on overflow and underflow.
* sysdeps/ieee754/ldbl-128ibm/s_nexttoward.c: Include <errno.h>.
(__nexttoward): Set errno on overflow and underflow.
* sysdeps/ieee754/ldbl-128ibm/s_nexttowardf.c: Include <errno.h>.
(__nexttowardf): Set errno on overflow and underflow.
* sysdeps/ieee754/ldbl-96/s_nexttoward.c: Include <errno.h>.
(__nexttoward): Set errno on overflow and underflow.
* sysdeps/ieee754/ldbl-96/s_nexttowardf.c: Include <errno.h>.
(__nexttowardf): Set errno on overflow and underflow.
* sysdeps/ieee754/ldbl-opt/s_nexttowardfd.c: Include <errno.h>.
(__nldbl_nexttowardf): Set errno on overflow and underflow.
* sysdeps/m68k/m680x0/fpu/s_nextafterl.c: Include <errno.h>.
(__nextafterl): Set errno on overflow and underflow.
* math/libm-test.inc (nextafter_test_data): Do not allow errno
setting to be missing on overflow. Add more tests.
(nexttoward_test_data): Likewise.
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.
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.
The ldbl-96 version of lroundl is incorrect for systems with 64-bit
long when the argument's absolute value is just below a power of 2,
2^32 or more, and rounds up to the next integer; in such cases, it
returns 0. The problem is incrementing the high part of the mantissa
loses the high bit of the value (which is not an issue for any other
floating-point format, and is handled specially in lround when the bit
corresponding to 0.5 was in the high part rather than the low part).
This patch fixes this in a similar way to that used in llroundl:
storing the high part in an unsigned long variable before incrementing
it, so problems cannot occur in the case when this code is reachable.
I improved test coverage for both lround and llround by making them
use the same test inputs (appropriately conditioned on the size of
long in the lround case) - complete with the same comments, to make
comparison as easy as possible. (This test coverage improvement was
how I found the lroundl bug.)
Tested for x86_64 and x86.
[BZ #19071]
* sysdeps/ieee754/ldbl-96/s_lroundl.c (__lroundl): Use unsigned
long int variable to store possibly incremented high part of
mantissa.
* math/libm-test.inc (lround_test_data): Add tests used for
llround. Use [LONG_MAX > 0x7fffffff] consistently as condition
for tests requiring 64-bit long. Do not condition tests on
[TEST_FLOAT] unnecessarily.
(llround_test_data): Add tests used for lround. Add another
expectation for the "inexact" exception. Do not condition tests
on [TEST_FLOAT] unnecessarily.
ISO C requires overflowing results from nexttoward to be the
appropriate infinity independent of the rounding mode, but some
implementations use a rounding-mode-dependent result (this is the same
issue as was fixed for nextafter in bug 16677). This patch fixes the
problem by making the nexttoward implementations discard the result
from the floating-point computation that forced an overflow exception
and then return the infinity previously computed with integer
arithmetic.
Tested for x86_64, x86, mips64 and powerpc.
[BZ #19059]
* math/s_nexttowardf.c (__nexttowardf): Do not return value from
overflowing computation.
* sysdeps/i386/fpu/s_nexttoward.c (__nexttoward): Likewise.
* sysdeps/i386/fpu/s_nexttowardf.c (__nexttowardf): Likewise.
* sysdeps/ieee754/ldbl-128/s_nexttoward.c (__nexttoward):
Likewise.
* sysdeps/ieee754/ldbl-128/s_nexttowardf.c (__nexttowardf):
Likewise.
* sysdeps/ieee754/ldbl-128ibm/s_nexttoward.c (__nexttoward):
Likewise.
* sysdeps/ieee754/ldbl-128ibm/s_nexttowardf.c (__nexttowardf):
Likewise.
* sysdeps/ieee754/ldbl-96/s_nexttoward.c (__nexttoward): Likewise.
* sysdeps/ieee754/ldbl-96/s_nexttowardf.c (__nexttowardf):
Likewise.
* sysdeps/ieee754/ldbl-opt/s_nexttowardfd.c (__nldbl_nexttowardf):
Likewise.
* math/libm-test.inc (nexttoward_test_data): Add more tests.
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, 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.
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.
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.
As with other spots in the code, GCC 4.8 unnecessarily complains about
an uninitialized variable in tanl calcs, so this patch disables. With
it, the library and sees the usual set of test passes.
* sysdeps/ieee754/ldbl-96/k_tanl.c: Include <libc-internal.h>.
(__kernel_tanl): Ignore uninitialized warnings around use of SIGN.
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.
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.
If you remove the "override CFLAGS += -Wno-uninitialized" in
math/Makefile, you get errors from lgamma implementations of the form:
../sysdeps/ieee754/dbl-64/e_lgamma_r.c: In function '__ieee754_lgamma_r':
../sysdeps/ieee754/dbl-64/e_lgamma_r.c:297:13: error: 'nadj' may be used uninitialized in this function [-Werror=maybe-uninitialized]
if(hx<0) r = nadj - r;
This is one of the standard kinds of false positive uninitialized
warnings: nadj is set under a certain condition, and then later used
under the same condition. This patch uses DIAG_* macros to suppress
the warning on the use of nadj. The ldbl-128 / ldbl-128ibm
implementation has a substantially different structure that avoids
this issue.
Tested for x86_64. (In fact this patch eliminates the need for that
-Wno-uninitialized on x86_64, but I want to test on more architectures
before removing it.)
* sysdeps/ieee754/dbl-64/e_lgamma_r.c: Include <libc-internal.h>.
(__ieee754_lgamma_r): Ignore uninitialized warnings around use of
NADJ.
* sysdeps/ieee754/flt-32/e_lgammaf_r.c: Include <libc-internal.h>.
(__ieee754_lgammaf_r): Ignore uninitialized warnings around use of
NADJ.
* sysdeps/ieee754/ldbl-96/e_lgammal_r.c: Include <libc-internal.h>.
(__ieee754_lgammal_r): Ignore uninitialized warnings around use of
NADJ.
ldbl-96 remquol wrongly handles the case where the first argument is
finite and the second infinite, because the check for the second
argument being a NaN fails to disregard the explicit high mantissa bit
and so wrongly interprets an infinity as being a NaN. This patch
fixes this by masking off that bit, and improves test coverage for
both remainder and remquo (various cases were missing tests, or, as in
the case of the bug, were tested only for one of the two functions).
Tested for x86_64 and x86.
[BZ #18244]
* sysdeps/ieee754/ldbl-96/s_remquol.c (__remquol): Ignore explicit
high mantissa bit when testing whether P is a NaN.
* math/libm-test.inc (remainder_test_data): Add more tests.
(remquo_test_data): Likewise.
Similar to various other bugs in this area, some atanh 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. (No change in this regard is needed
for the i386 implementation; special handling to force underflows in
these cases will only be needed there when the spurious underflows,
bug 18049, get fixed.)
Tested for x86_64, x86, powerpc and mips64.
[BZ #16352]
* sysdeps/i386/fpu/e_atanh.S (dbl_min): New object.
(__ieee754_atanh): Force underflow exception for results with
small absolute value.
* sysdeps/i386/fpu/e_atanhf.S (flt_min): New object.
(__ieee754_atanhf): Force underflow exception for results with
small absolute value.
* sysdeps/ieee754/dbl-64/e_atanh.c: Include <float.h>.
(__ieee754_atanh): Force underflow exception for results with
small absolute value.
* sysdeps/ieee754/flt-32/e_atanhf.c: Include <float.h>.
(__ieee754_atanhf): Force underflow exception for results with
small absolute value.
* sysdeps/ieee754/ldbl-128/e_atanhl.c: Include <float.h>.
(__ieee754_atanhl): Force underflow exception for results with
small absolute value.
* sysdeps/ieee754/ldbl-128ibm/e_atanhl.c: Include <float.h>.
(__ieee754_atanhl): Force underflow exception for results with
small absolute value.
* sysdeps/ieee754/ldbl-96/e_atanhl.c: Include <float.h>.
(__ieee754_atanhl): Force underflow exception for results with
small absolute value.
* math/auto-libm-test-in: Do not allow missing underflow
exceptions from atanh.
* math/auto-libm-test-out: Regenerated.
The threshold in ldbl-96 atanhl for when to return the argument,
0x1p-28, is a bit too big, and that in ldbl-128ibm atanhl is much too
big (the relevant condition being x^3/3 being < 0.5ulp of x),
resulting in errors a bit above the limits of those considered
acceptable in glibc in the ldbl-96 case, and in large errors in the
ldbl-128ibm case. This patch changes those implementations to use
more appropriate thresholds and adds tests around the thresholds for
various formats.
Tested for x86_64, x86 and powerpc. x86_64 and x86 ulps updated
accordingly.
[BZ #18046]
[BZ #18047]
* sysdeps/ieee754/ldbl-128ibm/e_atanhl.c (__ieee754_atanhl): Use
0x1p-56L as threshold for just returning the argument.
* sysdeps/ieee754/ldbl-96/e_atanhl.c (__ieee754_atanhl): Use
0x1p-32L as threshold for just returning the argument.
* math/auto-libm-test-in: Add more tests of atanh.
* math/auto-libm-test-out: Regenerated.
* sysdeps/i386/fpu/libm-test-ulps: Update.
* sysdeps/x86_64/fpu/libm-test-ulp: Likewise.
math/Makefile currently has:
# The fdlibm code generates a lot of these warnings but is otherwise clean.
override CFLAGS += -Wno-uninitialized
This is of course undesirable; warnings should be disabled as narrowly
as possible. To remove this override, we need to fix files that
generate such warnings, or put warning-disabling pragmas in them.
This patch does so for Bessel function implementations, one of the
cases that have the warnings if the override is removed. The warnings
arise because functions set pointer variables p and q only for certain
values of the function argument, then use them unconditionally. As
the static functions in question only get called for arguments that
satisfy the last condition in the if/else chain, the natural fix is to
change the last "else if" to just "else", which this patch does. (The
ldbl-128 / ldbl-128ibm implementation of these functions is
substantially different and looks like it already does use "else" in
the last case in the nearest corresponding code.)
Tested for x86_64 and x86.
* sysdeps/ieee754/dbl-64/e_j0.c (pzero): Change last case for
setting p and q from "else if" to "else".
(qzero): Likewise.
* sysdeps/ieee754/dbl-64/e_j1.c (pone): Likewise.
(qone): Likewise.
* sysdeps/ieee754/flt-32/e_j0f.c (pzerof): Likewise.
(qzerof): Likewise.
* sysdeps/ieee754/flt-32/e_j1f.c (ponef): Likewise.
(qonef): Likewise.
* sysdeps/ieee754/ldbl-96/e_j0l.c (pzero): Likewise.
(qzero): Likewise.
* sysdeps/ieee754/ldbl-96/e_j1l.c (pone): Likewise.
(qone): Likewise.
Similar to various other bugs in this area, some asin 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, powerpc and mips64.
[BZ #16351]
* sysdeps/i386/fpu/e_asin.S (dbl_min): New object.
(MO): New macro.
(__ieee754_asin): Force underflow exception for results with small
absolute value.
* sysdeps/i386/fpu/e_asinf.S (flt_min): New object.
(MO): New macro.
(__ieee754_asinf): Force underflow exception for results with
small absolute value.
* sysdeps/ieee754/dbl-64/e_asin.c: Include <float.h> and <math.h>.
(__ieee754_asin): Force underflow exception for results with small
absolute value.
* sysdeps/ieee754/flt-32/e_asinf.c: Include <float.h>.
(__ieee754_asinf): Force underflow exception for results with
small absolute value.
* sysdeps/ieee754/ldbl-128/e_asinl.c: Include <float.h>.
(__ieee754_asinl): Force underflow exception for results with
small absolute value.
* sysdeps/ieee754/ldbl-128ibm/e_asinl.c: Include <float.h>.
(__ieee754_asinl): Force underflow exception for results with
small absolute value.
* sysdeps/ieee754/ldbl-96/e_asinl.c: Include <float.h>.
(__ieee754_asinl): Force underflow exception for results with
small absolute value.
* sysdeps/x86_64/fpu/multiarch/e_asin.c [HAVE_FMA4_SUPPORT]:
Include <math.h>.
* math/auto-libm-test-in: Do not mark underflow exceptions as
possibly missing for bug 16351.
* math/auto-libm-test-out: Regenerated.
Various remquo implementations produce a zero remainder with the wrong
sign (a zero remainder should always have the sign of the first
argument, as specified in IEEE 754) in round-downward mode, resulting
from the sign of 0 - 0. This patch checks for zero results and fixes
their sign accordingly.
Tested for x86_64, x86, mips64 and powerpc.
[BZ #17987]
* sysdeps/ieee754/dbl-64/s_remquo.c (__remquo): Ensure sign of
zero result does not depend on the sign resulting from
subtraction.
* sysdeps/ieee754/dbl-64/wordsize-64/s_remquo.c (__remquo):
Likewise.
* sysdeps/ieee754/flt-32/s_remquof.c (__remquof): Likewise.
* sysdeps/ieee754/ldbl-128/s_remquol.c (__remquol): Likewise.
* sysdeps/ieee754/ldbl-128ibm/s_remquol.c (__remquol): Likewise.
* sysdeps/ieee754/ldbl-96/s_remquol.c (__remquol): Likewise.
* math/libm-test.inc (remquo_test_data): Add more tests.
Various remquo implementations, when computing the last three bits of
the quotient, have spurious overflows when 4 times the second argument
to remquo overflows. These overflows can in turn cause bad results in
rounding modes where that overflow results in a finite value. This
patch adds tests to avoid the problem multiplications in cases where
they would overflow, similar to those that control an earlier
multiplication by 8.
Tested for x86_64, x86, mips64 and powerpc.
[BZ #17978]
* sysdeps/ieee754/dbl-64/s_remquo.c (__remquo): Do not form
products 4 * y and 2 * y where those would overflow.
* sysdeps/ieee754/dbl-64/wordsize-64/s_remquo.c (__remquo):
Likewise.
* sysdeps/ieee754/flt-32/s_remquof.c (__remquof): Likewise.
* sysdeps/ieee754/ldbl-128/s_remquol.c (__remquol): Likewise.
* sysdeps/ieee754/ldbl-128ibm/s_remquol.c (__remquol): Likewise.
* sysdeps/ieee754/ldbl-96/s_remquol.c (__remquol): Likewise.
* math/libm-test.inc (remquo_test_data): Add more tests.
This patch makes sincos set errno to EDOM when passed an infinity,
similarly to sin and cos.
Tested for x86_64, x86, powerpc and mips64. I don't know if the
architecture-specific implementations for ia64 and m68k might need
corresponding fixes.
2015-02-11 Joseph Myers <joseph@codesourcery.com>
[BZ #15467]
* sysdeps/ieee754/dbl-64/s_sincos.c: Include <errno.h>.
(__sincos): Set errno to EDOM for infinite argument.
* sysdeps/ieee754/flt-32/s_sincosf.c: Include <errno.h>.
(SINCOSF_FUNC): Set errno to EDOM for infinite argument.
* sysdeps/ieee754/ldbl-128/s_sincosl.c: Include <errno.h>.
(__sincosl): Set errno to EDOM for infinite argument.
* sysdeps/ieee754/ldbl-128ibm/s_sincosl.c: Include <errno.h>.
(__sincosl): Set errno to EDOM for infinite argument.
* sysdeps/ieee754/ldbl-96/s_sincosl.c: Include <errno.h>.
(__sincosl): Set errno to EDOM for infinite argument.
* math/libm-test.inc (sincos_test_data): Test errno setting.
The ldbl-96 implementation of scalblnl (used for x86_64 and ia64) uses
a condition k <= -63 to determine when a standard underflowing result
tiny*__copysignl(tiny,x) should be returned. However, that condition
corresponds to values with exponent -16446 or less, and in the case of
-16446, the correct result for round-to-nearest depends on whether the
value is exactly 0x1p-16446 (half the least subnormal) or more than
that. This patch fixes the bug by changing the condition to k <= -64
and accordingly adjusting the exponent by 64 not 63 when converting to
a normal value.
Tested for x86_64.
[BZ #17803]
* sysdeps/ieee754/ldbl-96/s_scalblnl.c (twom63): Rename to
twom64. Adjust value to 0x1p-64L.
(__scalblnl): Only return standard underflowing result for K <=
-64 not K <= -63; adjust exponent for underflowing result by 64
not 63.
* math/libm-test.inc (scalbn_test_data): Add more tests.
(scalbln_test_data): Likewise.
The ldbl-96 implementation of scalblnl (used for x86_64 and ia64) is
incorrect for subnormal arguments (this is a separate bug from bug
17803, which is about underflowing results). There are two problems
with the adjustments of subnormal arguments: the "two63" variable
multiplied by is actually 0x1p52L not 0x1p63L, so is insufficient to
make values normal, and then GET_LDOUBLE_EXP(es,x), used to extract
the new exponent, extracts it into a variable that isn't used, while
the value taken to by the new exponent is wrongly taken from the high
part of the mantissa before the adjustment (hx). This patch fixes
both those problems and adds appropriate tests.
Tested for x86_64.
[BZ #17834]
* sysdeps/ieee754/ldbl-96/s_scalblnl.c (two63): Change value to
0x1p63L.
(__scalblnl): Get new exponent of adjusted subnormal value from ES
not HX.
* math/libm-test.inc (scalbn_test_data): Add more tests.
(scalbln_test_data): Likewise.
This patch fixes bugs 16561 and 16562, bad results of yn in overflow
cases in non-default rounding modes, both because an intermediate
overflow in the recurrence does not get detected if the result is not
an infinity and because an overflowing result may occur in the wrong
sign. The fix is to set FE_TONEAREST mode internally for the parts of
the function where such overflows can occur (which includes the call
to y1 - where yn is used to compute a Bessel function of order -1,
negating the result of y1 isn't correct for overflowing results in
directed rounding modes) and then compute an overflowing value in the
original rounding mode if the to-nearest result was an infinity.
Tested x86_64 and x86 and ulps updated accordingly. Also tested for
mips64 and powerpc32 to test the ldbl-128 and ldbl-128ibm changes.
(The tests for these bugs were added in my previous y1 patch, so the
only thing this patch has to do with the testsuite is enable yn
testing in all rounding modes.)
[BZ #16561]
[BZ #16562]
* sysdeps/ieee754/dbl-64/e_jn.c: Include <float.h>.
(__ieee754_yn): Set FE_TONEAREST mode internally and then
recompute overflowing results in original rounding mode.
* sysdeps/ieee754/flt-32/e_jnf.c: Include <float.h>.
(__ieee754_ynf): Set FE_TONEAREST mode internally and then
recompute overflowing results in original rounding mode.
* sysdeps/ieee754/ldbl-128/e_jnl.c: Include <float.h>.
(__ieee754_ynl): Set FE_TONEAREST mode internally and then
recompute overflowing results in original rounding mode.
* sysdeps/ieee754/ldbl-128ibm/e_jnl.c: Include <float.h>.
(__ieee754_ynl): Set FE_TONEAREST mode internally and then
recompute overflowing results in original rounding mode.
* sysdeps/ieee754/ldbl-96/e_jnl.c: Include <float.h>.
(__ieee754_ynl): Set FE_TONEAREST mode internally and then
recompute overflowing results in original rounding mode.
* sysdeps/i386/fpu/fenv_private.h [!__SSE2_MATH__]
(libc_feholdsetround_ctx): New macro.
* math/libm-test.inc (yn_test): Use ALL_RM_TEST.
* sysdeps/i386/fpu/libm-test-ulps: Update.
* sysdeps/x86_64/fpu/libm-test-ulps : Likewise.
This patch fixes bug 16354, spurious underflows from cosh when a tiny
argument is passed to expm1 and expm1 correctly underflows although
the final result of cosh should be 1. As noted in that bug, some
cases are latent because of expm1 implementations not raising
underflow (bug 16353), but all the implementations are fixed
similarly. They already contained checks for tiny arguments, but the
checks were too late to avoid underflow from expm1 (although they
would avoid underflow from subsequent squaring of the result of
expm1); they are moved before the expm1 calls.
The thresholds used for considering arguments tiny are not
particularly consistent in how they relate to the precision of the
floating-point format in question. They are, however, all sufficient
to ensure that the round-to-nearest result of cosh is indeed 1 below
the threshold (although sometimes they are smaller than necessary).
But the previous logic did not return 1, but the previously computed 1
+ expm1(abs(x)) value. And the thresholds in the ldbl-128 and
ldbl-128ibm code (0x1p-71L - I suspect 0x3f8b was intended in the code
instead of 0x3fb8 - and (roughly) 0x1p-55L) are not sufficient for
that value to be 1. So by moving the test for tiny arguments, and
consequently returning 1 directly now the expm1 value hasn't been
computed by that point, this patch also fixes bug 17061, the (large
number of ulps) inaccuracy for small arguments in those
implementations. Tests for that bug are duly added.
Tested x86_64 and x86 and ulps updated accordingly. Also tested for
mips64 and powerpc32 to validate the ldbl-128 and ldbl-128ibm changes.
[BZ #16354]
[BZ #17061]
* sysdeps/ieee754/dbl-64/e_cosh.c (__ieee754_cosh): Check for
small arguments before calling __expm1.
* sysdeps/ieee754/flt-32/e_coshf.c (__ieee754_coshf): Check for
small arguments before calling __expm1f.
* sysdeps/ieee754/ldbl-128/e_coshl.c (__ieee754_coshl): Check for
small arguments before calling __expm1l.
* sysdeps/ieee754/ldbl-128ibm/e_coshl.c (__ieee754_coshl):
Likewise.
* sysdeps/ieee754/ldbl-96/e_coshl.c (__ieee754_coshl): Likewise.
* math/auto-libm-test-in: Add more cosh tests. Do not allow
spurious underflow for some cosh tests.
* math/auto-libm-test-out: Regenerated.
* sysdeps/i386/fpu/libm-test-ulps: Update.
This patch fixes bug 17050, missing errno setting for y1 overflow (for
small positive arguments). An appropriate check is added for overflow
directly in the __ieee754_y1 implementation, similar to the check
present for yn (doing it there rather than in the wrapper also avoids
yn needing to repeat the check when called for order 1 or -1 and it
uses __ieee754_y1).
Tested x86_64 and x86; no ulps update needed. Also tested for mips64
to verify the ldbl-128 fix (the ldbl-128ibm code just #includes the
ldbl-128 file).
[BZ #17050]
* sysdeps/ieee754/dbl-64/e_j1.c: Include <errno.h>.
(__ieee754_y1): Set errno if return value overflows.
* sysdeps/ieee754/flt-32/e_j1f.c: Include <errno.h>.
(__ieee754_y1f): Set errno if return value overflows.
* sysdeps/ieee754/ldbl-128/e_j1l.c: Include <errno.h>.
(__ieee754_y1l): Set errno if return value overflows.
* sysdeps/ieee754/ldbl-96/e_j1l.c: Include <errno.h>.
(__ieee754_y1l): Set errno if return value overflows.
* math/auto-libm-test-in: Add more tests of y0, y1 and yn.
* math/auto-libm-test-out: Regenerated.
Errno is not set and the testcases will fail.
Now the scalbln-aliases are removed in i386/m68
and the wrappers are used when calling the scalbln-functions.
On ia64 only scalblnf has its own implementation.
For scalbln and scalblnl the ieee754/dbl-64 and ieee754/ldbl-96 are used, thus
the wrappers are needed, too.
According to C99 and C11 Annex F, acosh (1) should be +0 in all
rounding modes. However, some implementations in glibc wrongly return
-0 in round-downward mode (which is what you get if you end up
computing log1p (-0), via 1 - 1 being -0 in round-downward mode).
This patch fixes the problem implementations, by correcting the test
for an exact 1 value in the ldbl-96 implementation to allow for the
explicit high bit of the mantissa, and by inserting fabs instructions
in the i386 implementations; tests of acosh are duly converted to
ALL_RM_TEST. I believe all the other sysdeps/ieee754 implementations
are already OK (I haven't checked the ia64 versions, but if buggy then
that will be obvious from the results of test runs after this patch is
in).
Tested x86_64 and x86 and ulps updated accordingly.
[BZ #16927]
* sysdeps/i386/fpu/e_acosh.S (__ieee754_acosh): Use fabs on x-1
value.
* sysdeps/i386/fpu/e_acoshf.S (__ieee754_acoshf): Likewise.
* sysdeps/i386/fpu/e_acoshl.S (__ieee754_acoshl): Likewise.
* sysdeps/ieee754/ldbl-96/e_acoshl.c (__ieee754_acoshl): Correct
for explicit high bit of mantissa when testing for argument equal
to 1.
* math/libm-test.inc (acosh_test): Use ALL_RM_TEST.
* sysdeps/i386/fpu/libm-test-ulps: Update.
* sysdeps/x86_64/fpu/libm-test-ulps: Likewise.
Bug 16516 reports spurious underflows from erf (for all floating-point
types), when the result is close to underflowing but does not actually
underflow.
erf (x) is about (2/sqrt(pi))*x for x close to 0, so there are
subnormal arguments for which it does not underflow. The various
implementations do (x + efx*x) (for efx = 2/sqrt(pi) - 1), for greater
accuracy than if just using a single multiplication by an
approximation to 2/sqrt(pi) (effectively, this way there are a few
more bits in the approximation to 2/sqrt(pi)). This can introduce
underflows when efx*x underflows even though the final result does
not, so a scaled calculation with 8*efx is done in these cases - but 8
is not a big enough scale factor to avoid all such underflows. 16 is
(any underflows with a scale factor of 16 would only occur when the
final result underflows), so this patch changes the code to use that
factor. Rather than recomputing all the values of the efx8 variable,
it is removed, leaving it to the compiler's constant folding to
compute 16*efx. As such scaling can also lose underflows when the
final scaling down happens to be exact, appropriate checks are added
to ensure underflow exceptions occur when required in such cases.
Tested x86_64 and x86; no ulps updates needed. Also spot-checked for
powerpc32 and mips64 to verify the changes to the ldbl-128ibm and
ldbl-128 implementations.
[BZ #16516]
* sysdeps/ieee754/dbl-64/s_erf.c (efx8): Remove variable.
(__erf): Scale by 16 instead of 8 in potentially underflowing
case. Ensure exception if result actually underflows.
* sysdeps/ieee754/flt-32/s_erff.c (efx8): Remove variable.
(__erff): Scale by 16 instead of 8 in potentially underflowing
case. Ensure exception if result actually underflows.
* sysdeps/ieee754/ldbl-128/s_erfl.c: Include <float.h>.
(efx8): Remove variable.
(__erfl): Scale by 16 instead of 8 in potentially underflowing
case. Ensure exception if result actually underflows.
* sysdeps/ieee754/ldbl-128ibm/s_erfl.c: Include <float.h>.
(efx8): Remove variable.
(__erfl): Scale by 16 instead of 8 in potentially underflowing
case. Ensure exception if result actually underflows.
* sysdeps/ieee754/ldbl-96/s_erfl.c: Include <float.h>.
(efx8): Remove variable.
(__erfl): Scale by 16 instead of 8 in potentially underflowing
case. Ensure exception if result actually underflows.
* math/auto-libm-test-in: Add more tests of erf.
* math/auto-libm-test-out: Regenerated.
