As with various other issues of this kind, bug 16977 is log10 (1)
wrongly returning -0 rather than +0 in round-downward mode because of
an implementation effectively in terms of log1p (x - 1). This patch
fixes the issue in the same way used for log.
Tested x86_64 and x86 and ulps updated accordingly. Also tested for
mips64 to confirm a fix was needed for ldbl-128 and to validate that
fix (also applied to ldbl-128ibm since that version of logl is
essentially the same as the ldbl-128 one).
[BZ #16977]
* sysdeps/i386/fpu/e_log10.S (__ieee754_log10): Take absolute
value when x - 1 is zero.
* sysdeps/i386/fpu/e_log10f.S (__ieee754_log10f): Likewise.
* sysdeps/i386/fpu/e_log10l.S (__ieee754_log10l): Likewise.
* sysdeps/ieee754/ldbl-128/e_log10l.c (__ieee754_log10l): Return
0.0L for an argument of 1.0L.
* sysdeps/ieee754/ldbl-128ibm/e_log10l.c (__ieee754_log10l):
Likewise.
* sysdeps/x86_64/fpu/e_log10l.S (__ieee754_log10l): Take absolute
value when x - 1 is zero.
* math/libm-test.inc (log10_test): Use ALL_RM_TEST.
* sysdeps/i386/fpu/libm-test-ulps: Update.
* sysdeps/x86_64/fpu/libm-test-ulps: Likewise.
Bug 16564 is spurious overflow of log1pl (LDBL_MAX) in FE_UPWARD mode,
resulting from log1pl adding 1 to its argument (for arguments not
close to 0), which overflows in that mode. This patch fixes this by
avoiding adding 1 to large arguments (precisely what counts as large
depends on the floating-point format).
Tested x86_64 and x86, and spot-checked log1pl tests on mips64 and
powerpc64.
[BZ #16564]
* sysdeps/i386/fpu/s_log1pl.S (__log1pl): Do not add 1 to positive
arguments with exponent 65 or above.
* sysdeps/ieee754/ldbl-128/s_log1pl.c (__log1pl): Do not add 1 to
arguments 0x1p113L or above.
* sysdeps/ieee754/ldbl-128ibm/s_log1pl.c (__log1pl): Do not add 1
to arguments 0x1p107L or above.
* sysdeps/x86_64/fpu/s_log1pl.S (__log1pl): Do not add 1 to
positive arguments with exponent 65 or above.
* math/auto-libm-test-in: Add more tests of log1p.
* math/auto-libm-test-out: Regenerated.
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.
According to ISO C Annex F, log (1) should be +0 in all rounding
modes, but some implementations in glibc wrongly return -0 in
round-downward mode (mapping to log1p (x - 1) is problematic because 1
- 1 is -0 in round-downward mode, and log1p (-0) is -0). This patch
fixes this. (It helps with some implementations of other functions
such as acosh, log2 and log10 that call out to log, but not enough to
enable all-rounding-modes testing for those functions without further
fixes to other implementations of them.)
Tested x86_64 and x86 and ulps updated accordingly, and did spot tests
for mips64 for the ldbl-128 fix, and i586 for the sysdeps/i386/fpu
implementations shadowed by those in sysdeps/i386/i686/fpu.
[BZ #16731]
* sysdeps/i386/fpu/e_log.S (__ieee754_log): Take absolute value
when x - 1 is zero.
* sysdeps/i386/fpu/e_logf.S (__ieee754_logf): Likewise.
* sysdeps/i386/fpu/e_logl.S (__ieee754_logl): Likewise.
* sysdeps/i386/i686/fpu/e_logl.S (__ieee754_logl): Likewise.
* sysdeps/ieee754/dbl-64/e_log.c (__ieee754_log): Return +0 when
argument is 1.
* sysdeps/ieee754/ldbl-128/e_logl.c (__ieee754_logl): Likewise.
* sysdeps/x86_64/fpu/e_logl.S: Take absolute value when x - 1 is
zero.
* math/libm-test.inc (log_test): Use ALL_RM_TEST.
* sysdeps/i386/fpu/libm-test-ulps: Update.
* sysdeps/x86_64/fpu/libm-test-ulps: Likewise.
ISO C requires the result of nextafter to be independent of the
rounding mode, even when underflow or overflow occurs. This patch
fixes the bug in various nextafter implementations that, having done
an overflowing computation to force an overflow exception (correct),
they then return the result of that computation rather than an
infinity computed some other way (incorrect, when the overflowing
result of arithmetic with that sign and rounding mode is finite but
the correct result is infinite) - generally by falling through to
existing code to return a value that in fact is correct for this case
(but was computed by an integer increment and so without generating
the exceptions required). Having fixed the bug, the previously
deferred conversion of nextafter testing in libm-test.inc to
ALL_RM_TEST is also included.
Tested x86_64 and x86; also spot-checked results of nextafter tests
for powerpc32 and mips64 to test the ldbl-128ibm and ldbl-128
changes. (The m68k change is untested.)
[BZ #16677]
* math/s_nextafter.c (__nextafter): Do not return value from
overflowing computation.
* sysdeps/i386/fpu/s_nextafterl.c (__nextafterl): Likewise.
* sysdeps/ieee754/flt-32/s_nextafterf.c (__nextafterf): Likewise.
* sysdeps/ieee754/ldbl-128/s_nextafterl.c (__nextafterl):
Likewise.
* sysdeps/ieee754/ldbl-128ibm/s_nextafterl.c (__nextafterl):
Likewise.
* sysdeps/m68k/m680x0/fpu/s_nextafterl.c (__nextafterl): Likewise.
* math/libm-test.inc (nextafter_test): Use ALL_RM_TEST.
This patch fixes bug 16400, spurious underflow exceptions for ldbl-128
/ ldbl-128ibm lgammal with small positive arguments, by just using
-__logl (x) as the result in the problem cases (similar to the
previous fix for problems with small negative arguments).
Tested powerpc32, and also tested on mips64 that this does not require
ulps regeneration for the ldbl-128 case.
* sysdeps/ieee754/ldbl-128/e_lgammal_r.c (__ieee754_lgammal_r):
Return -__logl (x) for small positive arguments without evaluating
a polynomial.
This patch fixes bug 16337, ldbl-128 lgammal spurious overflows for
small negative arguments (the arguments in question are already in the
testsuite). The implementation uses the reflection formula to compute
lgamma of negative x from lgamma of -x, effectively resulting in a
calculation -log(x^2) + log(-x); cancellation isn't problematic in
this case (bugs for problematic cancellation in lgamma are 2542, 2543,
2558), but the x^2 calculation can underflow (in which case there is
spurious logic to return an overflowing value - lgamma can only ever
correctly overflow for large positive arguments, though tgamma can
overflow for small arguments of either sign as well as large positive
arguments). The fix is simply to calculate the result directly with
logl when the argument is a small enough negative number.
Tested mips64.
* sysdeps/ieee754/ldbl-128/e_lgammal_r.c (__ieee754_lgammal_r):
Calculate results for small negative arguments directly rather
than using reflection formula with special underflow handling.
This patch fixes bug 16338, ldbl-128 logl not handling subnormals
(with consequent inaccuracy for lgammal as well). The fix is simply
to use __frexpl when determining the exponent, as done already in
log2l and log10l. Given the lack of testing of small arguments to any
of the log* functions, appropriate tests are added for all of them.
Tested x86_64 and x86 and ulps updated accordingly, and spot tests
also run for mips64 to confirm the ldbl-128 fix.
Note that while this fixes lgammal inaccuracy for small positive
arguments, I suspect that there will still be problems with spurious
underflows in that case.
* sysdeps/ieee754/ldbl-128/e_logl.c (__ieee754_logl): Use __frexpl
to determine exponent and adjust argument to have exponent of -1.
* math/auto-libm-test-in: Add more tests of log, log10, log1p and
log2.
* math/auto-libm-test-out: Regenerated.
* sysdeps/x86_64/fpu/libm-test-ulps: Update.
http://sourceware.org/ml/libc-alpha/2013-08/msg00081.html
This is the first of a series of patches to ban ieee854_long_double
and the ieee854_long_double macros when using IBM long double. union
ieee854_long_double just isn't correct for IBM long double, especially
when little-endian, and pretending it is OK has allowed a number of
bugs to remain undetected in sysdeps/ieee754/ldbl-128ibm/.
This changes the few places in generic code that use it.
* stdio-common/printf_size.c (__printf_size): Don't use
union ieee854_long_double in fpnum union.
* stdio-common/printf_fphex.c (__printf_fphex): Likewise. Use
signbit macro to retrieve sign from long double.
* stdio-common/printf_fp.c (___printf_fp): Use signbit macro to
retrieve sign from long double.
* sysdeps/ieee754/ldbl-128ibm/printf_fphex.c: Adjust for fpnum change.
* sysdeps/ieee754/ldbl-128/printf_fphex.c: Likewise.
* sysdeps/ieee754/ldbl-96/printf_fphex.c: Likewise.
* sysdeps/x86_64/fpu/printf_fphex.c: Likewise.
* math/test-misc.c (main): Don't use union ieee854_long_double.
ports/
* sysdeps/ia64/fpu/printf_fphex.c: Adjust for fpnum change.
With help from Joseph Myers.
* sysdeps/ieee754/ldbl-128/s_atanl.c (__atanl): Handle tiny and
very large arguments properly.
* math/libm-test.inc (atan_test): New tests.
(atan2_test): New tests.
* sysdeps/sparc/fpu/libm-test-ulps: Update.
* sysdeps/x86_64/fpu/libm-test-ulps: Update.
With help from Joseph Myers.
* sysdeps/ieee754/flt-32/e_j0f.c (__ieee754_y0f): Adjust tinyness
cutoff to 2**-13.
* sysdeps/ieee754/flt-32/e_j1f.c (__ieee754_y1f): Adjust tinyness
cutoff to 2**-25.
* sysdeps/ieee754/ldbl-128/e_j0l.c (U0): New constant.
( __ieee754_y0l): Avoid arithmetic underflow when 'x' is very
small.
* sysdeps/ieee754/ldbl-128/e_j1l.c (__ieee754_y1l): Likewise.
* math/libm-test.inc (y0_test): New tests.
(y1_test): New tests.
* sysdeps/i386/fpu/libm-test-ulps: Update.
* sysdeps/x86_64/fpu/libm-test-ulps: Update.
* sysdeps/sparc/fpu/libm-test-ulps: Update.
POSIX 2008 states that if the input for 'logb[f|l]' is a subnormal number
it should be treated as if it were normalized. This means the
implementation should calculate the log2 of the mantissa and add it to the
subnormal exponent (-126 for float and -1022 for double and IBM long
double). This patch takes care of that.
