From: Michele Pezzutti Date: Sun, 18 Nov 2018 18:32:26 +0000 (+0100) Subject: PR libstdc++/83566 - cyl_bessel_j returns wrong result for x>1000 X-Git-Url: https://git.libre-soc.org/?a=commitdiff_plain;h=bee39274cb9ffb5415ab50fa102673b0394d4f39;p=gcc.git PR libstdc++/83566 - cyl_bessel_j returns wrong result for x>1000 2018-11-16 Michele Pezzutti Edward Smith-Rowland <3dw4rd@verizon.net> PR libstdc++/83566 - cyl_bessel_j returns wrong result for x>1000 for high orders. * include/tr1/bessel_function.tcc: Perform no fewer than nu/2 iterations of the asymptotic series (nu is the Bessel order). * testsuite/tr1/5_numerical_facilities/special_functions/ 09_cyl_bessel_j/check_value.cc: Add tests at nu=100, 1000<=x<=2000. * testsuite/tr1/5_numerical_facilities/special_functions/ 11_cyl_neumann/check_value.cc: Ditto. * testsuite/special_functions/08_cyl_bessel_j/check_value.cc: Ditto. * testsuite/special_functions/10_cyl_neumann/check_value.cc: Ditto. Co-Authored-By: Edward Smith-Rowland <3dw4rd@verizon.net> From-SVN: r266252 --- diff --git a/libstdc++-v3/ChangeLog b/libstdc++-v3/ChangeLog index b3940b2853c..56c1199f381 100644 --- a/libstdc++-v3/ChangeLog +++ b/libstdc++-v3/ChangeLog @@ -1,3 +1,17 @@ +2018-11-18 Michele Pezzutti + Edward Smith-Rowland <3dw4rd@verizon.net> + + PR libstdc++/83566 - cyl_bessel_j returns wrong result for x>1000 + for high orders. + * include/tr1/bessel_function.tcc: Perform no fewer than nu/2 iterations + of the asymptotic series (nu is the Bessel order). + * testsuite/tr1/5_numerical_facilities/special_functions/ + 09_cyl_bessel_j/check_value.cc: Add tests at nu=100, 1000<=x<=2000. + * testsuite/tr1/5_numerical_facilities/special_functions/ + 11_cyl_neumann/check_value.cc: Ditto. + * testsuite/special_functions/08_cyl_bessel_j/check_value.cc: Ditto. + * testsuite/special_functions/10_cyl_neumann/check_value.cc: Ditto. + 2018-11-17 Jonathan Wakely Implement std::pmr::synchronized_pool_resource diff --git a/libstdc++-v3/include/tr1/bessel_function.tcc b/libstdc++-v3/include/tr1/bessel_function.tcc index 26c66cabe29..7b94b84a6ea 100644 --- a/libstdc++-v3/include/tr1/bessel_function.tcc +++ b/libstdc++-v3/include/tr1/bessel_function.tcc @@ -27,6 +27,10 @@ * Do not attempt to use it directly. @headername{tr1/cmath} */ +/* __cyl_bessel_jn_asymp adapted from GNU GSL version 2.4 specfunc/bessel_j.c + * Copyright (C) 1996-2003 Gerard Jungman + */ + // // ISO C++ 14882 TR1: 5.2 Special functions // @@ -358,24 +362,51 @@ namespace tr1 void __cyl_bessel_jn_asymp(_Tp __nu, _Tp __x, _Tp & __Jnu, _Tp & __Nnu) { - const _Tp __mu = _Tp(4) * __nu * __nu; - const _Tp __mum1 = __mu - _Tp(1); - const _Tp __mum9 = __mu - _Tp(9); - const _Tp __mum25 = __mu - _Tp(25); - const _Tp __mum49 = __mu - _Tp(49); - const _Tp __xx = _Tp(64) * __x * __x; - const _Tp __P = _Tp(1) - __mum1 * __mum9 / (_Tp(2) * __xx) - * (_Tp(1) - __mum25 * __mum49 / (_Tp(12) * __xx)); - const _Tp __Q = __mum1 / (_Tp(8) * __x) - * (_Tp(1) - __mum9 * __mum25 / (_Tp(6) * __xx)); + const _Tp __mu = _Tp(4) * __nu * __nu; + const _Tp __8x = _Tp(8) * __x; + + _Tp __P = _Tp(0); + _Tp __Q = _Tp(0); + + _Tp __k = _Tp(0); + _Tp __term = _Tp(1); + + int __epsP = 0; + int __epsQ = 0; + + _Tp __eps = std::numeric_limits<_Tp>::epsilon(); + + do + { + __term *= (__k == 0 + ? _Tp(1) + : -(__mu - (2 * __k - 1) * (2 * __k - 1)) / (__k * __8x)); + + __epsP = std::abs(__term) < __eps * std::abs(__P); + __P += __term; + + __k++; + + __term *= (__mu - (2 * __k - 1) * (2 * __k - 1)) / (__k * __8x); + __epsQ = std::abs(__term) < __eps * std::abs(__Q); + __Q += __term; + + if (__epsP && __epsQ && __k > (__nu / 2.)) + break; + + __k++; + } + while (__k < 1000); const _Tp __chi = __x - (__nu + _Tp(0.5L)) - * __numeric_constants<_Tp>::__pi_2(); + * __numeric_constants<_Tp>::__pi_2(); + const _Tp __c = std::cos(__chi); const _Tp __s = std::sin(__chi); const _Tp __coef = std::sqrt(_Tp(2) / (__numeric_constants<_Tp>::__pi() * __x)); + __Jnu = __coef * (__c * __P - __s * __Q); __Nnu = __coef * (__s * __P + __c * __Q); diff --git a/libstdc++-v3/testsuite/special_functions/08_cyl_bessel_j/check_value.cc b/libstdc++-v3/testsuite/special_functions/08_cyl_bessel_j/check_value.cc index 475c6dcaa4a..7900af3c959 100644 --- a/libstdc++-v3/testsuite/special_functions/08_cyl_bessel_j/check_value.cc +++ b/libstdc++-v3/testsuite/special_functions/08_cyl_bessel_j/check_value.cc @@ -698,6 +698,39 @@ data026[21] = }; const double toler026 = 1.0000000000000006e-11; +// Test data for nu=100.0000000000000000 +// max(|f - f_GSL|): 3.9438938226332709e-14 at index 19 +// max(|f - f_GSL| / |f_GSL|): 2.0193411077170867e-11 +// mean(f - f_GSL): 1.6682360684660055e-15 +// variance(f - f_GSL): 5.3274331668346898e-28 +// stddev(f - f_GSL): 2.3081232997469372e-14 +const testcase_cyl_bessel_j +data027[21] = +{ + { 1.1676135007789573e-02, 100.0000000000000000, 1000.0000000000000000, 0.0 }, + { -1.1699854778025796e-02, 100.0000000000000000, 1100.0000000000000000, 0.0 }, + { -2.2801483405083697e-02, 100.0000000000000000, 1200.0000000000000000, 0.0 }, + { -1.6973500787373915e-02, 100.0000000000000000, 1300.0000000000000000, 0.0 }, + { -1.4154528803481308e-03, 100.0000000000000000, 1400.0000000000000000, 0.0 }, + { 1.3333726584495232e-02, 100.0000000000000000, 1500.0000000000000000, 0.0 }, + { 1.9802562020148559e-02, 100.0000000000000000, 1600.0000000000000000, 0.0 }, + { 1.6129771279838816e-02, 100.0000000000000000, 1700.0000000000000000, 0.0 }, + { 5.3753369281536031e-03, 100.0000000000000000, 1800.0000000000000000, 0.0 }, + { -6.9238868725645785e-03, 100.0000000000000000, 1900.0000000000000000, 0.0 }, + { -1.5487871720069789e-02, 100.0000000000000000, 2000.0000000000000000, 0.0 }, + { -1.7275186717671070e-02, 100.0000000000000000, 2100.0000000000000000, 0.0 }, + { -1.2233030525173150e-02, 100.0000000000000000, 2200.0000000000000000, 0.0 }, + { -2.8518508672241900e-03, 100.0000000000000000, 2300.0000000000000000, 0.0 }, + { 7.0784372270289329e-03, 100.0000000000000000, 2400.0000000000000000, 0.0 }, + { 1.3955367586928166e-02, 100.0000000000000000, 2500.0000000000000000, 0.0 }, + { 1.5574059842493392e-02, 100.0000000000000000, 2600.0000000000000000, 0.0 }, + { 1.1718043044647556e-02, 100.0000000000000000, 2700.0000000000000000, 0.0 }, + { 4.0320953231285607e-03, 100.0000000000000000, 2800.0000000000000000, 0.0 }, + { -4.6895111783053977e-03, 100.0000000000000000, 2900.0000000000000000, 0.0 }, + { -1.1507715400035966e-02, 100.0000000000000000, 3000.0000000000000000, 0.0 }, +}; +const double toler027 = 1.0000000000000006e-10; + template void test(const testcase_cyl_bessel_j (&data)[Num], Ret toler) @@ -748,5 +781,6 @@ main() test(data024, toler024); test(data025, toler025); test(data026, toler026); + test(data027, toler027); return 0; } diff --git a/libstdc++-v3/testsuite/special_functions/10_cyl_neumann/check_value.cc b/libstdc++-v3/testsuite/special_functions/10_cyl_neumann/check_value.cc index da1250278d5..835ce8cd676 100644 --- a/libstdc++-v3/testsuite/special_functions/10_cyl_neumann/check_value.cc +++ b/libstdc++-v3/testsuite/special_functions/10_cyl_neumann/check_value.cc @@ -742,6 +742,39 @@ data028[20] = }; const double toler028 = 1.0000000000000006e-11; +// Test data for nu=100.0000000000000000 +// max(|f - f_GSL|): 3.9022387751663778e-14 at index 16 +// max(|f - f_GSL| / |f_GSL|): 2.4760677072012703e-11 +// mean(f - f_GSL): 3.6878362466971231e-16 +// variance(f - f_GSL): 5.0707962306468580e-28 +// stddev(f - f_GSL): 2.2518428521206487e-14 +const testcase_cyl_neumann +data029[21] = +{ + { -2.2438688257729954e-02, 100.0000000000000000, 1000.0000000000000000, 0.0 }, + { -2.1077595159819992e-02, 100.0000000000000000, 1100.0000000000000000, 0.0 }, + { -3.5299439206692585e-03, 100.0000000000000000, 1200.0000000000000000, 0.0 }, + { 1.4250019326536615e-02, 100.0000000000000000, 1300.0000000000000000, 0.0 }, + { 2.1304679089735663e-02, 100.0000000000000000, 1400.0000000000000000, 0.0 }, + { 1.5734395077905267e-02, 100.0000000000000000, 1500.0000000000000000, 0.0 }, + { 2.5544633636137774e-03, 100.0000000000000000, 1600.0000000000000000, 0.0 }, + { -1.0722045524849367e-02, 100.0000000000000000, 1700.0000000000000000, 0.0 }, + { -1.8036919243226864e-02, 100.0000000000000000, 1800.0000000000000000, 0.0 }, + { -1.6958415593079763e-02, 100.0000000000000000, 1900.0000000000000000, 0.0 }, + { -8.8788704566276667e-03, 100.0000000000000000, 2000.0000000000000000, 0.0 }, + { 2.2504407108413179e-03, 100.0000000000000000, 2100.0000000000000000, 0.0 }, + { 1.1833215246712251e-02, 100.0000000000000000, 2200.0000000000000000, 0.0 }, + { 1.6398784536343945e-02, 100.0000000000000000, 2300.0000000000000000, 0.0 }, + { 1.4675984403642338e-02, 100.0000000000000000, 2400.0000000000000000, 0.0 }, + { 7.7523920451654229e-03, 100.0000000000000000, 2500.0000000000000000, 0.0 }, + { -1.5759822576003489e-03, 100.0000000000000000, 2600.0000000000000000, 0.0 }, + { -9.9314877404787089e-03, 100.0000000000000000, 2700.0000000000000000, 0.0 }, + { -1.4534495161704743e-02, 100.0000000000000000, 2800.0000000000000000, 0.0 }, + { -1.4059273497237509e-02, 100.0000000000000000, 2900.0000000000000000, 0.0 }, + { -8.9385158149605185e-03, 100.0000000000000000, 3000.0000000000000000, 0.0 }, +}; +const double toler029 = 1.0000000000000006e-10; + template void test(const testcase_cyl_neumann (&data)[Num], Ret toler) @@ -794,5 +827,6 @@ main() test(data026, toler026); test(data027, toler027); test(data028, toler028); + test(data029, toler029); return 0; } diff --git a/libstdc++-v3/testsuite/tr1/5_numerical_facilities/special_functions/09_cyl_bessel_j/check_value.cc b/libstdc++-v3/testsuite/tr1/5_numerical_facilities/special_functions/09_cyl_bessel_j/check_value.cc index 06a2d192df0..ee587cbad86 100644 --- a/libstdc++-v3/testsuite/tr1/5_numerical_facilities/special_functions/09_cyl_bessel_j/check_value.cc +++ b/libstdc++-v3/testsuite/tr1/5_numerical_facilities/special_functions/09_cyl_bessel_j/check_value.cc @@ -698,6 +698,39 @@ data026[21] = }; const double toler026 = 1.0000000000000006e-11; +// Test data for nu=100.0000000000000000 +// max(|f - f_GSL|): 3.9438938226332709e-14 at index 19 +// max(|f - f_GSL| / |f_GSL|): 2.0193411077170867e-11 +// mean(f - f_GSL): 1.6682360684660055e-15 +// variance(f - f_GSL): 5.3274331668346898e-28 +// stddev(f - f_GSL): 2.3081232997469372e-14 +const testcase_cyl_bessel_j +data027[21] = +{ + { 1.1676135007789573e-02, 100.0000000000000000, 1000.0000000000000000, 0.0 }, + { -1.1699854778025796e-02, 100.0000000000000000, 1100.0000000000000000, 0.0 }, + { -2.2801483405083697e-02, 100.0000000000000000, 1200.0000000000000000, 0.0 }, + { -1.6973500787373915e-02, 100.0000000000000000, 1300.0000000000000000, 0.0 }, + { -1.4154528803481308e-03, 100.0000000000000000, 1400.0000000000000000, 0.0 }, + { 1.3333726584495232e-02, 100.0000000000000000, 1500.0000000000000000, 0.0 }, + { 1.9802562020148559e-02, 100.0000000000000000, 1600.0000000000000000, 0.0 }, + { 1.6129771279838816e-02, 100.0000000000000000, 1700.0000000000000000, 0.0 }, + { 5.3753369281536031e-03, 100.0000000000000000, 1800.0000000000000000, 0.0 }, + { -6.9238868725645785e-03, 100.0000000000000000, 1900.0000000000000000, 0.0 }, + { -1.5487871720069789e-02, 100.0000000000000000, 2000.0000000000000000, 0.0 }, + { -1.7275186717671070e-02, 100.0000000000000000, 2100.0000000000000000, 0.0 }, + { -1.2233030525173150e-02, 100.0000000000000000, 2200.0000000000000000, 0.0 }, + { -2.8518508672241900e-03, 100.0000000000000000, 2300.0000000000000000, 0.0 }, + { 7.0784372270289329e-03, 100.0000000000000000, 2400.0000000000000000, 0.0 }, + { 1.3955367586928166e-02, 100.0000000000000000, 2500.0000000000000000, 0.0 }, + { 1.5574059842493392e-02, 100.0000000000000000, 2600.0000000000000000, 0.0 }, + { 1.1718043044647556e-02, 100.0000000000000000, 2700.0000000000000000, 0.0 }, + { 4.0320953231285607e-03, 100.0000000000000000, 2800.0000000000000000, 0.0 }, + { -4.6895111783053977e-03, 100.0000000000000000, 2900.0000000000000000, 0.0 }, + { -1.1507715400035966e-02, 100.0000000000000000, 3000.0000000000000000, 0.0 }, +}; +const double toler027 = 1.0000000000000006e-10; + template void test(const testcase_cyl_bessel_j (&data)[Num], Ret toler) @@ -748,5 +781,6 @@ main() test(data024, toler024); test(data025, toler025); test(data026, toler026); + test(data027, toler027); return 0; } diff --git a/libstdc++-v3/testsuite/tr1/5_numerical_facilities/special_functions/11_cyl_neumann/check_value.cc b/libstdc++-v3/testsuite/tr1/5_numerical_facilities/special_functions/11_cyl_neumann/check_value.cc index 3c0cbe7b977..d73f6a604fe 100644 --- a/libstdc++-v3/testsuite/tr1/5_numerical_facilities/special_functions/11_cyl_neumann/check_value.cc +++ b/libstdc++-v3/testsuite/tr1/5_numerical_facilities/special_functions/11_cyl_neumann/check_value.cc @@ -742,6 +742,39 @@ data028[20] = }; const double toler028 = 1.0000000000000006e-11; +// Test data for nu=100.0000000000000000 +// max(|f - f_GSL|): 3.9022387751663778e-14 at index 16 +// max(|f - f_GSL| / |f_GSL|): 2.4760677072012703e-11 +// mean(f - f_GSL): 3.6878362466971231e-16 +// variance(f - f_GSL): 5.0707962306468580e-28 +// stddev(f - f_GSL): 2.2518428521206487e-14 +const testcase_cyl_neumann +data029[21] = +{ + { -2.2438688257729954e-02, 100.0000000000000000, 1000.0000000000000000, 0.0 }, + { -2.1077595159819992e-02, 100.0000000000000000, 1100.0000000000000000, 0.0 }, + { -3.5299439206692585e-03, 100.0000000000000000, 1200.0000000000000000, 0.0 }, + { 1.4250019326536615e-02, 100.0000000000000000, 1300.0000000000000000, 0.0 }, + { 2.1304679089735663e-02, 100.0000000000000000, 1400.0000000000000000, 0.0 }, + { 1.5734395077905267e-02, 100.0000000000000000, 1500.0000000000000000, 0.0 }, + { 2.5544633636137774e-03, 100.0000000000000000, 1600.0000000000000000, 0.0 }, + { -1.0722045524849367e-02, 100.0000000000000000, 1700.0000000000000000, 0.0 }, + { -1.8036919243226864e-02, 100.0000000000000000, 1800.0000000000000000, 0.0 }, + { -1.6958415593079763e-02, 100.0000000000000000, 1900.0000000000000000, 0.0 }, + { -8.8788704566276667e-03, 100.0000000000000000, 2000.0000000000000000, 0.0 }, + { 2.2504407108413179e-03, 100.0000000000000000, 2100.0000000000000000, 0.0 }, + { 1.1833215246712251e-02, 100.0000000000000000, 2200.0000000000000000, 0.0 }, + { 1.6398784536343945e-02, 100.0000000000000000, 2300.0000000000000000, 0.0 }, + { 1.4675984403642338e-02, 100.0000000000000000, 2400.0000000000000000, 0.0 }, + { 7.7523920451654229e-03, 100.0000000000000000, 2500.0000000000000000, 0.0 }, + { -1.5759822576003489e-03, 100.0000000000000000, 2600.0000000000000000, 0.0 }, + { -9.9314877404787089e-03, 100.0000000000000000, 2700.0000000000000000, 0.0 }, + { -1.4534495161704743e-02, 100.0000000000000000, 2800.0000000000000000, 0.0 }, + { -1.4059273497237509e-02, 100.0000000000000000, 2900.0000000000000000, 0.0 }, + { -8.9385158149605185e-03, 100.0000000000000000, 3000.0000000000000000, 0.0 }, +}; +const double toler029 = 1.0000000000000006e-10; + template void test(const testcase_cyl_neumann (&data)[Num], Ret toler) @@ -794,5 +827,6 @@ main() test(data026, toler026); test(data027, toler027); test(data028, toler028); + test(data029, toler029); return 0; }