https://gcc.gnu.org/g:62bc11d461f6e5a150c6c0593a0cec5851b76b19
commit r17-2310-g62bc11d461f6e5a150c6c0593a0cec5851b76b19 Author: Jonathan Wakely <[email protected]> Date: Fri Jun 5 15:51:23 2026 +0100 libstdc++: Fix Cylindrical Bessel functions for infinity [PR119143] The result should be zero at infinity. libstdc++-v3/ChangeLog: PR libstdc++/119143 * include/tr1/bessel_function.tcc (__cyl_bessel_jn_asymp): Return zero for infinite x. * testsuite/special_functions/08_cyl_bessel_j/check_inf.cc: New test. * testsuite/special_functions/10_cyl_neumann/check_inf.cc: New test. * testsuite/tr1/5_numerical_facilities/special_functions/11_cyl_neumann/check_inf.cc: New test. * testsuite/tr1/5_numerical_facilities/special_functions/09_cyl_bessel_j/check_inf.cc: New test. Reviewed-by: Tomasz Kamiński <[email protected]> Diff: --- libstdc++-v3/include/tr1/bessel_function.tcc | 6 ++++ .../special_functions/08_cyl_bessel_j/check_inf.cc | 38 ++++++++++++++++++++++ .../special_functions/10_cyl_neumann/check_inf.cc | 38 ++++++++++++++++++++++ .../special_functions/09_cyl_bessel_j/check_inf.cc | 36 ++++++++++++++++++++ .../special_functions/11_cyl_neumann/check_inf.cc | 37 +++++++++++++++++++++ 5 files changed, 155 insertions(+) diff --git a/libstdc++-v3/include/tr1/bessel_function.tcc b/libstdc++-v3/include/tr1/bessel_function.tcc index 7e06d98ce3c6..bba792bd89d3 100644 --- a/libstdc++-v3/include/tr1/bessel_function.tcc +++ b/libstdc++-v3/include/tr1/bessel_function.tcc @@ -362,6 +362,12 @@ namespace tr1 void __cyl_bessel_jn_asymp(_Tp __nu, _Tp __x, _Tp & __Jnu, _Tp & __Nnu) { + if (__builtin_expect(__builtin_isinf(__x), 0)) + { + __Jnu = __Nnu = _Tp(0); + return; + } + const _Tp __mu = _Tp(4) * __nu * __nu; const _Tp __8x = _Tp(8) * __x; diff --git a/libstdc++-v3/testsuite/special_functions/08_cyl_bessel_j/check_inf.cc b/libstdc++-v3/testsuite/special_functions/08_cyl_bessel_j/check_inf.cc new file mode 100644 index 000000000000..50c1f3098297 --- /dev/null +++ b/libstdc++-v3/testsuite/special_functions/08_cyl_bessel_j/check_inf.cc @@ -0,0 +1,38 @@ +// { dg-do run { target c++11 } } +// { dg-require-c-std "" } +// { dg-options "-D__STDCPP_WANT_MATH_SPEC_FUNCS__" } +// { dg-add-options ieee } + +// ISO/IEC IS 29124 8.1.8 cyl_bessel_j +// J_nu(+inf) == 0 for all finite nu >= 0 + +#include <cmath> +#include <limits> +#include <testsuite_hooks.h> + +void +test01() +{ + float xinf_f = std::numeric_limits<float>::infinity(); + double xinf_d = std::numeric_limits<double>::infinity(); + long double xinf_l = std::numeric_limits<long double>::infinity(); + + VERIFY(std::cyl_bessel_j(0.0F, xinf_f) == 0.0F); + VERIFY(std::cyl_bessel_jf(0.0F, xinf_f) == 0.0F); + VERIFY(std::cyl_bessel_j(0.0, xinf_d) == 0.0); + VERIFY(std::cyl_bessel_j(0.0L, xinf_l) == 0.0L); + VERIFY(std::cyl_bessel_jl(0.0L, xinf_l) == 0.0L); + + VERIFY(std::cyl_bessel_j(1.0F, xinf_f) == 0.0F); + VERIFY(std::cyl_bessel_j(1.0, xinf_d) == 0.0); + VERIFY(std::cyl_bessel_j(1.0L, xinf_l) == 0.0L); + + VERIFY(std::cyl_bessel_j(2.5, xinf_d) == 0.0); +} + +int +main() +{ + test01(); + return 0; +} diff --git a/libstdc++-v3/testsuite/special_functions/10_cyl_neumann/check_inf.cc b/libstdc++-v3/testsuite/special_functions/10_cyl_neumann/check_inf.cc new file mode 100644 index 000000000000..716dc79b1aab --- /dev/null +++ b/libstdc++-v3/testsuite/special_functions/10_cyl_neumann/check_inf.cc @@ -0,0 +1,38 @@ +// { dg-do run { target c++11 } } +// { dg-require-c-std "" } +// { dg-options "-D__STDCPP_WANT_MATH_SPEC_FUNCS__" } +// { dg-add-options ieee } + +// ISO/IEC IS 29124 8.1.10 cyl_neumann +// N_nu(+inf) == 0 for all finite nu >= 0 + +#include <cmath> +#include <limits> +#include <testsuite_hooks.h> + +void +test01() +{ + float xinf_f = std::numeric_limits<float>::infinity(); + double xinf_d = std::numeric_limits<double>::infinity(); + long double xinf_l = std::numeric_limits<long double>::infinity(); + + VERIFY(std::cyl_neumann(0.0F, xinf_f) == 0.0F); + VERIFY(std::cyl_neumannf(0.0F, xinf_f) == 0.0F); + VERIFY(std::cyl_neumann(0.0, xinf_d) == 0.0); + VERIFY(std::cyl_neumann(0.0L, xinf_l) == 0.0L); + VERIFY(std::cyl_neumannl(0.0L, xinf_l) == 0.0L); + + VERIFY(std::cyl_neumann(1.0F, xinf_f) == 0.0F); + VERIFY(std::cyl_neumann(1.0, xinf_d) == 0.0); + VERIFY(std::cyl_neumann(1.0L, xinf_l) == 0.0L); + + VERIFY(std::cyl_neumann(2.5, xinf_d) == 0.0); +} + +int +main() +{ + test01(); + return 0; +} diff --git a/libstdc++-v3/testsuite/tr1/5_numerical_facilities/special_functions/09_cyl_bessel_j/check_inf.cc b/libstdc++-v3/testsuite/tr1/5_numerical_facilities/special_functions/09_cyl_bessel_j/check_inf.cc new file mode 100644 index 000000000000..d05b7b8c9e3b --- /dev/null +++ b/libstdc++-v3/testsuite/tr1/5_numerical_facilities/special_functions/09_cyl_bessel_j/check_inf.cc @@ -0,0 +1,36 @@ +// { dg-require-c-std "" } +// { dg-add-options ieee } + +// TR1 5.2.1.9 cyl_bessel_j +// J_nu(+inf) == 0 for all finite nu >= 0 + +#include <tr1/cmath> +#include <limits> +#include <testsuite_hooks.h> + +void +test01() +{ + float xinf_f = std::numeric_limits<float>::infinity(); + double xinf_d = std::numeric_limits<double>::infinity(); + long double xinf_l = std::numeric_limits<long double>::infinity(); + + VERIFY(std::tr1::cyl_bessel_j(0.0F, xinf_f) == 0.0F); + VERIFY(std::tr1::cyl_bessel_jf(0.0F, xinf_f) == 0.0F); + VERIFY(std::tr1::cyl_bessel_j(0.0, xinf_d) == 0.0); + VERIFY(std::tr1::cyl_bessel_j(0.0L, xinf_l) == 0.0L); + VERIFY(std::tr1::cyl_bessel_jl(0.0L, xinf_l) == 0.0L); + + VERIFY(std::tr1::cyl_bessel_j(1.0F, xinf_f) == 0.0F); + VERIFY(std::tr1::cyl_bessel_j(1.0, xinf_d) == 0.0); + VERIFY(std::tr1::cyl_bessel_j(1.0L, xinf_l) == 0.0L); + + VERIFY(std::tr1::cyl_bessel_j(2.5, xinf_d) == 0.0); +} + +int +main() +{ + test01(); + return 0; +} diff --git a/libstdc++-v3/testsuite/tr1/5_numerical_facilities/special_functions/11_cyl_neumann/check_inf.cc b/libstdc++-v3/testsuite/tr1/5_numerical_facilities/special_functions/11_cyl_neumann/check_inf.cc new file mode 100644 index 000000000000..9294cab64db0 --- /dev/null +++ b/libstdc++-v3/testsuite/tr1/5_numerical_facilities/special_functions/11_cyl_neumann/check_inf.cc @@ -0,0 +1,37 @@ +// { dg-require-c-std "" } +// { dg-add-options ieee } + +// TR1 5.2.1.11 cyl_neumann +// N_nu(+inf) == 0 for all finite nu >= 0 + +#include <tr1/cmath> +#include <limits> +#include <testsuite_hooks.h> + +void +test01() +{ + float xinf_f = std::numeric_limits<float>::infinity(); + double xinf_d = std::numeric_limits<double>::infinity(); + long double xinf_l = std::numeric_limits<long double>::infinity(); + + VERIFY(std::tr1::cyl_neumann(0.0F, xinf_f) == 0.0F); + VERIFY(std::tr1::cyl_neumannf(0.0F, xinf_f) == 0.0F); + VERIFY(std::tr1::cyl_neumann(0.0, xinf_d) == 0.0); + VERIFY(std::tr1::cyl_neumann(0.0L, xinf_l) == 0.0L); + VERIFY(std::tr1::cyl_neumannl(0.0L, xinf_l) == 0.0L); + + VERIFY(std::tr1::cyl_neumann(1.0F, xinf_f) == 0.0F); + VERIFY(std::tr1::cyl_neumann(1.0, xinf_d) == 0.0); + VERIFY(std::tr1::cyl_neumann(1.0L, xinf_l) == 0.0L); + + VERIFY(std::tr1::cyl_neumann(2.5, xinf_d) == 0.0); +} + +int +main() +{ + test01(); + return 0; +} +
