summary refs log tree commit diff
path: root/sysdeps/ieee754/ldbl-96
Commit message (Collapse)AuthorAgeFilesLines
* Fix yn overflow handling in non-default rounding modes (bug 16561, bug 16562).Joseph Myers2014-06-271-63/+74
| | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | 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 cosh spurious underflows from expm1 (bug 16354), inaccurate results near ↵Joseph Myers2014-06-231-1/+1
| | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | 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-231-1/+5
| | | | | | | | | | | | | | | | | | | | | | | | | 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.
* [BZ #6803] Set errno for scalbln, scalbnStefan Liebler2014-06-201-1/+0
| | | | | | | | | | | 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.
* 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-141-3/+11
| | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | 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.
* Use glibc_likely instead __builtin_expect.Ondřej Bílka2014-02-106-27/+27
|
* Update copyright notices with scripts/update-copyrightsAllan McRae2014-01-0126-26/+26
|
* Fix hypot handling of subnormals (bug 16316, bug 16330).Joseph Myers2013-12-171-0/+11
|
* Fix Bessel function error handling (bug 6807, bug 15901).Joseph Myers2013-12-041-1/+2
|
* Fix erfc errno setting on underflow (bug 6786).Joseph Myers2013-12-031-2/+11
|
* Remove unused ldbl-96 functions (bug 15004).Joseph Myers2013-11-2816-1046/+0
|
* PowerPC floating point little-endian [1 of 15]Alan Modra2013-10-041-3/+5
| | | | | | | | | | | | | | | | | | | | | | | | | | 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.
* [BZ #15522] strtod ("nan(N)") returning a sNaN in some casesThomas Schwinge2013-08-291-5/+4
|
* Fix cbrtl for ldbl-96Andreas Schwab2013-08-131-1/+1
|
* Fix leading whitespaces.Ondrej Bilka2013-06-061-1/+1
|
* Remove trailing whitespace.Joseph Myers2013-06-051-2/+2
|
* Fix ldbl-96 hypotl of subnormals (bug 15529).Joseph Myers2013-05-241-1/+1
|
* Fix nearbyint scheduling of arithmetic past fesetenv (bug 15490).Joseph Myers2013-05-191-0/+2
|
* Improve tgamma accuracy (bugs 2546, 2560, 5159, 15426).Joseph Myers2013-05-083-11/+253
|
* New <math.h> macro named issignaling to check for a signaling NaN (sNaN).Thomas Schwinge2013-04-021-0/+43
| | | | It is based on draft TS 18661 and currently enabled as a GNU extension.
* Fix y1l spurious overflows for ldbl-96 (bug 15283).Joseph Myers2013-03-161-1/+1
|
* Update copyright notices with scripts/update-copyrights.Joseph Myers2013-01-0225-25/+25
|
* Use hex float 64-bit values in ldbl-96 asinl (bug 14803).Joseph Myers2012-11-281-3/+6
|
* Fix fma underflows with small x * y (bug 14793).Joseph Myers2012-11-061-18/+27
|
* Fix fma overflow results outside round-to-nearest mode (bug 14797).Joseph Myers2012-11-041-4/+6
|
* Make fma use of Dekker and Knuth algorithms use round-to-nearest (bug 14796).Joseph Myers2012-11-032-4/+32
|
* Fix fma (a, b, c) for small a * b (bugs 14784, 14785).Joseph Myers2012-11-011-4/+32
|
* Fix fma underflow exceptions in after-rounding edge cases.Joseph Myers2012-10-311-0/+12
|
* Fix fma missing underflows and bad results for some subnormal results (bugs ↵Joseph Myers2012-10-301-14/+8
| | | | 14152, 14783).
* Fix sign of inexact zero return from fma (bug 14645).Joseph Myers2012-10-011-0/+5
|
* Fix sign of exact zero return from fma (bug 14638).Joseph Myers2012-09-292-1/+10
|
* Fix inaccuracy of clog, clog10 near |z| = 1 (bug 13629).Joseph Myers2012-09-252-0/+150
|
* Quash warning in s_sincosl.Marek Polacek2012-08-171-2/+2
|
* Set up errno properly for yn.Marek Polacek2012-07-251-0/+4
|
* Fix ynl return value with LDBL_MIN.Marek Polacek2012-07-121-1/+2
|
* Ensure additions are not scheduled after fetestexcept in fmaf and fmal.Joseph Myers2012-06-011-1/+4
|
* Optimize handling of denormals in logb/logbf/logblAndreas Schwab2012-05-261-4/+6
|
* Fix for logb/logbf/logbl (bugs 13954/13955/13956)Adhemerval Zanella2012-05-101-14/+21
| | | | | | | | 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.
* Fix missing exceptions from exp (bugs 13787, 13922, 14036).Joseph Myers2012-05-051-18/+5
|
* Fix nexttoward bugs (bugs 2550, 2570).Joseph Myers2012-05-012-16/+4
|
* Finish ilogb changesAndreas Jaeger2012-04-181-2/+1
| | | | | | | | | | | | | | | | | | [BZ# 6794] * sysdeps/ieee754/ldbl-96/s_ilogbl.c: Moved to ... * sysdeps/ieee754/ldbl-96/e_ilogbl.c: ... here. Rename __ilogbl to __ieee754_ilogbl and remove weak_alias. * sysdeps/ieee754/ldbl-128/s_ilogbl.c: Moved to ... * sysdeps/ieee754/ldbl-128/e_ilogbl.c: ... here. Rename __ilogbl to __ieee754_ilogbl and remove weak_alias. * sysdeps/ieee754/ldbl-64-128/s_ilogbl.c: Moved to ... * sysdeps/ieee754/ldbl-64-128/e_ilogbl.c: ... here. * sysdeps/sparc/sparc64/soft-fp/s_ilogbl.c: Moved to ... * sysdeps/sparc/sparc64/soft-fp/e_ilogbl.c: ... here. Rename __ilogbl to __ieee754_ilogbl and remove weak_alias.
* Avoid overflows from long double functions using __kernel_standard.Joseph Myers2012-03-281-2/+2
|
* Fix tan, tanl for large inputs.Joseph Myers2012-03-162-3/+145
|
* Include program generating __sincosl_table in comment.Joseph Myers2012-03-161-2/+41
|
* Implement ldbl-96 sinl / cosl / sincosl (bug 13851).Joseph Myers2012-03-166-2/+931
|
* Use <> for math.h and math_private.h everywhere.Richard Henderson2012-03-0948-83/+83
| | | | Entire tree edited via find | grep | sed.
* Fix up long double fphex.Marek Polacek2012-03-061-1/+2
|
* Fix scalbn, scalbln integer overflow.Joseph Myers2012-03-022-6/+10
|
* Use non-signaling floating-point comparisons in math functions.Aurelien Jarno2012-02-191-3/+3
|