aboutsummaryrefslogtreecommitdiff
path: root/sysdeps/ieee754/flt-32/e_j0f.c
AgeCommit message (Collapse)Author
2022-01-10math: Fix float conversion regressions with gcc-12 [BZ #28713]Szabolcs Nagy
Converting double precision constants to float is now affected by the runtime dynamic rounding mode instead of being evaluated at compile time with default rounding mode (except static object initializers). This can change the computed result and cause performance regression. The known correctness issues (increased ulp errors) are already fixed, this patch fixes remaining cases of unnecessary runtime conversions. Add float M_* macros to math.h as new GNU extension API. To avoid conversions the new M_* macros are used and instead of casting double literals to float, use float literals (only required if the conversion is inexact). The patch was tested on aarch64 where the following symbols had new spurious conversion instructions that got fixed: __clog10f __gammaf_r_finite@GLIBC_2.17 __j0f_finite@GLIBC_2.17 __j1f_finite@GLIBC_2.17 __jnf_finite@GLIBC_2.17 __kernel_casinhf __lgamma_negf __log1pf __y0f_finite@GLIBC_2.17 __y1f_finite@GLIBC_2.17 cacosf cacoshf casinhf catanf catanhf clogf gammaf_positive Fixes bug 28713. Reviewed-by: Paul Zimmermann <Paul.Zimmermann@inria.fr>
2021-10-05Fixed inaccuracy of j0f (BZ #28185)Paul Zimmermann
The largest errors over the full binary32 range are after this patch (on x86_64): RNDN: libm wrong by up to 9.00e+00 ulp(s) [9] for x=0x1.04c39cp+6 RNDZ: libm wrong by up to 9.00e+00 ulp(s) [9] for x=0x1.04c39cp+6 RNDU: libm wrong by up to 9.00e+00 ulp(s) [9] for x=0x1.04c39cp+6 RNDD: libm wrong by up to 8.98e+00 ulp(s) [9] for x=0x1.4b7066p+7 Inputs that were yielding huge errors have been added to "make check". Reviewed-by: Adhemeral Zanella <adhemerval.zanella@linaro.org>
2021-09-03Remove "Contributed by" linesSiddhesh Poyarekar
We stopped adding "Contributed by" or similar lines in sources in 2012 in favour of git logs and keeping the Contributors section of the glibc manual up to date. Removing these lines makes the license header a bit more consistent across files and also removes the possibility of error in attribution when license blocks or files are copied across since the contributed-by lines don't actually reflect reality in those cases. Move all "Contributed by" and similar lines (Written by, Test by, etc.) into a new file CONTRIBUTED-BY to retain record of these contributions. These contributors are also mentioned in manual/contrib.texi, so we just maintain this additional record as a courtesy to the earlier developers. The following scripts were used to filter a list of files to edit in place and to clean up the CONTRIBUTED-BY file respectively. These were not added to the glibc sources because they're not expected to be of any use in future given that this is a one time task: https://gist.github.com/siddhesh/b5ecac94eabfd72ed2916d6d8157e7dc https://gist.github.com/siddhesh/15ea1f5e435ace9774f485030695ee02 Reviewed-by: Carlos O'Donell <carlos@redhat.com>
2021-04-02Fix the inaccuracy of j0f/j1f/y0f/y1f [BZ #14469, #14470, #14471, #14472]Paul Zimmermann
For j0f/j1f/y0f/y1f, the largest error for all binary32 inputs is reduced to at most 9 ulps for all rounding modes. The new code is enabled only when there is a cancellation at the very end of the j0f/j1f/y0f/y1f computation, or for very large inputs, thus should not give any visible slowdown on average. Two different algorithms are used: * around the first 64 zeros of j0/j1/y0/y1, approximation polynomials of degree 3 are used, computed using the Sollya tool (https://www.sollya.org/) * for large inputs, an asymptotic formula from [1] is used [1] Fast and Accurate Bessel Function Computation, John Harrison, Proceedings of Arith 19, 2009. Inputs yielding the new largest errors are added to auto-libm-test-in, and ulps are regenerated for various targets (thanks Adhemerval Zanella). Tested on x86_64 with --disable-multi-arch and on powerpc64le-linux-gnu. Reviewed-by: Adhemerval Zanella <adhemerval.zanella@linaro.org>
2020-08-07math: Fix inaccuracy of j0f for x >= 2^127 when sin(x)+cos(x) is tinyPaul Zimmermann
Checked on x86_64-linux-gnu and i686-linux-gnu.
2020-02-14Adjust thresholds in Bessel function implementations (bug 14469).Joseph Myers
A recent discussion in bug 14469 notes that a threshold in float Bessel function implementations, used to determine when to use a simpler implementation approach, results in substantially inaccurate results. As I discussed in <https://sourceware.org/ml/libc-alpha/2013-03/msg00345.html>, a heuristic argument suggests 2^(S+P) as the right order of magnitude for a suitable threshold, where S is the number of significand bits in the floating-point type and P is the number of significant bits in the representation of the floating-point type, and the float and ldbl-96 implementations use thresholds that are too small. Some threshold does need using, there or elsewhere in the implementation, to avoid spurious underflow and overflow for large arguments. This patch sets the thresholds in the affected implementations to more heuristically justifiable values. Results will still be inaccurate close to zeroes of the functions (thus this patch does *not* fix any of the bugs for Bessel function inaccuracy); fixing that would require a different implementation approach, likely along the lines described in <http://www.cl.cam.ac.uk/~jrh13/papers/bessel.ps.gz>. So the justification for a change such as this would be statistical rather than based on particular tests that had excessive errors and no longer do so (no doubt such tests could be found, but would probably be too fragile to add to the testsuite, as liable to give large errors again from very small implementation changes or even from compiler changes). See <https://sourceware.org/ml/libc-alpha/2020-02/msg00638.html> for such statistics of the resulting improvements for float functions. Tested (glibc testsuite) for x86_64.
2020-01-03Add libm_alias_finite for _finite symbolsWilco Dijkstra
This patch adds a new macro, libm_alias_finite, to define all _finite symbol. It sets all _finite symbol as compat symbol based on its first version (obtained from the definition at built generated first-versions.h). The <fn>f128_finite symbols were introduced in GLIBC 2.26 and so need special treatment in code that is shared between long double and float128. It is done by adding a list, similar to internal symbol redifinition, on sysdeps/ieee754/float128/float128_private.h. Alpha also needs some tricky changes to ensure we still emit 2 compat symbols for sqrt(f). Passes buildmanyglibc. Co-authored-by: Adhemerval Zanella <adhemerval.zanella@linaro.org> Reviewed-by: Siddhesh Poyarekar <siddhesh@sourceware.org>
2018-05-11Do not include math-barriers.h in math_private.h.Joseph Myers
This patch continues the math_private.h cleanup by stopping math_private.h from including math-barriers.h and making the users of the barrier macros include the latter header directly. No attempt is made to remove any math_private.h includes that are now unused, except in strtod_l.c where that is done to avoid line number changes in assertions, so that installed stripped shared libraries can be compared before and after the patch. (I think the floating-point environment support in math_private.h should also move out - some architectures already have fenv_private.h as an architecture-internal header included from their math_private.h - and after moving that out might be a better time to identify unused math_private.h includes.) Tested for x86_64 and x86, and tested with build-many-glibcs.py that installed stripped shared libraries are unchanged by the patch. * sysdeps/generic/math_private.h: Do not include <math-barriers.h>. * stdlib/strtod_l.c: Include <math-barriers.h> instead of <math_private.h>. * math/fromfp.h: Include <math-barriers.h>. * math/math-narrow.h: Likewise. * math/s_nextafter.c: Likewise. * math/s_nexttowardf.c: Likewise. * sysdeps/aarch64/fpu/s_llrint.c: Likewise. * sysdeps/aarch64/fpu/s_llrintf.c: Likewise. * sysdeps/aarch64/fpu/s_lrint.c: Likewise. * sysdeps/aarch64/fpu/s_lrintf.c: Likewise. * sysdeps/i386/fpu/s_nextafterl.c: Likewise. * sysdeps/i386/fpu/s_nexttoward.c: Likewise. * sysdeps/i386/fpu/s_nexttowardf.c: Likewise. * sysdeps/ieee754/dbl-64/e_atan2.c: Likewise. * sysdeps/ieee754/dbl-64/e_atanh.c: Likewise. * sysdeps/ieee754/dbl-64/e_exp.c: Likewise. * sysdeps/ieee754/dbl-64/e_exp2.c: Likewise. * sysdeps/ieee754/dbl-64/e_j0.c: Likewise. * sysdeps/ieee754/dbl-64/e_sqrt.c: Likewise. * sysdeps/ieee754/dbl-64/s_expm1.c: Likewise. * sysdeps/ieee754/dbl-64/s_fma.c: Likewise. * sysdeps/ieee754/dbl-64/s_fmaf.c: Likewise. * sysdeps/ieee754/dbl-64/s_log1p.c: Likewise. * sysdeps/ieee754/dbl-64/s_nearbyint.c: Likewise. * sysdeps/ieee754/dbl-64/wordsize-64/s_nearbyint.c: Likewise. * sysdeps/ieee754/flt-32/e_atanhf.c: Likewise. * sysdeps/ieee754/flt-32/e_j0f.c: Likewise. * sysdeps/ieee754/flt-32/s_expm1f.c: Likewise. * sysdeps/ieee754/flt-32/s_log1pf.c: Likewise. * sysdeps/ieee754/flt-32/s_nearbyintf.c: Likewise. * sysdeps/ieee754/flt-32/s_nextafterf.c: Likewise. * sysdeps/ieee754/k_standardl.c: Likewise. * sysdeps/ieee754/ldbl-128/e_asinl.c: Likewise. * sysdeps/ieee754/ldbl-128/e_expl.c: Likewise. * sysdeps/ieee754/ldbl-128/e_powl.c: Likewise. * sysdeps/ieee754/ldbl-128/s_fmal.c: Likewise. * sysdeps/ieee754/ldbl-128/s_nearbyintl.c: Likewise. * sysdeps/ieee754/ldbl-128/s_nextafterl.c: Likewise. * sysdeps/ieee754/ldbl-128/s_nexttoward.c: Likewise. * sysdeps/ieee754/ldbl-128/s_nexttowardf.c: Likewise. * sysdeps/ieee754/ldbl-128ibm/e_asinl.c: Likewise. * sysdeps/ieee754/ldbl-128ibm/s_fmal.c: Likewise. * sysdeps/ieee754/ldbl-128ibm/s_nextafterl.c: Likewise. * sysdeps/ieee754/ldbl-128ibm/s_nexttoward.c: Likewise. * sysdeps/ieee754/ldbl-128ibm/s_nexttowardf.c: Likewise. * sysdeps/ieee754/ldbl-128ibm/s_rintl.c: Likewise. * sysdeps/ieee754/ldbl-96/e_atanhl.c: Likewise. * sysdeps/ieee754/ldbl-96/e_j0l.c: Likewise. * sysdeps/ieee754/ldbl-96/s_fma.c: Likewise. * sysdeps/ieee754/ldbl-96/s_fmal.c: Likewise. * sysdeps/ieee754/ldbl-96/s_nexttoward.c: Likewise. * sysdeps/ieee754/ldbl-96/s_nexttowardf.c: Likewise. * sysdeps/ieee754/ldbl-opt/s_nexttowardfd.c: Likewise. * sysdeps/m68k/m680x0/fpu/s_nextafterl.c: Likewise.
2018-03-15Rename all __ieee754_sqrt(f/l) calls to sqrt(f/l)Wilco Dijkstra
Use sqrt(f/l) to enable inlining by GCC - if inlining doesn't happen, the asm redirect ensures we will still call __ieee754_sqrt(f/l). * sysdeps/ieee754/dbl-64/e_acosh.c (__ieee754_acosh): Use sqrt. * sysdeps/ieee754/dbl-64/e_gamma_r.c (gamma_positive): Likewise. * sysdeps/ieee754/dbl-64/e_hypot.c (__ieee754_hypot): Likewise. * sysdeps/ieee754/dbl-64/e_j0.c (__ieee754_j0): Likewise. * sysdeps/ieee754/dbl-64/e_j1.c (__ieee754_j1): Likewise. * sysdeps/ieee754/dbl-64/e_jn.c (__ieee754_jn): Likewise. * sysdeps/ieee754/dbl-64/s_asinh.c (__asinh): Likewise. * sysdeps/ieee754/dbl-64/wordsize-64/e_acosh.c (__ieee754_acosh): Likewise. * sysdeps/ieee754/flt-32/e_acosf.c (__ieee754_acosf): Likewise. * sysdeps/ieee754/flt-32/e_acoshf.c (__ieee754_acoshf): Likewise. * sysdeps/ieee754/flt-32/e_asinf.c (__ieee754_asinf): Likewise. * sysdeps/ieee754/flt-32/e_gammaf_r.c (gammaf_positive): Likewise. * sysdeps/ieee754/flt-32/e_hypotf.c (__ieee754_hypotf): Likewise. * sysdeps/ieee754/flt-32/e_j0f.c (__ieee754_j0f): Likewise. * sysdeps/ieee754/flt-32/e_j1f.c (__ieee754_j1f): Likewise. * sysdeps/ieee754/flt-32/e_powf.c (__ieee754_powf): Likewise. * sysdeps/ieee754/flt-32/s_asinhf.c (__asinhf): Likewise. * sysdeps/ieee754/ldbl-128/e_acoshl.c (__ieee754_acoshl): Use sqrtl. * sysdeps/ieee754/ldbl-128/e_acosl.c (__ieee754_acosl): Likewise. * sysdeps/ieee754/ldbl-128/e_asinl.c (__ieee754_asinl): Likewise. * sysdeps/ieee754/ldbl-128/e_gammal_r.c (gammal_positive): Likewise. * sysdeps/ieee754/ldbl-128/e_hypotl.c (__ieee754_hypotl): Likewise. * sysdeps/ieee754/ldbl-128/e_j0l.c (__ieee754_j0l): Likewise. * sysdeps/ieee754/ldbl-128/e_j1l.c (__ieee754_j1l): Likewise. * sysdeps/ieee754/ldbl-128/e_jnl.c (__ieee754_jnl): Likewise. * sysdeps/ieee754/ldbl-128/e_powl.c (__ieee754_powl): Likewise. * sysdeps/ieee754/ldbl-128/s_asinhl.c (__ieee754_asinhl): Likewise. * sysdeps/ieee754/ldbl-128ibm/e_acoshl.c (__ieee754_acoshl): Likewise. * sysdeps/ieee754/ldbl-128ibm/e_acosl.c (__ieee754_acosl): Likewise. * sysdeps/ieee754/ldbl-128ibm/e_asinl.c (__ieee754_asinl): Likewise. * sysdeps/ieee754/ldbl-128ibm/e_gammal_r.c (gammal_positive): Likewise. * sysdeps/ieee754/ldbl-128ibm/e_hypotl.c (__ieee754_hypotl): Likewise. * sysdeps/ieee754/ldbl-128ibm/e_j0l.c (__ieee754_j0l): Likewise. * sysdeps/ieee754/ldbl-128ibm/e_j1l.c (__ieee754_j1l): Likewise * sysdeps/ieee754/ldbl-128ibm/e_jnl.c (__ieee754_jnl): Likewise. * sysdeps/ieee754/ldbl-128ibm/e_powl.c (__ieee754_powl): Likewise. * sysdeps/ieee754/ldbl-128ibm/s_asinhl.c (__ieee754_asinhl): Likewise. * sysdeps/ieee754/ldbl-96/e_acoshl.c (__ieee754_acoshl): Use sqrtl. * sysdeps/ieee754/ldbl-96/e_asinl.c (__ieee754_asinl): Likewise. * sysdeps/ieee754/ldbl-96/e_gammal_r.c (gammal_positive): Likewise. * sysdeps/ieee754/ldbl-96/e_hypotl.c (__ieee754_hypotl): Likewise. * sysdeps/ieee754/ldbl-96/e_j0l.c (__ieee754_j0l): Likewise. * sysdeps/ieee754/ldbl-96/e_j1l.c (__ieee754_j1l): Likewise. * sysdeps/ieee754/ldbl-96/e_jnl.c (__ieee754_jnl): Likewise. * sysdeps/ieee754/ldbl-96/s_asinhl.c (__ieee754_asinhl): Likewise. * sysdeps/m68k/m680x0/fpu/e_pow.c (__ieee754_pow): Likewise. * sysdeps/powerpc/fpu/e_hypot.c (__ieee754_hypot): Likewise. * sysdeps/powerpc/fpu/e_hypotf.c (__ieee754_hypotf): Likewise.
2017-02-15Fix y0 and y1 exception handling for zero input [BZ #21134]Gabriel F. T. Gomes
The Bessel functions of the second type (Yn) should raise the "divide by zero" exception when input is zero (both positive and negative). Current code gives the right output, but fails to set the exception. This error is exposed for float, double, and long double when linking with -lieee. Without this flag, the error is not exposed, because the wrappers for these functions, which use __kernel_standard functionality, set the exception as expected. Tested for powerpc64le. [BZ #21134] * sysdeps/ieee754/dbl-64/e_j0.c (__ieee754_y0): Raise the "divide by zero" exception when the input is zero. * sysdeps/ieee754/dbl-64/e_j1.c (__ieee754_y1): Likewise. * sysdeps/ieee754/flt-32/e_j0f.c (__ieee754_y0f): Likewise. * sysdeps/ieee754/flt-32/e_j1f.c (__ieee754_y1f): Likewise. * sysdeps/ieee754/ldbl-128/e_j0l.c (__ieee754_y0l): Likewise. * sysdeps/ieee754/ldbl-128/e_j1l.c (__ieee754_y1l): Likewise.
2015-02-26Avoid uninitialized warnings in Bessel functions.Joseph Myers
math/Makefile currently has: # The fdlibm code generates a lot of these warnings but is otherwise clean. override CFLAGS += -Wno-uninitialized This is of course undesirable; warnings should be disabled as narrowly as possible. To remove this override, we need to fix files that generate such warnings, or put warning-disabling pragmas in them. This patch does so for Bessel function implementations, one of the cases that have the warnings if the override is removed. The warnings arise because functions set pointer variables p and q only for certain values of the function argument, then use them unconditionally. As the static functions in question only get called for arguments that satisfy the last condition in the if/else chain, the natural fix is to change the last "else if" to just "else", which this patch does. (The ldbl-128 / ldbl-128ibm implementation of these functions is substantially different and looks like it already does use "else" in the last case in the nearest corresponding code.) Tested for x86_64 and x86. * sysdeps/ieee754/dbl-64/e_j0.c (pzero): Change last case for setting p and q from "else if" to "else". (qzero): Likewise. * sysdeps/ieee754/dbl-64/e_j1.c (pone): Likewise. (qone): Likewise. * sysdeps/ieee754/flt-32/e_j0f.c (pzerof): Likewise. (qzerof): Likewise. * sysdeps/ieee754/flt-32/e_j1f.c (ponef): Likewise. (qonef): Likewise. * sysdeps/ieee754/ldbl-96/e_j0l.c (pzero): Likewise. (qzero): Likewise. * sysdeps/ieee754/ldbl-96/e_j1l.c (pone): Likewise. (qone): Likewise.
2012-11-18Correct tinyness handling in long-double and float y0/y1.David S. Miller
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.
2012-03-09Use <> for math.h and math_private.h everywhere.Richard Henderson
Entire tree edited via find | grep | sed.
2011-10-26Fix uses of math_force_evalAndreas Schwab
2011-10-25Use math_force_eval in more placesUlrich Drepper
2011-10-12Optimize libmUlrich Drepper
libm is now somewhat integrated with gcc's -ffinite-math-only option and lots of the wrapper functions have been optimized.
2003-12-28Update.Andreas Jaeger
* sysdeps/ieee754/dbl-64/e_j0.c (__ieee754_y0): Raise only overflow for 0 as argument. Raise Invalid exception for negative args. * sysdeps/ieee754/dbl-64/e_jn.c (__ieee754_yn): Likewise. * sysdeps/ieee754/dbl-64/e_j1.c (__ieee754_y0): Likewise. * sysdeps/ieee754/ldb-128/e_jnl.c (__ieee754_ynl): Likewise. * sysdeps/ieee754/ldb-128/e_j0l.c (__ieee754_y0l): Likewise. * sysdeps/ieee754/ldb-128/e_j1l.c (__ieee754_y1l): Likewise. * sysdeps/ieee754/ldb-96/e_jnl.c (__ieee754_ynl): Likewise. * sysdeps/ieee754/ldb-96/e_j0l.c (__ieee754_y0l): Likewise. * sysdeps/ieee754/ldb-96/e_j1l.c (__ieee754_y1l): Likewise. * sysdeps/ieee754/flt-32/e_jnf.c (__ieee754_ynf): Likewise. * sysdeps/ieee754/flt-32/e_j0f.c (__ieee754_y0f): Likewise. * sysdeps/ieee754/flt-32/e_j1f.c (__ieee754_y1f): Likewise. * math/libm-test.inc (yn_test): Expect invalid exception for negative arguments. (y0_test): Likewise. (y1_test): Likewise.
2001-02-13Little optimization, use sincos.Ulrich Drepper
2001-02-13Update.Ulrich Drepper
2001-02-12 Ulrich Drepper <drepper@redhat.com> * sysdeps/dbl-64/e_j0.c: Little optimization, use sincos. * sysdeps/flt-32/e_j0f.c: Likewise. * sysdeps/ldbl-96/e_j0l.c: Likewise.
2000-10-25Update.Ulrich Drepper
2000-10-25 Ulrich Drepper <drepper@redhat.com> * sysdeps/ieee754/dbl-64/e_jn.c: Use __ieee754_sqrt instead of __sqrt. * sysdeps/ieee754/dbl-64/e_j1.c: Likewise. * sysdeps/ieee754/dbl-64/e_j0.c: Likewise. * sysdeps/ieee754/flt-32/e_j1f.c: Likewise. * sysdeps/ieee754/flt-32/e_j0f.c: Likewise.
1999-07-14Update.Ulrich Drepper