about summary refs log tree commit diff
path: root/sysdeps/ieee754
Commit message (Collapse)AuthorAgeFilesLines
* Fix ldbl-128 expm1l spurious underflow (bug 16539).Joseph Myers2014-06-301-0/+5
| | | | | | | | | | | | | | | This patch fixes spurious underflows from ldbl-128 expm1l, as reported in <https://sourceware.org/ml/libc-alpha/2014-06/msg00835.html> and exposed by the tests added for such a bug in the x86 / x86-64 version. The bug and fix are essentially the same, so no separate bug is filed in Bugzilla. Tested for mips64. [BZ #16539] * sysdeps/ieee754/ldbl-128/s_expm1l.c: Include <float.h>. (__expm1l): Return argument unchanged when small but not subnormal.
* Fix ldbl-128 powl sign of result in overflow / underflow cases (bug 17097).Joseph Myers2014-06-291-13/+13
| | | | | | | | | | | | | | | | | | | | | | | | This patch fixes bug 17097, ldbl-128 powl producing overflowing / underflowing results with positive sign when the result should have been negative. This was shown up by the tests in non-default rounding modes added by my patch for bug 16315, but isn't actually limited to non-default rounding modes: rather, when rounding to nearest the wrappers produced a result with the correct sign and so always hid the bug unless -lieee was used to disable the wrappers. The problem is that in the cases where Y is large enough that the result overflows or underflows for X not very close to 1, but not large enough to overflow or underflow for all X != +/- 1 (in the latter case Y is always an even integer), a positive overflowing / underflowing result is always returned, rather than one with the correct sign. This patch moves the relevant part of computation of the sign earlier and returns a result of the correct sign. Tested for mips64. [BZ #17097] * sysdeps/ieee754/ldbl-128/e_powl.c (__ieee754_powl): Return result with correct sign in case of exponents that produce overflow except for X very close to 1.
* Fix yn overflow handling in non-default rounding modes (bug 16561, bug 16562).Joseph Myers2014-06-275-248/+304
| | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | 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.
* Fix exp10 spurious underflows (bug 16560).Joseph Myers2014-06-253-0/+6
| | | | | | | | | | | | | | | | | | | | | | | | | This patch fixes spurious underflows from exp10 for arguments near 0 (part of bug 16560; that bug also includes spurious underflows from exp2, which are not fixed by this patch). The problem is underflows in the internal computation converting the exp10 argument to arguments for exp (with extra precision), and the fix is simply to return 1 early for arguments near enough to 0 (just as arguments with large enough magnitude have their own overflow / underflow logic at the start of the function). Tested x86_64 and x86 and ulps updated accordingly; also tested for powerpc32 and mips64 to validate the ldbl-128ibm and ldbl-128 changes. [BZ #16560] * sysdeps/ieee754/dbl-64/e_exp10.c (__ieee754_exp10): Return 1 for arguments close to 0. * sysdeps/ieee754/ldbl-128/e_exp10l.c (__ieee754_exp10l): Likewise. * sysdeps/ieee754/ldbl-128ibm/e_exp10l.c (__ieee754_exp10l): Likewise. * math/auto-libm-test-in: Add more tests of exp10. * math/auto-libm-test-out: Regenerated. * sysdeps/x86_64/fpu/libm-test-ulps: Update.
* Fix ldbl-128 erfl spurious underflows (bug 16287).Joseph Myers2014-06-241-0/+2
| | | | | | | | | | | | | | | | | | | | This patch fixes bug 16287, spurious underflows from ldbl-128 erfl arising from it calling erfcl for arguments with absolute value at least 1.0, although for large positive arguments erfcl correctly underflows but erfl shouldn't. The fix is simply to avoid calling erfcl, and just return 1, for arguments above a cut-off large enough that erfl correctly rounds to-nearest as 1 but not so large that erfcl underflows. Tested mips64. Also tested x86_64 and x86 to confirm the new tests (taken from the tests of erfc) don't cause any problems there; no ulps updates needed. [BZ #16287] * sysdeps/ieee754/ldbl-128/s_erfl.c (__erfl): Return 1 without calling __erfcl for arguments at least 16. * math/auto-libm-test-in: Add more tests of erf. * math/auto-libm-test-out: Regenerated.
* Fix cosh spurious underflows from expm1 (bug 16354), inaccurate results near ↵Joseph Myers2014-06-235-7/+7
| | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | 0 (bug 17061). 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.
* Set errno for y1 overflow (bug 17050).Joseph Myers2014-06-234-4/+22
| | | | | | | | | | | | | | | | | | | | | | | | | 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.
* Fix pow overflow in non-default rounding modes (bug 16315).Joseph Myers2014-06-231-20/+41
| | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | This patch fixes bug 16315, bad pow handling of overflow/underflow in non-default rounding modes. Tests of pow are duly converted to ALL_RM_TEST to run all tests in all rounding modes. There are two main issues here. First, various implementations compute a negative result by negating a positive result, but this yields inappropriate overflow / underflow values for directed rounding, so either overflow / underflow results need recomputing in the correct sign, or the relevant overflowing / underflowing operation needs to be made to have a result of the correct sign. Second, the dbl-64 implementation sets FE_TONEAREST internally; in the overflow / underflow case, the result needs recomputing in the original rounding mode. Tested x86_64 and x86 and ulps updated accordingly. [BZ #16315] * sysdeps/i386/fpu/e_pow.S (__ieee754_pow): Ensure possibly overflowing or underflowing operations take place with sign of result. * sysdeps/i386/fpu/e_powf.S (__ieee754_powf): Likewise. * sysdeps/i386/fpu/e_powl.S (__ieee754_powl): Likewise. * sysdeps/ieee754/dbl-64/e_pow.c: Include <math.h>. (__ieee754_pow): Recompute overflowing and underflowing results in original rounding mode. * sysdeps/x86/fpu/powl_helper.c: Include <stdbool.h>. (__powl_helper): Allow negative argument X and scale negated value as needed. Avoid passing value outside [-1, 1] to f2xm1. * sysdeps/x86_64/fpu/e_powl.S (__ieee754_powl): Ensure possibly overflowing or underflowing operations take place with sign of result. * sysdeps/x86_64/fpu/multiarch/e_pow.c [HAVE_FMA4_SUPPORT]: Include <math.h>. * math/auto-libm-test-in: Add more tests of pow. * math/auto-libm-test-out: Regenerated. * math/libm-test.inc (pow_test): Use ALL_RM_TEST. (pow_tonearest_test_data): Remove. (pow_test_tonearest): Likewise. (pow_towardzero_test_data): Likewise. (pow_test_towardzero): Likewise. (pow_downward_test_data): Likewise. (pow_test_downward): Likewise. (pow_upward_test_data): Likewise. (pow_test_upward): Likewise. (main): Don't call removed functions. * sysdeps/i386/fpu/libm-test-ulps: Update. * sysdeps/x86_64/fpu/libm-test-ulps: Likewise.
* [BZ #6803] Set errno for scalbln, scalbnStefan Liebler2014-06-209-13/+50
| | | | | | | | | | | 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.
* PowerPC: Fix nearbyintl failure for few inputsRajalakshmi Srinivasaraghavan2014-06-171-0/+5
| | | | | | | | This patch fixes few failures in nearbyintl() where the fraction part is close to 0.5.i The new tests added report few extra failures in nearbyint_downward and nearbyint_towardzero which is a known issue. Fixes #17031.
* Fix log2 (1) in round-downward mode (bug 17042).Joseph Myers2014-06-102-0/+6
| | | | | | | | | | | | | | | | | | | | | | | | | | | As with other issues of this kind, bug 17042 is log2 (1) wrongly returning -0 instead of +0 in round-downward mode because of implementations effectively in terms of log1p (x - 1). This patch fixes the issue in the same way used for log and log10. 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 log2l is essentially the same as the ldbl-128 one). [BZ #17042] * sysdeps/i386/fpu/e_log2.S (__ieee754_log2): Take absolete value when x - 1 is zero. * sysdeps/i386/fpu/e_log2f.S (__ieee754_log2f): Likewise. * sysdeps/i386/fpu/e_log2l.S (__ieee754_log2l): Likewise. * sysdeps/ieee754/ldbl-128/e_log2l.c (__ieee754_log2l): Return 0.0L for an argument of 1.0L. * sysdeps/ieee754/ldbl-128ibm/e_log2l.c (__ieee754_log2l): Likewise. * sysdeps/x86_64/fpu/e_log2l.S (__ieee754_log2l): Take absolute value when x - 1 is zero. * math/libm-test.inc (log2_test): Use ALL_RM_TEST. * sysdeps/i386/fpu/libm-test-ulps: Update. * sysdeps/x86_64/fpu/libm-test-ulps: Likewise.
* Fix excessive ULP for y1_upward (0x2p+0) in test-float and test-ifloat.David S. Miller2014-05-311-0/+1
| | | | | * sysdeps/ieee754/flt-32/e_j1f.c (__ieee754_y1f): Force computations to occur in round to nearest mode when |x| >= 2.0
* Fix log10 (1) in round-downward mode (bug 16977).Joseph Myers2014-05-232-0/+6
| | | | | | | | | | | | | | | | | | | | | | | | | | | 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.
* Fix log1pl (LDBL_MAX) in FE_UPWARD mode (bug 16564).Joseph Myers2014-05-142-2/+8
| | | | | | | | | | | | | | | | | | | | | | | 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.
* Fix acosh (1) in round-downward mode (bug 16927).Joseph Myers2014-05-141-1/+1
| | | | | | | | | | | | | | | | | | | | | | | | | | | | | 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.
* Fix erf underflow handling near 0 (bug 16516).Joseph Myers2014-05-145-17/+54
| | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | 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.
* [BZ #16823] Fix log1pl returning wrong infinity signStefan Liebler2014-04-293-3/+3
|
* [BZ #16824] Fix failing y1 due to too large ulps in downward/upward rounding ↵Stefan Liebler2014-04-161-0/+1
| | | | mode.
* Correct IBM long double frexpl.Alan Modra2014-04-161-44/+102
| | | | | | | | | | | Besides fixing the bugzilla, this also fixes corner-cases where the high and low double differ greatly in magnitude, and handles a denormal input without resorting to a fp rescale. [BZ #16740] [BZ #16619] * sysdeps/ieee754/ldbl-128ibm/s_frexpl.c (__frexpl): Rewrite. * math/libm-test.inc (frexp_test_data): Add tests.
* Fix catan, catanh, __ieee754_logf in round-downward mode (bug 16799, bug 16800).Joseph Myers2014-04-021-1/+1
| | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | This patch fixes incorrect results from catan and catanh of certain special inputs in round-downward mode (bug 16799), and incorrect results of __ieee754_logf (+/-0) in round-downward mode (bug 16800) that show up through catan/catanh when tested in all rounding modes, but not directly in the testing for logf because the bug gets hidden by the wrappers. Both bugs involve a zero that should be +0 being -0 instead: one computed as (1-x)*(1+x) in the catan/catanh case, and one as (x-x) in the logf case. The fixes ensure positive zero is used. Testing of catan and catanh in all rounding modes is duly enabled. I expect there are various other bugs in special cases in __ieee754_* functions that are normally hidden by the wrappers but would show up for testing with -lieee (or in future with -fno-math-errno if we replace -lieee and _LIB_VERSION with compile-time redirection to new *_noerrno symbol names). Tested x86_64 and x86 and ulps updated accordingly. [BZ #16799] [BZ #16800] * math/s_catan.c (__catan): Avoid passing -0 denominator to atan2 with 0 numerator. * math/s_catanf.c (__catanf): Likewise. * math/s_catanh.c (__catanh): Likewise. * math/s_catanhf.c (__catanhf): Likewise. * math/s_catanhl.c (__catanhl): Likewise. * math/s_catanl.c (__catanl): Likewise. * sysdeps/ieee754/flt-32/e_logf.c (__ieee754_logf): Always divide by positive zero when computing -Inf result. * math/libm-test.inc (catan_test): Use ALL_RM_TEST. (catanh_test): Likewise. * sysdeps/i386/fpu/libm-test-ulps: Update. * sysdeps/x86_64/fpu/libm-test-ulps: Likewise.
* Correct IBM long double nextafterl.Alan Modra2014-04-021-16/+33
| | | | | | | | | | | | Fix for values near a power of two, and some tidies. [BZ #16739] * sysdeps/ieee754/ldbl-128ibm/s_nextafterl.c (__nextafterl): Correct output when value is near a power of two. Use int64_t for lx and remove casts. Use decimal rather than hex exponent constants. Don't use long double multiplication when double will suffice. * math/libm-test.inc (nextafter_test_data): Add tests. * NEWS: Add 16739 and 16786 to bug list.
* Fix implicit __isinf declarations in exp.Joseph Myers2014-03-241-0/+1
| | | | | | | | | | | | | | | My recent exp patch introduced warnings about implicit __isinf declarations in exp because e_exp.c didn't include <math.h>. This patch fixes this. Because <math.h> can't be included after <math_private.h> (because of macro definitions of __nan*), it was necessary to put an include in sysdeps/x86_64/fpu/multiarch/e_exp.c as well. Tested x86_64. * sysdeps/ieee754/dbl-64/e_exp.c: Include <math.h>. * sysdeps/x86_64/fpu/multiarch/e_exp.c [HAVE_FMA4_SUPPORT || HAVE_AVX_SUPPORT]: Likewise.
* Fix dbl-64 exp overflow/underflow in non-default rounding modes (bug 16284).Joseph Myers2014-03-241-143/+153
| | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | The dbl-64 version of exp needs round-to-nearest mode for its internal computations, but that has the consequence of inappropriate overflowing and underflowing results in other rounding modes. This patch fixes this by recomputing the relevant results in cases where the round-to-nearest result overflows to infinity or underflows to zero (most of the diffs are actually just consequent reindentation). Tests are enabled in all rounding modes for complex functions using exp - but not for cexp because it turns out there are bugs causing spurious underflows for cexp for some tests, which will need to be fixed separately (I suspect ccos ccosh csin csinh ctan ctanh have similar bugs, just not shown by the present set of test inputs). Tested x86_64 and x86 and ulps updated accordingly. [BZ #16284] * sysdeps/ieee754/dbl-64/e_exp.c (__ieee754_exp): Use original rounding mode to recompute results that overflow to infinity or underflow to zero. * math/auto-libm-test-in: Don't mark tests as expected to fail for bug 16284. * math/auto-libm-test-out: Regenerated. * math/libm-test.inc (ccos_test): Use ALL_RM_TEST. (ccosh_test): Likewise. (csin_test_data): Use plus_oflow. (csin_test): Use ALL_RM_TEST. (csinh_test_data): Use plus_oflow. (csinh_test): Use ALL_RM_TEST. * sysdeps/i386/fpu/libm-test-ulps: Update. * sysdeps/x86_64/fpu/libm-test-ulps: Likewise.
* Fix log (1) in round-downward mode (bug 16731).Joseph Myers2014-03-212-0/+6
| | | | | | | | | | | | | | | | | | | | | | | | | | | | | | 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.
* Fix nextafter overflow in non-default rounding modes (bug 16677).Joseph Myers2014-03-113-9/+16
| | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | 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.
* Fix sign of input to bsloww1 (BZ #16623)Siddhesh Poyarekar2014-02-271-6/+10
| | | | | | | | In 84ba214c, I removed some redundant sign computations and in the process, I incorrectly got rid of a temporary variable, thus passing the absolute value of the input to bsloww1. This caused #16623. This fix undoes the incorrect change.
* BZ #16447: Fix ldbl-128 expl implementation.Andreas Krebbel2014-02-111-1/+1
| | | | | Extend the range of numbers handled via unsafe mode. Add expl testcase and regenerate ULPs for s390.
* Use glibc_likely instead __builtin_expect.Ondřej Bílka2014-02-1042-129/+129
|
* [BZ #16427] Fix ldbl-128 exp overflows.Andreas Krebbel2014-01-151-13/+5
| | | | | Invoke the non-IEEE handling only for numbers special also in the IEEE case. This aligns the exp handling with the other ldbl variants.
* PowerPC: remove wrong truncl implementation for PowerPC64Adhemerval Zanella2014-01-081-3/+0
| | | | | | | | | | | | | | | The truncl assembly implementation (sysdeps/powerpc/powerpc64/fpu/s_truncl.S) returns wrong results for some inputs where first double is a exact integer and the precision is determined by second long double. Checking on implementation comments and history, I am very confident the assembly implementation was based on a version before commit 5c68d401698a58cf7da150d9cce769fa6679ba5f that fixes BZ#2423 (Errors in long double (ldbl-128ibm) rounding functions in glibc-2.4). By just removing the implementation and make the build select sysdeps/ieee754/ldbl-128ibm/s_truncl.c instead it fixes tgammal issues regarding wrong result sign.
* Fix ldbl-128ibm expm1l on large arguments (bug 16408).Joseph Myers2014-01-081-0/+2
| | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | This patch fixes bug 16408, ldbl-128ibm expm1l returning NaN for some large arguments. The basic problem is that the approach of converting the exponent to the form n * log(2) + y, where -0.5 <= y <= 0.5, then computing 2^n * expm1(y) + (2^n - 1) falls over when 2^n overflows (starting slightly before the point where expm1 overflows, when y is negative and n is the least integer for which 2^n overflows). The ldbl-128 code, and the x86/x86_64 code, make expm1l fall back to expl for large positive arguments to avoid this issue. This patch makes the ldbl-128ibm code do the same. (The problem appears for the particular argument in the testsuite because the ldbl-128ibm code also uses an overflow threshold that's for ldbl-128 and is too big for ldbl-128ibm, but the problem described applies for large non-overflowing cases as well, although during the freeze is not a suitable time for making the expm1 tests cover cases close to overflow more thoroughly.) This leaves some code for large positive arguments in expm1l that is now dead. To keep the code for ldbl-128 and ldbl-128ibm similar, and to avoid unnecessary changes during the freeze, the patch doesn't remove it; instead I propose to file a bug in Bugzilla as a reminder that this code (for overflow, including errno setting, and for arguments of +Inf) is no longer needed and should be removed from both those expm1l implementations. Tested powerpc32. * sysdeps/ieee754/ldbl-128ibm/s_expm1l.c (__expm1l): Use __expl for large positive arguments.
* Fix ldbl-128ibm coshl spurious overflows (bug 16407).Joseph Myers2014-01-071-1/+1
| | | | | | | | | | | | | | | | | | | | | This patch fixes bug 16407, spurious overflows from ldbl-128ibm coshl. The implementation assumed that a high part (reinterpreted as an integer) of the absolute value of the argument of 0x408633ce8fb9f87dLL or more meant overflow, but the actual threshold has high part 0x408633ce8fb9f87eLL (and a negative low part). The patch adjusts the threshold accordingly. sinhl probably has the same issue, but I didn't get that far in adding tests of special cases (such as just below and above overflow) before the freeze and during the freeze is not a suitable time to add them (as they'd require ulps to be regenerated again), so I'm not changing that function for now; when I add more tests of special cases, we'll discover whether sinhl indeed has this problem. Tested powerpc32. * sysdeps/ieee754/ldbl-128ibm/e_coshl.c (__ieee754_coshl): Increase overflow threshold.
* Fix ldbl-128 / ldbl-128ibm lgammal spurious underflow (bug 16400).Joseph Myers2014-01-061-1/+3
| | | | | | | | | | | | | | 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.
* Fix ldbl-128ibm logl inaccuracy (bug 16386).Joseph Myers2014-01-021-0/+8
| | | | | | | | | | | | This patch fixes bug 16386, ldbl-128ibm logl inaccuracy (with consequent inaccuracy for lgammal) for arguments where the high double is subnormal, which showed up while attempting to regenerate ulps for powerpc-nofpu for 2.19. The problem here is logic failing to allow for subnormals when calculating the exponent of the argument. Tested for powerpc-nofpu. * sysdeps/ieee754/ldbl-128ibm/e_logl.c (__ieee754_logl): Adjust numbers with subnormal high part when calculating exponent.
* Fix ldbl-128ibm asinhl inaccuracy (bug 16385).Joseph Myers2014-01-021-1/+1
| | | | | | | | | | | This patch fixes bug 16385, ldbl-128ibm asinhl inaccuracy, which showed up while attempting to regenerate ulps for powerpc-nofpu for 2.19. The problem here was use of fabs instead of fabsl meaning large arguments were reduced to the precision of double. Tested for powerpc-nofpu. * sysdeps/ieee754/ldbl-128ibm/s_asinhl.c (__asinhl): Use fabsl not fabs.
* Fix ldbl-128ibm acoshl inaccuracy (bug 16384).Joseph Myers2014-01-021-2/+2
| | | | | | | | | | | | | | | This patch fixes bug 16384, ldbl-128ibm acoshl inaccuracy, which showed up while attempting to regenerate ulps for powerpc-nofpu for 2.19. There were two separate problems, use of __log1p instead of __log1pl and an insufficiently accurate constant value for log 2 (which this patch replaces by use of M_LN2l), each of which could cause substantial inaccuracy in affected cases. Tested for powerpc-nofpu. * sysdeps/ieee754/ldbl-128ibm/e_acoshl.c (ln2): Initialize with M_LN2l. (__ieee754_acoshl): Use __log1pl not __log1p.
* Update copyright notices with scripts/update-copyrightsAllan McRae2014-01-01205-205/+205
|
* Fix typo in csloww()Allan McRae2013-12-271-1/+1
| | | | | An incorrect variable name was used during the refactoring done in commit 4aafb73c.
* Fix ldbl-128 lgammal for small negative arguments (bug 16337).Joseph Myers2013-12-221-2/+2
| | | | | | | | | | | | | | | | | | | | | 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.
* Consolidate code to compute sin and cos from lookup tablesSiddhesh Poyarekar2013-12-201-232/+129
| | | | | | | | | | | This patch consolidates the multiple copies of code that looks up sin and cos of a number from the lookup table and computes the final value, into static functions. This does not have a noticeable performance impact since the functions are inlined by gcc. There is further scope for consolidation in the functions but they cause a more noticable impact on performance (>5%) due to which I have held back on them.
* Remove more redundant computations in s_sin.cSiddhesh Poyarekar2013-12-201-42/+33
| | | | | | | | Removed more redundant computations in the slow paths of the sin and cos functions. The notable change is the passing of the most significant bits of X to the slow functions to check if X is positive so that just the absolute value of x can be passed and the repeated ABS() operation is avoided.
* Remove redundant arguments in reduce_and_computeSiddhesh Poyarekar2013-12-201-4/+4
| | | | | The A and DA arguments in reduce_and_compute are useless and hence have been removed.
* Remove some redundant computations in s_sin.cSiddhesh Poyarekar2013-12-201-24/+15
| | | | | | | | There are multiple points in the code where the absolute value of a number is computed multiple times or is computed even though the value can only be positive. This change removes those redundant computations. Tested on x86_64 to verify that there were no regressions in the testsuite.
* Compile e_sqrt.c with -ffp-contract=off.Marcus Shawcroft2013-12-181-0/+1
|
* Fix ldbl-128 logl for subnormals (bug 16338).Joseph Myers2013-12-181-3/+2
| | | | | | | | | | | | | | | | | | | | | | 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.
* Fix dbl-64 hypot spurious underflows (bug 16314).Joseph Myers2013-12-171-1/+1
|
* Fix hypot handling of subnormals (bug 16316, bug 16330).Joseph Myers2013-12-173-0/+33
|
* Minor code cleanup in s_sin.cSiddhesh Poyarekar2013-12-161-20/+13
| | | | | | | - Remove redundant mynumber union definitions - Clean up a clumsy ternary operator - Rename TAYLOR_SINCOS to TAYLOR_SIN since we're only expanding the sin Taylor series in it.
* Consolidate definition of constant t22Siddhesh Poyarekar2013-12-111-9/+3
|
* Use double constants instead of the struct numberSiddhesh Poyarekar2013-12-112-136/+106
|