Message ID | 5227844d-a96b-b587-a599-51f698718c09@tiscali.it |
---|---|
State | New |
Headers | show |
Series | Fix Bug 83566 - cyl_bessel_j returns wrong result for x>1000 for high orders | expand |
On 12/31/2017 09:38 PM, Michele Pezzutti wrote: > Hi. > > This patch intends to fix Bug 83566 - cyl_bessel_j returns wrong > result for x>1000 for high orders. > Seehttps://gcc.gnu.org/bugzilla/show_bug.cgi?id=83566 forissue > description. > > * libstdc++-v3/include/tr1/bessel_function.tcc > Series expansion in __cyl_bessel_jn_asymp() shall not be > truncated at the first terms. > At least nu/2 terms shall be added, in order to have a > guaranteed error bound. > > * > libstdc++-v3/testsuite/tr1/5_numerical_facilities/special_functions/09_cyl_bessel_j/check_value.cc > Testcase for x > 1000 added. > > * > libstdc++-v3/testsuite/tr1/5_numerical_facilities/special_functions/11_cyl_neumann/check_value.cc > Testcase for x > 1000 added. > > > diff --git a/libstdc++-v3/include/tr1/bessel_function.tcc > b/libstdc++-v3/include/tr1/bessel_function.tcc > index 7ac733d..852a31e 100644 > --- a/libstdc++-v3/include/tr1/bessel_function.tcc > +++ b/libstdc++-v3/include/tr1/bessel_function.tcc > @@ -359,15 +359,32 @@ namespace tr1 > __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 __8x = _Tp(8) * __x; > + > + _Tp __P = _Tp(1); > + _Tp __Q = _Tp(0); > + > + _Tp __eps = std::numeric_limits<_Tp>::epsilon(); > + > + _Tp __term = _Tp(1); > + _Tp __epsP = _Tp(0); > + _Tp __epsQ = _Tp(0); > + > + unsigned long __2k; > + for (unsigned long k = 1; k < 1000; k+=2) { > + __2k = 2 * k; > + > + __term *= (__mu - _Tp((__2k - 1) * (__2k - 1))) / (k * __8x); > + __Q += __term; > + __epsQ = std::abs(__term) < std::abs(__eps * __Q); > + > + __term *= -(__mu - _Tp((__2k + 1) * (__2k + 1)))/ (_Tp(k + > 1) * __8x); > + __P += __term; > + __epsP = std::abs(__term) < std::abs(__eps * __P); > + > + if (__epsP && __epsQ && k > __nu / 2.) > + break; > + } > > const _Tp __chi = __x - (__nu + _Tp(0.5L)) > * __numeric_constants<_Tp>::__pi_2(); > > > 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 26f4dd3..e340b78 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,26 @@ data026[21] = > }; > const double toler026 = 1.0000000000000006e-11; > > +// Test data for nu=100.00000000000000. > +// max(|f - f_GSL|): 2.5857788132910287e-14 > +// max(|f - f_GSL| / |f_GSL|): 1.6767662425535933e-11 > +const testcase_cyl_bessel_j<double> > +data027[11] = > +{ > + { 0.0116761350077845, 100.0000000000000000, 1000.0000000000000000 }, > + {-0.0116998547780258, 100.0000000000000000, 1100.0000000000000000 }, > + {-0.0228014834050837, 100.0000000000000000, 1200.0000000000000000 }, > + {-0.0169735007873739, 100.0000000000000000, 1300.0000000000000000 }, > + {-0.0014154528803530, 100.0000000000000000, 1400.0000000000000000 }, > + { 0.0133337265844988, 100.0000000000000000, 1500.0000000000000000 }, > + { 0.0198025620201474, 100.0000000000000000, 1600.0000000000000000 }, > + { 0.0161297712798388, 100.0000000000000000, 1700.0000000000000000 }, > + { 0.0053753369281577, 100.0000000000000000, 1800.0000000000000000 }, > + {-0.0069238868725646, 100.0000000000000000, 1900.0000000000000000 }, > + {-0.0154878717200738, 100.0000000000000000, 2000.0000000000000000 }, > +}; > +const double toler027 = 1.0000000000000006e-10; > + > template<typename Ret, unsigned int Num> > void > test(const testcase_cyl_bessel_j<Ret> (&data)[Num], Ret toler) > @@ -748,5 +768,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 5579149..9caf836 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,26 @@ data028[20] = > }; > const double toler028 = 1.0000000000000006e-11; > > +// Test data for nu=100.00000000000000. > +// max(|f - f_GSL|): 3.1049815496508870e-14 > +// max(|f - f_GSL| / |f_GSL|): 1.6767662425535933e-11 > +const testcase_cyl_neumann<double> > +data029[11] = > +{ > + {-0.0224386882577326, 100.0000000000000000, 1000.0000000000000000 }, > + {-0.0210775951598200, 100.0000000000000000, 1100.0000000000000000 }, > + {-0.0035299439206693, 100.0000000000000000, 1200.0000000000000000 }, > + { 0.0142500193265366, 100.0000000000000000, 1300.0000000000000000 }, > + { 0.0213046790897353, 100.0000000000000000, 1400.0000000000000000 }, > + { 0.0157343950779022, 100.0000000000000000, 1500.0000000000000000 }, > + { 0.0025544633636228, 100.0000000000000000, 1600.0000000000000000 }, > + {-0.0107220455248494, 100.0000000000000000, 1700.0000000000000000 }, > + {-0.0180369192432256, 100.0000000000000000, 1800.0000000000000000 }, > + {-0.0169584155930798, 100.0000000000000000, 1900.0000000000000000 }, > + {-0.0088788704566206, 100.0000000000000000, 2000.0000000000000000 }, > +}; > +const double toler029 = 1.0000000000000006e-10; > + > template<typename Ret, unsigned int Num> > void > test(const testcase_cyl_neumann<Ret> (&data)[Num], Ret toler) > @@ -794,5 +814,6 @@ main() > test(data026, toler026); > test(data027, toler027); > test(data028, toler028); > + test(data029, toler029); > return 0; > } > I like the patch. I have a similar one in the tr29124 branch. Anyway, I got held up and I think it's good to have new folks looking into this. It looks good except that you need to uglify k. Also, keep in mind that these series are asymptotic - they'll eventually blow up. You should watch the magnitude of sequential terms and bail out of the loop if either term starts growing. Ed
On 12/31/2017 09:38 PM, Michele Pezzutti wrote: > Hi. > > This patch intends to fix Bug 83566 - cyl_bessel_j returns wrong > result for x>1000 for high orders. > Seehttps://gcc.gnu.org/bugzilla/show_bug.cgi?id=83566 forissue > description. > > * libstdc++-v3/include/tr1/bessel_function.tcc > Series expansion in __cyl_bessel_jn_asymp() shall not be > truncated at the first terms. > At least nu/2 terms shall be added, in order to have a > guaranteed error bound. > > * > libstdc++-v3/testsuite/tr1/5_numerical_facilities/special_functions/09_cyl_bessel_j/check_value.cc > Testcase for x > 1000 added. > > * > libstdc++-v3/testsuite/tr1/5_numerical_facilities/special_functions/11_cyl_neumann/check_value.cc > Testcase for x > 1000 added. > > > diff --git a/libstdc++-v3/include/tr1/bessel_function.tcc > b/libstdc++-v3/include/tr1/bessel_function.tcc > index 7ac733d..852a31e 100644 > --- a/libstdc++-v3/include/tr1/bessel_function.tcc > +++ b/libstdc++-v3/include/tr1/bessel_function.tcc > @@ -359,15 +359,32 @@ namespace tr1 > __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 __8x = _Tp(8) * __x; > + > + _Tp __P = _Tp(1); > + _Tp __Q = _Tp(0); > + > + _Tp __eps = std::numeric_limits<_Tp>::epsilon(); > + > + _Tp __term = _Tp(1); > + _Tp __epsP = _Tp(0); > + _Tp __epsQ = _Tp(0); > + > + unsigned long __2k; > + for (unsigned long k = 1; k < 1000; k+=2) { > + __2k = 2 * k; > + > + __term *= (__mu - _Tp((__2k - 1) * (__2k - 1))) / (k * __8x); > + __Q += __term; > + __epsQ = std::abs(__term) < std::abs(__eps * __Q); > + > + __term *= -(__mu - _Tp((__2k + 1) * (__2k + 1)))/ (_Tp(k + > 1) * __8x); > + __P += __term; > + __epsP = std::abs(__term) < std::abs(__eps * __P); > + > + if (__epsP && __epsQ && k > __nu / 2.) > + break; > + } > > const _Tp __chi = __x - (__nu + _Tp(0.5L)) > * __numeric_constants<_Tp>::__pi_2(); > > > 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 26f4dd3..e340b78 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,26 @@ data026[21] = > }; > const double toler026 = 1.0000000000000006e-11; > > +// Test data for nu=100.00000000000000. > +// max(|f - f_GSL|): 2.5857788132910287e-14 > +// max(|f - f_GSL| / |f_GSL|): 1.6767662425535933e-11 > +const testcase_cyl_bessel_j<double> > +data027[11] = > +{ > + { 0.0116761350077845, 100.0000000000000000, 1000.0000000000000000 }, > + {-0.0116998547780258, 100.0000000000000000, 1100.0000000000000000 }, > + {-0.0228014834050837, 100.0000000000000000, 1200.0000000000000000 }, > + {-0.0169735007873739, 100.0000000000000000, 1300.0000000000000000 }, > + {-0.0014154528803530, 100.0000000000000000, 1400.0000000000000000 }, > + { 0.0133337265844988, 100.0000000000000000, 1500.0000000000000000 }, > + { 0.0198025620201474, 100.0000000000000000, 1600.0000000000000000 }, > + { 0.0161297712798388, 100.0000000000000000, 1700.0000000000000000 }, > + { 0.0053753369281577, 100.0000000000000000, 1800.0000000000000000 }, > + {-0.0069238868725646, 100.0000000000000000, 1900.0000000000000000 }, > + {-0.0154878717200738, 100.0000000000000000, 2000.0000000000000000 }, > +}; > +const double toler027 = 1.0000000000000006e-10; > + > template<typename Ret, unsigned int Num> > void > test(const testcase_cyl_bessel_j<Ret> (&data)[Num], Ret toler) > @@ -748,5 +768,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 5579149..9caf836 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,26 @@ data028[20] = > }; > const double toler028 = 1.0000000000000006e-11; > > +// Test data for nu=100.00000000000000. > +// max(|f - f_GSL|): 3.1049815496508870e-14 > +// max(|f - f_GSL| / |f_GSL|): 1.6767662425535933e-11 > +const testcase_cyl_neumann<double> > +data029[11] = > +{ > + {-0.0224386882577326, 100.0000000000000000, 1000.0000000000000000 }, > + {-0.0210775951598200, 100.0000000000000000, 1100.0000000000000000 }, > + {-0.0035299439206693, 100.0000000000000000, 1200.0000000000000000 }, > + { 0.0142500193265366, 100.0000000000000000, 1300.0000000000000000 }, > + { 0.0213046790897353, 100.0000000000000000, 1400.0000000000000000 }, > + { 0.0157343950779022, 100.0000000000000000, 1500.0000000000000000 }, > + { 0.0025544633636228, 100.0000000000000000, 1600.0000000000000000 }, > + {-0.0107220455248494, 100.0000000000000000, 1700.0000000000000000 }, > + {-0.0180369192432256, 100.0000000000000000, 1800.0000000000000000 }, > + {-0.0169584155930798, 100.0000000000000000, 1900.0000000000000000 }, > + {-0.0088788704566206, 100.0000000000000000, 2000.0000000000000000 }, > +}; > +const double toler029 = 1.0000000000000006e-10; > + > template<typename Ret, unsigned int Num> > void > test(const testcase_cyl_neumann<Ret> (&data)[Num], Ret toler) > @@ -794,5 +814,6 @@ main() > test(data026, toler026); > test(data027, toler027); > test(data028, toler028); > + test(data029, toler029); > return 0; > } > As a second project look at the location where the number of iterations is exceeded in the main bessel routine. Instead of throwing, just call the asymp routine. Ed
On 01/02/2018 02:28 AM, Ed Smith-Rowland wrote: > I like the patch. > > I have a similar one in the tr29124 branch. > > Anyway, I got held up and I think it's good to have new folks looking > into this. > > It looks good except that you need to uglify k. I looked at the GSL implementation, based on same reference, and their loop is cleaner. What about porting that implementation here? Possible? My implementation is also using one more term for P than for Q, which is discouraged in GSL, according to their comments. > Also, keep in mind that these series are asymptotic - they'll > eventually blow up. > > You should watch the magnitude of sequential terms and bail out of the > loop if either term starts growing. > > Ed Yes, you are right. As nu increases, the multiplying '__term' gets larger, and the result will lose accuracy. This cannot be used for large nu. I have not been able to figure out how to detect when it starts to lose accuracy though, other than empirically. GSL has no checks on terms and uses this method only under certain conditions, on nu and x. Otherwise, they use other methods.
On 12/31/2017 09:38 PM, Michele Pezzutti wrote: > Hi. > > This patch intends to fix Bug 83566 - cyl_bessel_j returns wrong > result for x>1000 for high orders. > Seehttps://gcc.gnu.org/bugzilla/show_bug.cgi?id=83566 forissue > description. > > * libstdc++-v3/include/tr1/bessel_function.tcc > Series expansion in __cyl_bessel_jn_asymp() shall not be > truncated at the first terms. > At least nu/2 terms shall be added, in order to have a > guaranteed error bound. > > * > libstdc++-v3/testsuite/tr1/5_numerical_facilities/special_functions/09_cyl_bessel_j/check_value.cc > Testcase for x > 1000 added. > > * > libstdc++-v3/testsuite/tr1/5_numerical_facilities/special_functions/11_cyl_neumann/check_value.cc > Testcase for x > 1000 added. > > > diff --git a/libstdc++-v3/include/tr1/bessel_function.tcc > b/libstdc++-v3/include/tr1/bessel_function.tcc > index 7ac733d..852a31e 100644 > --- a/libstdc++-v3/include/tr1/bessel_function.tcc > +++ b/libstdc++-v3/include/tr1/bessel_function.tcc > @@ -359,15 +359,32 @@ namespace tr1 > __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 __8x = _Tp(8) * __x; > + > + _Tp __P = _Tp(1); > + _Tp __Q = _Tp(0); > + > + _Tp __eps = std::numeric_limits<_Tp>::epsilon(); > + > + _Tp __term = _Tp(1); > + _Tp __epsP = _Tp(0); > + _Tp __epsQ = _Tp(0); > + > + unsigned long __2k; > + for (unsigned long k = 1; k < 1000; k+=2) { > + __2k = 2 * k; > + > + __term *= (__mu - _Tp((__2k - 1) * (__2k - 1))) / (k * __8x); > + __Q += __term; > + __epsQ = std::abs(__term) < std::abs(__eps * __Q); > + > + __term *= -(__mu - _Tp((__2k + 1) * (__2k + 1)))/ (_Tp(k + > 1) * __8x); > + __P += __term; > + __epsP = std::abs(__term) < std::abs(__eps * __P); > + > + if (__epsP && __epsQ && k > __nu / 2.) > + break; > + } > > const _Tp __chi = __x - (__nu + _Tp(0.5L)) > * __numeric_constants<_Tp>::__pi_2(); > > > 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 26f4dd3..e340b78 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,26 @@ data026[21] = > }; > const double toler026 = 1.0000000000000006e-11; > > +// Test data for nu=100.00000000000000. > +// max(|f - f_GSL|): 2.5857788132910287e-14 > +// max(|f - f_GSL| / |f_GSL|): 1.6767662425535933e-11 > +const testcase_cyl_bessel_j<double> > +data027[11] = > +{ > + { 0.0116761350077845, 100.0000000000000000, 1000.0000000000000000 }, > + {-0.0116998547780258, 100.0000000000000000, 1100.0000000000000000 }, > + {-0.0228014834050837, 100.0000000000000000, 1200.0000000000000000 }, > + {-0.0169735007873739, 100.0000000000000000, 1300.0000000000000000 }, > + {-0.0014154528803530, 100.0000000000000000, 1400.0000000000000000 }, > + { 0.0133337265844988, 100.0000000000000000, 1500.0000000000000000 }, > + { 0.0198025620201474, 100.0000000000000000, 1600.0000000000000000 }, > + { 0.0161297712798388, 100.0000000000000000, 1700.0000000000000000 }, > + { 0.0053753369281577, 100.0000000000000000, 1800.0000000000000000 }, > + {-0.0069238868725646, 100.0000000000000000, 1900.0000000000000000 }, > + {-0.0154878717200738, 100.0000000000000000, 2000.0000000000000000 }, > +}; > +const double toler027 = 1.0000000000000006e-10; > + > template<typename Ret, unsigned int Num> > void > test(const testcase_cyl_bessel_j<Ret> (&data)[Num], Ret toler) > @@ -748,5 +768,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 5579149..9caf836 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,26 @@ data028[20] = > }; > const double toler028 = 1.0000000000000006e-11; > > +// Test data for nu=100.00000000000000. > +// max(|f - f_GSL|): 3.1049815496508870e-14 > +// max(|f - f_GSL| / |f_GSL|): 1.6767662425535933e-11 > +const testcase_cyl_neumann<double> > +data029[11] = > +{ > + {-0.0224386882577326, 100.0000000000000000, 1000.0000000000000000 }, > + {-0.0210775951598200, 100.0000000000000000, 1100.0000000000000000 }, > + {-0.0035299439206693, 100.0000000000000000, 1200.0000000000000000 }, > + { 0.0142500193265366, 100.0000000000000000, 1300.0000000000000000 }, > + { 0.0213046790897353, 100.0000000000000000, 1400.0000000000000000 }, > + { 0.0157343950779022, 100.0000000000000000, 1500.0000000000000000 }, > + { 0.0025544633636228, 100.0000000000000000, 1600.0000000000000000 }, > + {-0.0107220455248494, 100.0000000000000000, 1700.0000000000000000 }, > + {-0.0180369192432256, 100.0000000000000000, 1800.0000000000000000 }, > + {-0.0169584155930798, 100.0000000000000000, 1900.0000000000000000 }, > + {-0.0088788704566206, 100.0000000000000000, 2000.0000000000000000 }, > +}; > +const double toler029 = 1.0000000000000006e-10; > + > template<typename Ret, unsigned int Num> > void > test(const testcase_cyl_neumann<Ret> (&data)[Num], Ret toler) > @@ -794,5 +814,6 @@ main() > test(data026, toler026); > test(data027, toler027); > test(data028, toler028); > + test(data029, toler029); > return 0; > } > OK, on *third* look summing up to k = nu/2 at minimum will a achieve the result of not blowing up the asymptotic series: nu^2 - (2k-1)^2. And it will do that without a check. This stopping criterion should work even near x=nu which would be the most difficult case. The sum could go further for larger x but let's just go with your termination criterion for now. Later, with some experimentation, we could sum up to nu/2 at a minimum *then* snoop forward until the terms start drifting up. Or we could just solve k_max for this case as a function of x. Also, we may never need these extras. Thanks for doing this. Ed My "tests" in my work area were not actually testing this case. Bah Humbug!
On 01/02/2018 05:59 PM, Ed Smith-Rowland wrote: > OK, > > on *third* look summing up to k = nu/2 at minimum will a achieve the > result of not blowing up the asymptotic series: > > nu^2 - (2k-1)^2. And it will do that without a check. > > This stopping criterion should work even near x=nu which would be the > most difficult case. The sum could go further for larger x but let's > just go with your termination criterion for now. Later, with some > experimentation, we could sum up to nu/2 at a minimum *then* snoop > forward until the terms start drifting up. Or we could just solve > k_max for this case as a function of x. Also, we may never need these > extras. For what I could see, 'term' stops to contributing to P, Q few steps before k = nu /2. > > Thanks for doing this. > > Ed > > My "tests" in my work area were not actually testing this case. Bah > Humbug! > >
On 01/02/2018 04:41 PM, Michele Pezzutti wrote: > On 01/02/2018 05:59 PM, Ed Smith-Rowland wrote: >> OK, >> >> on *third* look summing up to k = nu/2 at minimum will a achieve the >> result of not blowing up the asymptotic series: >> >> nu^2 - (2k-1)^2. And it will do that without a check. >> >> This stopping criterion should work even near x=nu which would be the >> most difficult case. The sum could go further for larger x but let's >> just go with your termination criterion for now. Later, with some >> experimentation, we could sum up to nu/2 at a minimum *then* snoop >> forward until the terms start drifting up. Or we could just solve >> k_max for this case as a function of x. Also, we may never need these >> extras. > For what I could see, 'term' stops to contributing to P, Q few steps > before k = nu /2. > Right. I can see these terms rising, sometimes quite high, then falling. In many experiments, the sums converge before the terms blow up. In my old impl, I gave up too easily. Perhaps, in general, one could wait until the terms started decreasing, *then* if subsequent terms start growing again exit loop. But it seems I never get to that point before sum convergence. On the good side it seems that if x >= 20 nu then you can go quite high in nu. nu = 1000 x = 20000 Jnu = 0.00540683536187055 Nnu = -0.00162387432151849 Jnua = 0.00540683536187055 Nnua = -0.00162387432151849 Jnua - Jnu = 0 Nnua - Nnu = 0 Jpnu = 0.00162170521572447 Npnu = 0.0054001182451475 Jpnua = 0.00162170521572447 Npnua = 0.0054001182451475 Jpnua - Jpnu = 0 Npnua - Npnu = 0 pi x Wronski / 2 = 1.00000022382762
Hi. On 01/02/2018 05:43 PM, Michele Pezzutti wrote: > > On 01/02/2018 02:28 AM, Ed Smith-Rowland wrote: >> I like the patch. >> >> I have a similar one in the tr29124 branch. >> >> Anyway, I got held up and I think it's good to have new folks looking >> into this. >> >> It looks good except that you need to uglify k. > > I looked at the GSL implementation, based on same reference, and their > loop is cleaner. What about porting that implementation here? Possible? > > My implementation is also using one more term for P than for Q, which > is discouraged in GSL, according to their comments. Ed, do you have any comment about this point? Regarding the 2nd line, after my observations, usually term stops contributing to P, Q before k > nu/2, so actually, an offset by one is most likely without consequences. GSL implementation is nevertheless more elegant. > >> Also, keep in mind that these series are asymptotic - they'll >> eventually blow up. >> >> You should watch the magnitude of sequential terms and bail out of >> the loop if either term starts growing. >> >> Ed > Yes, you are right. As nu increases, the multiplying '__term' gets > larger, and the result will lose accuracy. > This cannot be used for large nu. > I have not been able to figure out how to detect when it starts to > lose accuracy though, other than empirically. > GSL has no checks on terms and uses this method only under certain > conditions, on nu and x. > Otherwise, they use other methods.
On 01/03/2018 02:49 PM, Michele Pezzutti wrote: > > Hi. > > On 01/02/2018 05:43 PM, Michele Pezzutti wrote: >> >> On 01/02/2018 02:28 AM, Ed Smith-Rowland wrote: >>> I like the patch. >>> >>> I have a similar one in the tr29124 branch. >>> >>> Anyway, I got held up and I think it's good to have new folks >>> looking into this. >>> >>> It looks good except that you need to uglify k. >> >> I looked at the GSL implementation, based on same reference, and >> their loop is cleaner. What about porting that implementation here? >> Possible? >> >> My implementation is also using one more term for P than for Q, which >> is discouraged in GSL, according to their comments. > > Ed, do you have any comment about this point? > Regarding the 2nd line, after my observations, usually term stops > contributing to P, Q before k > nu/2, so actually, an offset by one is > most likely without consequences. > GSL implementation is nevertheless more elegant. > I've seen their implementation. It's tight. Feel free to port it. I assume this is it: int gsl_sf_bessel_Jnu_asympx_e(const double nu, const double x, gsl_sf_result * result); You do want to give Q or b one more term so as to make that last factor as small as possible. They're right. In principal, these ideas should be ported to I,K but I think (and IIRC GSL agrees) these under,over-flow before they have much effect. I think I put the full series in there. They could use a similar clean-up. But let's get this in there first. Ed
On 01/04/2018 06:16 PM, Ed Smith-Rowland wrote: > On 01/03/2018 02:49 PM, Michele Pezzutti wrote: >> >> Hi. >> >> On 01/02/2018 05:43 PM, Michele Pezzutti wrote: >>> >>> On 01/02/2018 02:28 AM, Ed Smith-Rowland wrote: >>>> I like the patch. >>>> >>>> I have a similar one in the tr29124 branch. >>>> >>>> Anyway, I got held up and I think it's good to have new folks >>>> looking into this. >>>> >>>> It looks good except that you need to uglify k. >>> >>> I looked at the GSL implementation, based on same reference, and >>> their loop is cleaner. What about porting that implementation here? >>> Possible? >>> >>> My implementation is also using one more term for P than for Q, >>> which is discouraged in GSL, according to their comments. >> >> Ed, do you have any comment about this point? >> Regarding the 2nd line, after my observations, usually term stops >> contributing to P, Q before k > nu/2, so actually, an offset by one >> is most likely without consequences. >> GSL implementation is nevertheless more elegant. >> > I've seen their implementation. It's tight. Feel free to port it. > I assume this is it: > int > gsl_sf_bessel_Jnu_asympx_e(const double nu, const double x, > gsl_sf_result * result); > > You do want to give Q or b one more term so as to make that last > factor as small as possible. They're right. > Here it is. diff --git a/libstdc++-v3/include/tr1/bessel_function.tcc b/libstdc++-v3/include/tr1/bessel_function.tcc index 7ac733d..7842350 100644 --- a/libstdc++-v3/include/tr1/bessel_function.tcc +++ b/libstdc++-v3/include/tr1/bessel_function.tcc @@ -353,21 +353,47 @@ namespace tr1 * @param __x The argument of the Bessel functions. * @param __Jnu The output Bessel function of the first kind. * @param __Nnu The output Neumann function (Bessel function of the second kind). + * + * Adapted for libstdc++ from GNU GSL version 2.4 specfunc/bessel_j.c + * Copyright (C) 1996,1997,1998,1999,2000,2001,2002,2003 Gerard Jungman */ template <typename _Tp> 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 __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) < std::abs(__eps * __P); + __P += __term; + + k++; + + __term *= (__mu - (2 * k - 1) * (2 * k - 1)) / (k * __8x); + __epsQ = std::abs(__term) < std::abs(__eps * __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(); > In principal, these ideas should be ported to I,K but I think (and > IIRC GSL agrees) these under,over-flow before they have much effect. > I think I put the full series in there. They could use a similar > clean-up. > > But let's get this in there first. > > Ed > >
And the test cases 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 26f4dd3..e340b78 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,26 @@ data026[21] = }; const double toler026 = 1.0000000000000006e-11; +// Test data for nu=100.00000000000000. +// max(|f - f_GSL|): 2.5857788132910287e-14 +// max(|f - f_GSL| / |f_GSL|): 1.6767662425535933e-11 +const testcase_cyl_bessel_j<double> +data027[11] = +{ + { 0.0116761350077845, 100.0000000000000000, 1000.0000000000000000 }, + {-0.0116998547780258, 100.0000000000000000, 1100.0000000000000000 }, + {-0.0228014834050837, 100.0000000000000000, 1200.0000000000000000 }, + {-0.0169735007873739, 100.0000000000000000, 1300.0000000000000000 }, + {-0.0014154528803530, 100.0000000000000000, 1400.0000000000000000 }, + { 0.0133337265844988, 100.0000000000000000, 1500.0000000000000000 }, + { 0.0198025620201474, 100.0000000000000000, 1600.0000000000000000 }, + { 0.0161297712798388, 100.0000000000000000, 1700.0000000000000000 }, + { 0.0053753369281577, 100.0000000000000000, 1800.0000000000000000 }, + {-0.0069238868725646, 100.0000000000000000, 1900.0000000000000000 }, + {-0.0154878717200738, 100.0000000000000000, 2000.0000000000000000 }, +}; +const double toler027 = 1.0000000000000006e-10; + template<typename Ret, unsigned int Num> void test(const testcase_cyl_bessel_j<Ret> (&data)[Num], Ret toler) @@ -748,5 +768,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 5579149..957aacd 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,26 @@ data028[20] = }; const double toler028 = 1.0000000000000006e-11; +// Test data for nu=100.00000000000000. +// max(|f - f_GSL|): 3.1049815496508870e-14 +// max(|f - f_GSL| / |f_GSL|): 8.4272302674970308e-12 +const testcase_cyl_neumann<double> +data029[11] = +{ + {-0.0224386882577326, 100.0000000000000000, 1000.0000000000000000 }, + {-0.0210775951598200, 100.0000000000000000, 1100.0000000000000000 }, + {-0.0035299439206693, 100.0000000000000000, 1200.0000000000000000 }, + { 0.0142500193265366, 100.0000000000000000, 1300.0000000000000000 }, + { 0.0213046790897353, 100.0000000000000000, 1400.0000000000000000 }, + { 0.0157343950779022, 100.0000000000000000, 1500.0000000000000000 }, + { 0.0025544633636228, 100.0000000000000000, 1600.0000000000000000 }, + {-0.0107220455248494, 100.0000000000000000, 1700.0000000000000000 }, + {-0.0180369192432256, 100.0000000000000000, 1800.0000000000000000 }, + {-0.0169584155930798, 100.0000000000000000, 1900.0000000000000000 }, + {-0.0088788704566206, 100.0000000000000000, 2000.0000000000000000 }, +}; +const double toler029 = 1.0000000000000006e-11; + template<typename Ret, unsigned int Num> void test(const testcase_cyl_neumann<Ret> (&data)[Num], Ret toler) @@ -794,5 +814,6 @@ main() test(data026, toler026); test(data027, toler027); test(data028, toler028); + test(data029, toler029); return 0; } On 01/04/2018 09:54 PM, Michele Pezzutti wrote: > > > On 01/04/2018 06:16 PM, Ed Smith-Rowland wrote: >> I've seen their implementation. It's tight. Feel free to port it. >> I assume this is it: >> int >> gsl_sf_bessel_Jnu_asympx_e(const double nu, const double x, >> gsl_sf_result * result); >> >> You do want to give Q or b one more term so as to make that last >> factor as small as possible. They're right. >> > > Here it is. > > diff --git a/libstdc++-v3/include/tr1/bessel_function.tcc > b/libstdc++-v3/include/tr1/bessel_function.tcc > index 7ac733d..7842350 100644 > --- a/libstdc++-v3/include/tr1/bessel_function.tcc > +++ b/libstdc++-v3/include/tr1/bessel_function.tcc > @@ -353,21 +353,47 @@ namespace tr1 > * @param __x The argument of the Bessel functions. > * @param __Jnu The output Bessel function of the first kind. > * @param __Nnu The output Neumann function (Bessel function > of the second kind). > + * > + * Adapted for libstdc++ from GNU GSL version 2.4 > specfunc/bessel_j.c > + * Copyright (C) 1996,1997,1998,1999,2000,2001,2002,2003 Gerard > Jungman > */ > template <typename _Tp> > 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 __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) < std::abs(__eps * __P); > + __P += __term; > + > + k++; > + > + __term *= (__mu - (2 * k - 1) * (2 * k - 1)) / (k * __8x); > + __epsQ = std::abs(__term) < std::abs(__eps * __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(); > > > >> In principal, these ideas should be ported to I,K but I think (and >> IIRC GSL agrees) these under,over-flow before they have much effect. >> I think I put the full series in there. They could use a similar >> clean-up. >> >> But let's get this in there first. >> >> Ed >> >> >
On 01/04/2018 03:54 PM, Michele Pezzutti wrote: > > > On 01/04/2018 06:16 PM, Ed Smith-Rowland wrote: >> On 01/03/2018 02:49 PM, Michele Pezzutti wrote: >>> >>> Hi. >>> >>> On 01/02/2018 05:43 PM, Michele Pezzutti wrote: >>>> >>>> On 01/02/2018 02:28 AM, Ed Smith-Rowland wrote: >>>>> I like the patch. >>>>> >>>>> I have a similar one in the tr29124 branch. >>>>> >>>>> Anyway, I got held up and I think it's good to have new folks >>>>> looking into this. >>>>> >>>>> It looks good except that you need to uglify k. >>>> >>>> I looked at the GSL implementation, based on same reference, and >>>> their loop is cleaner. What about porting that implementation here? >>>> Possible? >>>> >>>> My implementation is also using one more term for P than for Q, >>>> which is discouraged in GSL, according to their comments. >>> >>> Ed, do you have any comment about this point? >>> Regarding the 2nd line, after my observations, usually term stops >>> contributing to P, Q before k > nu/2, so actually, an offset by one >>> is most likely without consequences. >>> GSL implementation is nevertheless more elegant. >>> >> I've seen their implementation. It's tight. Feel free to port it. >> I assume this is it: >> int >> gsl_sf_bessel_Jnu_asympx_e(const double nu, const double x, >> gsl_sf_result * result); >> >> You do want to give Q or b one more term so as to make that last >> factor as small as possible. They're right. >> > > Here it is. > > diff --git a/libstdc++-v3/include/tr1/bessel_function.tcc > b/libstdc++-v3/include/tr1/bessel_function.tcc > index 7ac733d..7842350 100644 > --- a/libstdc++-v3/include/tr1/bessel_function.tcc > +++ b/libstdc++-v3/include/tr1/bessel_function.tcc > @@ -353,21 +353,47 @@ namespace tr1 > * @param __x The argument of the Bessel functions. > * @param __Jnu The output Bessel function of the first kind. > * @param __Nnu The output Neumann function (Bessel function > of the second kind). > + * > + * Adapted for libstdc++ from GNU GSL version 2.4 > specfunc/bessel_j.c > + * Copyright (C) 1996,1997,1998,1999,2000,2001,2002,2003 Gerard > Jungman > */ > template <typename _Tp> > 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 __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) < std::abs(__eps * __P); > + __P += __term; > + > + k++; > + > + __term *= (__mu - (2 * k - 1) * (2 * k - 1)) / (k * __8x); > + __epsQ = std::abs(__term) < std::abs(__eps * __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(); > > > >> In principal, these ideas should be ported to I,K but I think (and >> IIRC GSL agrees) these under,over-flow before they have much effect. >> I think I put the full series in there. They could use a similar >> clean-up. >> >> But let's get this in there first. >> >> Ed >> >> > > This looks good to me. The only nit is that you have to uglify the k loop variable. (I hope we can un-uglify our variables when modules come!) Do you have copyright assignment in place and all that? Ed
On 05/01/18 09:54 -0500, Ed Smith-Rowland wrote: >On 01/04/2018 03:54 PM, Michele Pezzutti wrote: >> >> >>On 01/04/2018 06:16 PM, Ed Smith-Rowland wrote: >>>On 01/03/2018 02:49 PM, Michele Pezzutti wrote: >>>> >>>>Hi. >>>> >>>>On 01/02/2018 05:43 PM, Michele Pezzutti wrote: >>>>> >>>>>On 01/02/2018 02:28 AM, Ed Smith-Rowland wrote: >>>>>>I like the patch. >>>>>> >>>>>>I have a similar one in the tr29124 branch. >>>>>> >>>>>>Anyway, I got held up and I think it's good to have new >>>>>>folks looking into this. >>>>>> >>>>>>It looks good except that you need to uglify k. >>>>> >>>>>I looked at the GSL implementation, based on same reference, >>>>>and their loop is cleaner. What about porting that >>>>>implementation here? Possible? >>>>> >>>>>My implementation is also using one more term for P than for >>>>>Q, which is discouraged in GSL, according to their comments. >>>> >>>>Ed, do you have any comment about this point? >>>>Regarding the 2nd line, after my observations, usually term >>>>stops contributing to P, Q before k > nu/2, so actually, an >>>>offset by one is most likely without consequences. >>>>GSL implementation is nevertheless more elegant. >>>> >>>I've seen their implementation. It's tight. Feel free to port it. >>>I assume this is it: >>>int >>>gsl_sf_bessel_Jnu_asympx_e(const double nu, const double x, >>>gsl_sf_result * result); >>> >>>You do want to give Q or b one more term so as to make that last >>>factor as small as possible. They're right. >>> >> >>Here it is. >> >>diff --git a/libstdc++-v3/include/tr1/bessel_function.tcc >>b/libstdc++-v3/include/tr1/bessel_function.tcc >>index 7ac733d..7842350 100644 >>--- a/libstdc++-v3/include/tr1/bessel_function.tcc >>+++ b/libstdc++-v3/include/tr1/bessel_function.tcc >>@@ -353,21 +353,47 @@ namespace tr1 >> * @param __x The argument of the Bessel functions. >> * @param __Jnu The output Bessel function of the first kind. >> * @param __Nnu The output Neumann function (Bessel >>function of the second kind). >>+ * >>+ * Adapted for libstdc++ from GNU GSL version 2.4 >>specfunc/bessel_j.c >>+ * Copyright (C) 1996,1997,1998,1999,2000,2001,2002,2003 >>Gerard Jungman >> */ >> template <typename _Tp> >> 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 __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) < std::abs(__eps * __P); >>+ __P += __term; >>+ >>+ k++; >>+ >>+ __term *= (__mu - (2 * k - 1) * (2 * k - 1)) / (k * __8x); >>+ __epsQ = std::abs(__term) < std::abs(__eps * __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(); >> >> >> >>>In principal, these ideas should be ported to I,K but I think (and >>>IIRC GSL agrees) these under,over-flow before they have much >>>effect. >>>I think I put the full series in there. They could use a similar >>>clean-up. >>> >>>But let's get this in there first. >>> >>>Ed >>> >>> >> >> >This looks good to me. The only nit is that you have to uglify the k >loop variable. Right, it needs to be __k. I'm also not sure we can put the copyright notice there, rather than at the top of the file. We can collapse the years to 1996-2003 so I suggest adding to the top of the file something like: /* __cyl_bessel_jn_asymp adapted from GNU GSL version 2.4 specfunc/bessel_j.c * Copyright (C) 1996-2003 Gerard Jungman */ >(I hope we can un-uglify our variables when modules come!) > >Do you have copyright assignment in place and all that? That appears to be underway. Thanks for working on this.
On 01/05/2018 05:50 PM, Jonathan Wakely wrote: > On 05/01/18 09:54 -0500, Ed Smith-Rowland wrote: >> This looks good to me. The only nit is that you have to uglify the k >> loop variable. > > Right, it needs to be __k. > > I'm also not sure we can put the copyright notice there, rather than > at the top of the file. We can collapse the years to 1996-2003 so I > suggest adding to the top of the file something like: > > /* __cyl_bessel_jn_asymp adapted from GNU GSL version 2.4 > specfunc/bessel_j.c > * Copyright (C) 1996-2003 Gerard Jungman > */ > diff --git a/libstdc++-v3/include/tr1/bessel_function.tcc b/libstdc++-v3/include/tr1/bessel_function.tcc index 7ac733d..6f5ff34 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,16 +362,40 @@ 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) < std::abs(__eps * __P); + __P += __term; + + __k++; + + __term *= (__mu - (2 * __k - 1) * (2 * __k - 1)) / (__k * __8x); + __epsQ = std::abs(__term) < std::abs(__eps * __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(); >> (I hope we can un-uglify our variables when modules come!) >> >> Do you have copyright assignment in place and all that? > > That appears to be underway. No news since I sent questionnaire though. > > Thanks for working on this. > Welcome.
Hi, On 05/01/2018 23:46, Michele Pezzutti wrote: > + __term *= (__k == 0) ? _Tp(1) : -(__mu - (2 * __k - 1) * (2 * __k - 1)) > + / (__k * __8x); In such cases, per the Coding Standards, you want an outer level of parentheses wrapping the whole right side expression. See toward the end of https://www.gnu.org/prep/standards/html_node/Formatting.html#Formatting. I see that many other places will need fixing, at some point - this rather straightforward rule is often overlooked leading to brittle formatting of many expressions, probably because it's really obvious in practice only together with Emacs, I'm not sure. Also - this kind of stylistic nitpicking is partially personal taste - the parentheses around (__k == 0) seem definitely redundant to me, and I don't think you would find many examples of that in our code. About the Assignement, please be patient. For example, used to be the case that when RMS was traveling couldn't process Assignments, for example. It can be slow for many reasons, it's 100% normal. Paolo..
Formatting fixed. diff --git a/libstdc++-v3/include/tr1/bessel_function.tcc b/libstdc++-v3/include/tr1/bessel_function.tcc index 7ac733d..5f8fc9f 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,16 +362,42 @@ 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) < std::abs(__eps * __P); + __P += __term; + + __k++; + + __term *= (__mu - (2 * __k - 1) * (2 * __k - 1)) / (__k * __8x); + __epsQ = std::abs(__term) < std::abs(__eps * __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(); On 01/06/2018 10:23 AM, Paolo Carlini wrote: > Hi, > > On 05/01/2018 23:46, Michele Pezzutti wrote: >> + __term *= (__k == 0) ? _Tp(1) : -(__mu - (2 * __k - 1) * (2 * __k - >> 1)) >> + / (__k * __8x); > In such cases, per the Coding Standards, you want an outer level of > parentheses wrapping the whole right side expression. See toward the > end of > https://www.gnu.org/prep/standards/html_node/Formatting.html#Formatting. > I see that many other places will need fixing, at some point - this > rather straightforward rule is often overlooked leading to brittle > formatting of many expressions, probably because it's really obvious > in practice only together with Emacs, I'm not sure. > > Also - this kind of stylistic nitpicking is partially personal taste - > the parentheses around (__k == 0) seem definitely redundant to me, and > I don't think you would find many examples of that in our code. > > About the Assignement, please be patient. For example, used to be the > case that when RMS was traveling couldn't process Assignments, for > example. It can be slow for many reasons, it's 100% normal. > > Paolo.. > >
On 01/08/2018 02:08 PM, Michele Pezzutti wrote: > Formatting fixed. > > diff --git a/libstdc++-v3/include/tr1/bessel_function.tcc > b/libstdc++-v3/include/tr1/bessel_function.tcc > index 7ac733d..5f8fc9f 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,16 +362,42 @@ 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) < std::abs(__eps * __P); > + __P += __term; > + > + __k++; > + > + __term *= (__mu - (2 * __k - 1) * (2 * __k - 1)) / (__k * > __8x); > + __epsQ = std::abs(__term) < std::abs(__eps * __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(); > > On 01/06/2018 10:23 AM, Paolo Carlini wrote: >> Hi, >> >> On 05/01/2018 23:46, Michele Pezzutti wrote: >>> + __term *= (__k == 0) ? _Tp(1) : -(__mu - (2 * __k - 1) * (2 * __k >>> - 1)) >>> + / (__k * __8x); >> In such cases, per the Coding Standards, you want an outer level of >> parentheses wrapping the whole right side expression. See toward the >> end of >> https://www.gnu.org/prep/standards/html_node/Formatting.html#Formatting. >> I see that many other places will need fixing, at some point - this >> rather straightforward rule is often overlooked leading to brittle >> formatting of many expressions, probably because it's really obvious >> in practice only together with Emacs, I'm not sure. >> >> Also - this kind of stylistic nitpicking is partially personal taste >> - the parentheses around (__k == 0) seem definitely redundant to me, >> and I don't think you would find many examples of that in our code. >> >> About the Assignement, please be patient. For example, used to be the >> case that when RMS was traveling couldn't process Assignments, for >> example. It can be slow for many reasons, it's 100% normal. >> >> Paolo.. >> >> > > Hi all, I'd like us to get this in. Michele, how is the Copyright assignment going? I'd like to push it down through brach/gcc-6 also. Ed
Hi. The assignment/disclaimer process with the FSF is currently complete. Michele Il 06.05.2018 03:34 Ed Smith-Rowland ha scritto: > On 01/08/2018 02:08 PM, Michele Pezzutti wrote: > >> Formatting fixed. diff --git a/libstdc++-v3/include/tr1/bessel_function.tcc b/libstdc++-v3/include/tr1/bessel_function.tcc index 7ac733d..5f8fc9f 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,16 +362,42 @@ 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::epsilon(); + + do + { + __term *= (__k == 0 + ? _Tp(1) + : -(__mu - (2 * __k - 1) * (2 * __k - 1)) / (__k * __8x)); + + __epsP = std::abs(__term) < std::abs(__eps * __P); + __P += __term; + + __k++; + + __term *= (__mu - (2 * __k - 1) * (2 * __k - 1)) / (__k * __8x); + __epsQ = std::abs(__term) < std::abs(__eps * __Q); + __Q += __term; + + if (__epsP && __epsQ && __k > __nu / 2.) + break; + + __k++; + } + while (__k < 1000); + const _Tp __chi = __x - (__nu + _Tp(0.5L)) * __numeric_constants::__pi_2(); On 01/06/2018 10:23 AM, Paolo Carlini wrote: >> >>> Hi, On 05/01/2018 23:46, Michele Pezzutti wrote: >>> >>>> + __term *= (__k == 0) ? _Tp(1) : -(__mu - (2 * __k - 1) * (2 * __k - 1)) + / (__k * __8x); >>> In such cases, per the Coding Standards, you want an outer level of parentheses wrapping the whole right side expression. See toward the end of https://www.gnu.org/prep/standards/html_node/Formatting.html#Formatting [1]. I see that many other places will need fixing, at some point - this rather straightforward rule is often overlooked leading to brittle formatting of many expressions, probably because it's really obvious in practice only together with Emacs, I'm not sure. Also - this kind of stylistic nitpicking is partially personal taste - the parentheses around (__k == 0) seem definitely redundant to me, and I don't think you would find many examples of that in our code. About the Assignement, please be patient. For example, used to be the case that when RMS was traveling couldn't process Assignments, for example. It can be slow for many reasons, it's 100% normal. Paolo.. > > Hi all, > > I'd like us to get this in. > > Michele, how is the Copyright assignment going? > > I'd like to push it down through brach/gcc-6 also. > > Ed Con Mobile Open 6 GB hai 6 Giga, 600 minuti e 300 SMS per il tuo smartphone a 9€ al mese per sempre. Passa ora a Tiscali Mobile, il nostro mese è vero! http://tisca.li/Open6GB0318
Michele, So I think the patch you had plus tests testsuite/tr1/5_numerical_facilities/special_functions/09_cyl_bessel_j/check_value.cc testsuite/tr1/5_numerical_facilities/special_functions/11_cyl_neumann/check_value.cc plus I'd like the new tests also in testsuite/special_functions/08_cyl_bessel_j/check_value.cc testsuite/special_functions/10_cyl_neumann/check_value.cc for C++17. Put that up with a Changelog and see if we're all good. Thank you for going through this and other math issues. Ed
Hi. This patch intends to fix Bug 83566 - cyl_bessel_j returns wrong result for x>1000 for high orders. See https://gcc.gnu.org/bugzilla/show_bug.cgi?id=83566 for issue description. * libstdc++-v3/include/tr1/bessel_function.tcc Series expansion in __cyl_bessel_jn_asymp() shall not be truncated at the first terms. At least nu/2 terms shall be added, in order to have a guaranteed error bound. * libstdc++-v3/testsuite/tr1/5_numerical_facilities/special_functions/09_cyl_bessel_j/check_value.cc Testcase for x > 1000 added. * libstdc++-v3/testsuite/tr1/5_numerical_facilities/special_functions/11_cyl_neumann/check_value.cc Testcase for x > 1000 added. * libstdc++-v3/testsuite/special_functions/08_cyl_bessel_j/check_value.cc Testcase for x > 1000 added. * libstdc++-v3/testsuite/special_functions/10_cyl_neumann/check_value.cc Testcase for x > 1000 added. 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<double> +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<typename Ret, unsigned int Num> void test(const testcase_cyl_bessel_j<Ret> (&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<double> +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<typename Ret, unsigned int Num> void test(const testcase_cyl_neumann<Ret> (&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<double> +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<typename Ret, unsigned int Num> void test(const testcase_cyl_bessel_j<Ret> (&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<double> +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<typename Ret, unsigned int Num> void test(const testcase_cyl_neumann<Ret> (&data)[Num], Ret toler) @@ -794,5 +827,6 @@ main() test(data026, toler026); test(data027, toler027); test(data028, toler028); + test(data029, toler029); return 0; } Regards Michele On 06/05/2018 17:56, Ed Smith-Rowland wrote: > Michele, > > So I think the patch you had plus tests > testsuite/tr1/5_numerical_facilities/special_functions/09_cyl_bessel_j/check_value.cc > > testsuite/tr1/5_numerical_facilities/special_functions/11_cyl_neumann/check_value.cc > > plus I'd like the new tests also in > testsuite/special_functions/08_cyl_bessel_j/check_value.cc > testsuite/special_functions/10_cyl_neumann/check_value.cc > for C++17. > > Put that up with a Changelog and see if we're all good. > > Thank you for going through this and other math issues. > > Ed >
diff --git a/libstdc++-v3/include/tr1/bessel_function.tcc b/libstdc++-v3/include/tr1/bessel_function.tcc index 7ac733d..852a31e 100644 --- a/libstdc++-v3/include/tr1/bessel_function.tcc +++ b/libstdc++-v3/include/tr1/bessel_function.tcc @@ -359,15 +359,32 @@ namespace tr1 __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 __8x = _Tp(8) * __x; + + _Tp __P = _Tp(1); + _Tp __Q = _Tp(0); + + _Tp __eps = std::numeric_limits<_Tp>::epsilon(); + + _Tp __term = _Tp(1); + _Tp __epsP = _Tp(0); + _Tp __epsQ = _Tp(0); + + unsigned long __2k; + for (unsigned long k = 1; k < 1000; k+=2) { + __2k = 2 * k; + + __term *= (__mu - _Tp((__2k - 1) * (__2k - 1))) / (k * __8x); + __Q += __term; + __epsQ = std::abs(__term) < std::abs(__eps * __Q); + + __term *= -(__mu - _Tp((__2k + 1) * (__2k + 1)))/ (_Tp(k + 1) * __8x); + __P += __term; + __epsP = std::abs(__term) < std::abs(__eps * __P); + + if (__epsP && __epsQ && k > __nu / 2.) + break; + } const _Tp __chi = __x - (__nu + _Tp(0.5L))