Fix Bug 83566 - cyl_bessel_j returns wrong result for x>1000 for high orders

Michele Pezzutti mpezz@tiscali.it
Sun May 20 20:50:00 GMT 2018 

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
>

More information about the Libstdc++ mailing list