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.
I noticed that some of my recently added tests of lround and llround
wrongly expected the "inexact" exception to be absent for certain
within-range non-integer arguments. (It's unspecified whether this
exception is present or not for within-range non-integer arguments; it
mustn't be present for integer arguments and out-of-range arguments.)
This patch corrects those expectations.
Tested for x86_64 and x86.
* math/libm-test.inc (lround_test_data): Do not expect the absence
of "inexact" for some tests with non-integer arguments.
(llround_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 adds more tests of lrint, llrint, lround and llround, to
cover various standard special cases not previously covered, and more
tests of overflow.
Tested for x86_64 and x86.
* math/libm-test.inc (lrint_test_data): Add more tests.
(llrint_test_data): Likewise.
(lround_test_data): Likewise.
(llround_test_data): Likewise.
This patch makes lrint and llrint use the same test inputs in
libm-test.inc, appropriately conditioned on LONG_MAX in the lrint
case.
Tested for x86_64 and x86.
* math/libm-test.inc (lrint_test_data): Add tests used for llrint.
(llrint_test_data): Add tests used for lrint.
This patch adds more libm-test.inc expectations for the "inexact"
exception for scalb, in all cases except those with a non-integer
second argument (where results are unspecified by POSIX, so the
function does not count as fully determined and the spurious "inexact"
exceptions raised by the existing implementations alongside "invalid"
are OK).
Tested for x86_64 and x86.
* math/libm-test.inc (scalb_test_data): Add more expectations for
the "inexact" exception.
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.
The ldbl-128 / ldbl-128ibm implementation of lgamma has problems with
its handling of large arguments. It has an overflow threshold that is
correct only for ldbl-128, despite being used for both types - with
diagnostic control macros as a temporary measure to disable warnings
about that constant overflowing for ldbl-128ibm - and it has a
calculation that's roughly x * log(x) - x, resulting in overflows for
arguments that are roughly at most a factor 1/log(threshold) below the
overflow threshold.
This patch fixes both issues, using an overflow threshold appropriate
for the type in question and adding another case for large arguments
that avoids the possible intermediate overflow.
Tested for x86_64, x86, mips64 and powerpc.
[BZ #16347]
[BZ #19046]
* sysdeps/ieee754/ldbl-128/e_lgammal_r.c: Do not include
<libc-internal.h>.
(MAXLGM): Do not use diagnostic control macros.
[LDBL_MANT_DIG == 106] (MAXLGM): Change value to overflow
threshold for ldbl-128ibm.
(__ieee754_lgammal_r): For large arguments, multiply by log - 1
instead of multiplying by log then subtracting.
* math/auto-libm-test-in: Add more tests of lgamma.
* math/auto-libm-test-out: Regenerated.
When libm-test.inc prints the results of failing tests, the output can
be unhelpful for ldbl-128 and ldbl-128ibm because the precision used
is insufficient to distinguish values of those types, resulting in
reported values that look identical but differ by a large number of
ulps.
This patch changes it to use a precision appropriate for the type, for
both decimal and hex output (so output for float is more compact,
output for ldbl-128 and ldbl-128ibm is substantially wider). The
natural precision to use for decimal is given by the C11 <float.h>
macros such as FLT_DECIMAL_DIG. GCC's <float.h> only defines those in
C11 mode, so this patch uses the predefines such as
__FLT_DECIMAL_DIG__ (added in GCC 4.6) instead; if we move to building
with -std=gnu11 (or -std=gnu1x if we can't get rid of 4.6 support).
Tested for powerpc and mips64.
* math/libm-test.inc (TYPE_DECIMAL_DIG): New macro.
(TYPE_HEX_DIG): Likewise.
(print_float): Use TYPE_DECIMAL_DIG - 1 and TYPE_HEX_DIG - 1 as
precisions when printing floating-point numbers.
(check_float_internal): Likewise.
The i386 versions of acoshf and acosh raise a spurious "invalid"
exception for an argument that is a quiet NaN with the sign bit set.
The integer arithmetic to detect arguments < 1 also detects -NaN, and
then the computation 0 / 0 in that case raises the exception. This
patch fixes this by using (x - x) / (x - x) as the computation in that
case instead, which will always raise the exception for non-NaN
arguments reaching that code, but not for quiet NaN arguments.
Tested for x86_64 and x86.
[BZ #19032]
* sysdeps/i386/fpu/e_acosh.S (__ieee754_acosh): For arguments < 1,
compute result as (x - x) / (x - x) not as 0 / 0.
* sysdeps/i386/fpu/e_acoshf.S (__ieee754_acoshf): Likewise.
* math/libm-test.inc (acosh_test_data): Add another test of acosh.
This patch improves test coverage of the real libm functions [a-e]*,
ensuring that special cases and ranges of input values of potential
significance (such as close to overflow and underflow thresholds) are
more systematically covered.
This is a followup to
<https://sourceware.org/ml/libc-alpha/2013-12/msg00757.html> which
covered [a-c]* (however, I found more weaknesses in the coverage of
those functions when preparing this patch, hence the additional tests
being added for them here).
Addition of a test for acosh (-qNaN) is temporarily deferred, to be
included as part of a fix for bug 19032 which was discovered in the
course of adding these tests (and which illustrates the use of testing
-qNaN as well as +qNaN as input even to functions for which the sign
of a NaN isn't meant to be significant).
Tested for x86_64 and x86.
* math/auto-libm-test-in: Add more tests of acos, acosh, asin,
atan, atan2, atanh, cbrt, cos, cosh, erf, erfc, exp, exp10, exp2
and expm1.
* math/auto-libm-test-out: Regenerated.
* math/libm-test.inc (acos_test_data): Add more tests.
(asin_test_data): Likewise.
(asinh_test_data): Likewise.
(atan_test_data): Likewise.
(atanh_test_data): Likewise.
(atan2_test_data): Likewise.
(cbrt_test_data): Likewise.
(ceil_test_data): Likewise.
(copysign_test_data): Likewise.
(cos_test_data): Likewise.
(cosh_test_data): Likewise.
(erf_test_data): Likewise.
(erfc_test_data): Likewise.
(exp_test_data): Likewise.
(exp10_test_data): Likewise.
(exp2_test_data): Likewise.
(expm1_test_data): Likewise.
* sysdeps/x86_64/fpu/libm-test-ulps: Update.
This patch makes math/libm-test.inc more consistent regarding
including expectations for errno setting and "inexact" exceptions
where expected test results are given manually. Mostly this is a
matter of including ERRNO_UNCHANGED in expectations, but there are
also some cases where expectations regarding "inexact" were missing
for exactly determined functions (especially in cases where some other
exception was expected and it should also have been expected that
"inexact" was not set with that other exception), and one case for pow
where the NO_INEXACT_EXCEPTION expectation should not have been there
(the rule about not having "inexact" exceptions for NaN arguments is
only when those NaN arguments produce NaN results).
I deferred making such changes for complex functions and scalb.
Tested for x86_64 and x86.
* math/libm-test.inc (acos_test_data): Refine expectations for
errno and "inexact" exceptions.
(acosh_test_data): Likewise.
(asin_test_data): Likewise.
(asinh_test_data): Likewise.
(atan_test_data): Likewise.
(atanh_test_data): Likewise.
(atan2_test_data): Likewise.
(cbrt_test_data): Likewise.
(ceil_test_data): Likewise.
(copysign_test_data): Likewise.
(cosh_test_data): Likewise.
(erf_test_data): Likewise.
(erfc_test_data): Likewise.
(exp_test_data): Likewise.
(exp10_test_data): Likewise.
(exp2_test_data): Likewise.
(expm1_test_data): Likewise.
(fabs_test_data): Likewise.
(floor_test_data): Likewise.
(fma_test_data): Likewise.
(fmax_test_data): Likewise.
(fmin_test_data): Likewise.
(fmod_test_data): Likewise.
(fpclassify_test_data): Likewise.
(frexp_test_data): Likewise.
(hypot_test_data): Likewise.
(ilogb_test_data): Likewise.
(isgreater_test_data): Likewise.
(isgreaterequal_test_data): Likewise.
(isinf_test_data): Likewise.
(isless_test_data): Likewise.
(islessequal_test_data): Likewise.
(islessgreater_test_data): Likewise.
(isnan_test_data): Likewise.
(isnormal_test_data): Likewise.
(issignaling_test_data): Likewise.
(isunordered_test_data): Likewise.
(j0_test_data): Likewise.
(j1_test_data): Likewise.
(jn_test_data): Likewise.
(lgamma_test_data): Likewise.
(lrint_test_data): Likewise.
(llrint_test_data): Likewise.
(log_test_data): Likewise.
(log10_test_data): Likewise.
(log1p_test_data): Likewise.
(log2_test_data): Likewise.
(logb_test_data): Likewise.
(lround_test_data): Likewise.
(llround_test_data): Likewise.
(modf_test_data): Likewise.
(nearbyint_test_data): Likewise.
(nextafter_test_data): Likewise.
(nexttoward_test_data): Likewise.
(pow_test_data): Likewise.
(remainder_test_data): Likewise.
(remquo_test_data): Likewise.
(rint_test_data): Likewise.
(round_test_data): Likewise.
(signbit_test_data): Likewise.
(sinh_test_data): Likewise.
(sqrt_test_data): Likewise.
(tanh_test_data): Likewise.
(tgamma_test_data): Likewise.
(trunc_test_data): Likewise.
(y0_test_data): Likewise.
(y1_test_data): Likewise.
(yn_test_data): Likewise.
(significand_test_data): Likewise.
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.
The flt-32 version of powf can be inaccurate because of bugs in the
extra-precision calculation of (x-1)/(x+1) or (x-1.5)/(x+1.5) as part
of calculating log(x) with extra precision: a constant used (as part
of adding 1 or 1.5 through integer arithmetic) is incorrect, and then
the code fails to mask a computed high part before using it in
arithmetic that relies on s_h*t_h being exactly representable. This
patch fixes these bugs.
Tested for x86_64 and x86. x86_64 ulps for powf removed and
regenerated to reflect reduced ulps from the increased accuracy for
existing tests.
[BZ #18956]
* sysdeps/ieee754/flt-32/e_powf.c (__ieee754_powf): Add 0x00400000
not 0x0040000 for high bit of mantissa. Mask with 0xfffff000 when
extracting high part.
* math/auto-libm-test-in: Add another test of pow.
* math/auto-libm-test-out: Regenerated.
* sysdeps/x86_64/fpu/libm-test-ulps: Update.
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.
built-ins when available. Since going through the PLT is expensive for these small functions,
inlining results in major speedups (about 7x on Cortex-A57 for isinf). The GCC built-ins are not
correct if signalling NaN support is required, and thus are turned off in that case (see GCC bug
66462). The test-snan.c tests sNaNs and so must be explicitly built with -fsignaling-nans.
2015-09-18 Wilco Dijkstra <wdijkstr@arm.com>
[BZ #15367]
[BZ #17441]
* math/Makefile: Build test-snan.c with -fsignaling-nans.
* math/math.h (fpclassify): Use __builtin_fpclassify when
available. (signbit): Use __builtin_signbit(f/l).
(isfinite): Use__builtin_isfinite. (isnormal): Use
__builtin_isnormal. (isnan): Use __builtin_isnan.
(isinf): Use __builtin_isinf_sign.
C99/C11 Annex G specifies the sign of the zero part of the result of
ctan (x +/- i * Inf) and ctanh (+/-Inf + i * y). This patch fixes glibc
to follow that specification, along the lines I described in my review
of Andreas's previous patch for this issue
<https://sourceware.org/ml/libc-alpha/2014-08/msg00142.html>.
Tested for x86_64.
2015-09-17 Joseph Myers <joseph@codesourcery.com>
Andreas Schwab <schwab@suse.de>
[BZ #17118]
* math/s_ctan.c (__ctan): Determine sign of zero real part of
result when imaginary part of argument is infinite using sine and
cosine.
* math/s_ctanf.c (__ctanf): Likewise.
* math/s_ctanl.c (__ctanl): Likewise.
* math/s_ctanh.c (__ctanh): Determine sign of zero imaginary part
of result when real part of argument is infinite using sine and
cosine.
* math/s_ctanhf.c (__ctanhf): Likewise.
* math/s_ctanhl.c (__ctanhl): Likewise.
* math/libm-test.inc (ctan_test_data): Add more tests of ctan.
(ctanh_test_data): Add more tests of ctanh.
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.
math.h incorrectly declares various functions for XSI POSIX 2001 and
2008 editions. gamma was removed in the 2001 edition but is still
declared, along with gammaf and gammal which were never standard
functions. isnan is still declared as a function, along with isnanf
and isnanl which were never standard functions, although in 2001 the
function was replaced by the type-generic macro. scalbf and scalbl
are declared although never standard, and scalb was removed in the
2008 edition but is still declared. The scalb type-generic macro in
tgmath.h shouldn't be present for any POSIX version, since POSIX never
had such a type-generic macro.
This patch disables all those declarations in the relevant cases (as a
minimal fix, it leaves them enabled for __USE_MISC). For the matter
of declaring scalb but not scalbf or scalbl for the 2001 edition, a
new macro __MATH_DECLARING_DOUBLE is added, defined by math.h around
includes of bits/mathcalls.h, for bits/mathcalls.h to use to test
which type's functions are being declared.
Tested for x86_64 and x86 (testsuite, and that installed stripped
shared libraries are unchanged by the patch).
[BZ #18967]
* math/math.h (__MATH_DECLARING_DOUBLE): New macro. Define and
undefine around includes of <bits/mathcalls.h>.
* math/bits/mathcalls.h [!__USE_MISC && __USE_XOPEN2K] (isnan): Do
not declare function.
[!__USE_MISC && __USE_XOPEN2K] (gamma): Likewise.
[!__USE_MISC && (!__MATH_DECLARING_DOUBLE || __USE_XOPEN2K8)]
(scalb): Likewise.
* math/tgmath.h [!__USE_MISC && __USE_XOPEN_EXTENDED] (scalb): Do
not define macro.
* conform/Makefile (test-xfail-XOPEN2K/math.h/conform): Remove
variable.
(test-xfail-XOPEN2K/tgmath.h/conform): Likewise.
(test-xfail-XOPEN2K8/math.h/conform): Likewise.
(test-xfail-XOPEN2K8/tgmath.h/conform): Likewise.
Bug 16296 notes that fegetround is a pure function and should be
marked as such in fenv.h. This patch implements that.
Tested for x86_64 and x86 (testsuite, and that installed stripped
shared libraries are unchanged by this patch).
[BZ #16296]
* math/fenv.h (fegetround): Use __attribute_pure__.
* include/fenv.h (__fegetround): Likewise.
Similar to various other bugs in this area, ctan and ctanh can fail to
raise the underflow exception for some cases of results that are tiny
and inexact. This patch forces the exception in a similar way to
previous fixes.
Tested for x86_64 and x86.
[BZ #18595]
* math/s_ctan.c (__ctan): Force underflow exception for results
whose real or imaginary part has small absolute value.
* math/s_ctanf.c (__ctanf): Likewise.
* math/s_ctanh.c (__ctanh): Likewise.
* math/s_ctanhf.c (__ctanhf): Likewise.
* math/s_ctanhl.c (__ctanhl): Likewise.
* math/s_ctanl.c (__ctanl): Likewise.
* math/auto-libm-test-in: Do not allow missing underflow for ctan
and ctanh. Add more tests of ctan and ctanh.
On i386, the double version of exp10 can miss underflow exceptions if
the result is in the subnormal range for double but the last 11 bits
of the 64-bit extended-precision mantissa happen to be zero. This
patch forces the exception in a similar way to previous fixes.
As with the exp2 and exp fixes, the exp10f changes may in fact not be
needed to ensure underflow exceptions, but are included for
consistency and to fix the exp10 part of bug 18875 by ensuring that
excess range and precision is removed from underflowing return values.
Tested for x86_64 and x86.
[BZ #18875]
[BZ #18966]
* sysdeps/i386/fpu/e_exp10.S (dbl_min): New object.
(MO): New macro.
(__ieee754_exp10): For small results, force underflow exception
and remove excess range and precision from return value.
* sysdeps/i386/fpu/e_exp10f.S (flt_min): New object.
(MO): New macro.
(__ieee754_exp10f): For small results, force underflow exception
and remove excess range and precision from return value.
* math/auto-libm-test-in: Add more tests of exp10.
* math/auto-libm-test-out: Regenerated.
On i386, the double version of exp can miss underflow exceptions if
the result is in the subnormal range for double but the last 11 bits
of the 64-bit extended-precision mantissa happen to be zero. This
patch forces the exception in a similar way to previous fixes.
As with the exp2 fixes, the expf changes may in fact not be needed to
ensure underflow exceptions, but are included for consistency and to
fix the exp part of bug 18875 by ensuring that excess range and
precision is removed from underflowing return values.
Tested for x86_64 and x86.
[BZ #18875]
[BZ #18961]
* sysdeps/i386/fpu/e_exp.S (dbl_min): New object.
(MO): New macro.
(__ieee754_exp): For small results, force underflow exception and
remove excess range and precision from return value.
(__exp_finite): Likewise.
* sysdeps/i386/fpu/e_expf.S (flt_min): New object.
(MO): New macro.
(__ieee754_expf): For small results, force underflow exception and
remove excess range and precision from return value.
(__expf_finite): Likewise.
* math/auto-libm-test-in: Add more tests of exp.
* math/auto-libm-test-out: Regenerated.
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.
This patch adds more libm test inputs found through random test
generation to increase previously known ulps. This particular test
generation was run for mips64, so most of the increased ulps are for
ldbl-128 (float and double having been fairly well covered by such
testing for x86_64), but there's the odd ulps increase for other
formats.
Tested for x86_64, x86 and mips64.
* math/auto-libm-test-in: Add more tests of acos, acosh, asin,
asinh, atan, atan2, atanh, cabs, carg, cos, csqrt, erfc, exp,
exp10, exp2, log, log1p, log2, pow, sin, sincos, sinh, tan and
tanh.
* math/auto-libm-test-out: Regenerated.
* sysdeps/i386/fpu/libm-test-ulps: Update.
* sysdeps/mips/mips32/libm-test-ulps: Likewise.
* sysdeps/mips/mips64/libm-test-ulps: Likewise.
* sysdeps/x86_64/fpu/libm-test-ulps: Likewise.
The ldbl-128 / ldbl-128ibm implementation of lgammal converts (the
floor of minus) non-integer negative arguments to int to determine the
value of signgam. When those values are outside the range of int,
this produces spurious "invalid" exceptions and incorrect values of
signgam. This patch fixes this by instead determining signgam through
comparing half the integer in question to floor of half the integer.
Tested for mips64, x86_64 and x86.
[BZ #18952]
* sysdeps/ieee754/ldbl-128/e_lgammal_r.c (__ieee754_lgammal_r): Do
not convert non-integer negative arguments to int to determine the
value of signgam.
* math/auto-libm-test-in: Add more tests of lgamma.
* math/auto-libm-test-out: Regenerated.
This patch adds more libm test inputs found through random test
generation to increase observed ulps on x86_64.
Tested for x86_64 and x86.
* math/auto-libm-test-in: Add more tests of acosh, atanh, cbrt,
cosh, csqrt, erfc, expm1 and lgamma.
* math/auto-libm-test-out: Regenerated.
* sysdeps/i386/fpu/libm-test-ulps: Update.
* sysdeps/x86_64/fpu/libm-test-ulps: Likewise.
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.
The uninitialized variable warnings in math/ having been fixed for all
the supported floating-point formats, this patch removes the use of
-Wno-uninitialized there, continuing with the goal of avoiding -Wno-
options in makefiles as far as possible..
Tested for x86_64 and x86 (full build and testsuite runs), and for
powerpc and mips64 (verified that glibc builds without errors).
* math/Makefile (CFLAGS): Don't add -Wno-uninitialized.
The csqrt implementations in glibc can miss underflow exceptions when
the real or imaginary part of the result becomes tiny in the course of
scaling down (in particular, multiplication by 0.5) and that scaling
is exact although the relevant part of the mathematical result isn't.
This patch forces the exception in a similar way to previous fixes.
Tested for x86_64 and x86.
[BZ #18370]
* math/s_csqrt.c (__csqrt): Force underflow exception for results
whose real or imaginary part has small absolute value.
* math/s_csqrtf.c (__csqrtf): Likewise.
* math/s_csqrtl.c (__csqrtl): Likewise.
* math/auto-libm-test-in: Add more tests of csqrt.
* math/auto-libm-test-out: Regenerated.
* sysdeps/i386/fpu/libm-test-ulps: Update.
The csqrt functions scale up small arguments to avoid underflows when
calling hypot functions. However, even when hypot does not underflow,
a subsequent calculation of 0.5 * hypot can underflow. This patch
duly increases the threshold and scale factor to avoid such underflows
as well.
Tested for x86_64, x86 and mips64.
[BZ #18823]
* math/s_csqrt.c (__csqrt): Increase threshold and scale factor
for scaling up small arguments.
* math/s_csqrtf.c (__csqrtf): Likewise.
* math/s_csqrtl.c (__csqrtl): Likewise.
* math/auto-libm-test-in: Add more tests of csqrt.
* math/auto-libm-test-out: Regenerated.
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.
This patch adds more test inputs to various libm functions found
through random generation to have larger ulps errors than previously
listed in libm-test-ulp, on at least one of x86_64 and x86.
Tested for x86_64 and x86.
* math/auto-libm-test-in: Add more tests of acos, acosh, asin,
asinh, atan, atan2, atanh, cabs, cbrt, cosh, csqrt, erf, erfc,
exp, exp2, lgamma, log, log1p, log2, pow, sin, sincos, tan, tanh
and tgamma.
* 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, 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.
This patch adds more tests of various libm functions found through
random test generation to give increased ulps on 32-bit x86.
Tested for x86_64 and x86.
* math/auto-libm-test-in: Add more tests of acosh, asin, asinh,
atanh, cabs, carg, cbrt, cosh, csqrt, erf, erfc, exp, exp10,
expm1, hypot, log, log10, log1p, log2, pow, sinh, tan and tgamma.
* math/auto-libm-test-out: Regenerated.
* sysdeps/i386/fpu/libm-test-ulps: Update.
* sysdeps/x86_64/fpu/libm-test-ulps: Likewise.
ldbl-128ibm tanhl uses a too-small threshold to decide when to return
+/-1, resulting in large errors. This patch changes it to a more
appropriate threshold (the requirement is for 2*exp(-2|x|) to be small
in terms of ulps of 1).
Tested for x86_64, x86 and powerpc.
[BZ #18790]
* sysdeps/ieee754/ldbl-128ibm/s_tanhl.c (__tanhl): Increase
threshold for returning +/- 1.
* math/auto-libm-test-in: Add more tests of tanh.
* math/auto-libm-test-out: Regenerated.
* sysdeps/i386/fpu/libm-test-ulps: Update.
ldbl-128ibm sinhl uses a too-big threshold to decide when to return
the argument, resulting in large errors. This patch fixes it to use a
more appropriate threshold.
Tested for x86_64, x86 and powerpc.
[BZ #18789]
* sysdeps/ieee754/ldbl-128ibm/e_sinhl.c (__ieee754_sinhl): Use
smaller threshold for returning the argument.
* math/auto-libm-test-in: Add more tests of sinh.
* 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.
