This is the mail archive of the libstdc++@gcc.gnu.org mailing list for the libstdc++ project.


Index Nav: [Date Index] [Subject Index] [Author Index] [Thread Index]
Message Nav: [Date Prev] [Date Next] [Thread Prev] [Thread Next]
Other format: [Raw text]

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


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!



Index Nav: [Date Index] [Subject Index] [Author Index] [Thread Index]
Message Nav: [Date Prev] [Date Next] [Thread Prev] [Thread Next]