diff mbox series

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

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

Commit Message

Michele Pezzutti Jan. 1, 2018, 2:38 a.m. UTC 
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.


                              * __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;
  }

Comments

Ed Smith-Rowland Jan. 2, 2018, 1:28 a.m. UTC | #1 
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
Ed Smith-Rowland Jan. 2, 2018, 1:44 a.m. UTC | #2 
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
Michele Pezzutti Jan. 2, 2018, 4:43 p.m. UTC | #3 
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.
Ed Smith-Rowland Jan. 2, 2018, 4:59 p.m. UTC | #4 
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!
Michele Pezzutti Jan. 2, 2018, 9:41 p.m. UTC | #5 
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!
>
>
Ed Smith-Rowland Jan. 3, 2018, 5:34 a.m. UTC | #6 
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
Michele Pezzutti Jan. 3, 2018, 7:49 p.m. UTC | #7 
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.
Ed Smith-Rowland Jan. 4, 2018, 5:16 p.m. UTC | #8 
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
Michele Pezzutti Jan. 4, 2018, 8:54 p.m. UTC | #9 
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
>
>
Michele Pezzutti Jan. 4, 2018, 9:08 p.m. UTC | #10 
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
>>
>>
>
Ed Smith-Rowland Jan. 5, 2018, 2:54 p.m. UTC | #11 
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
Jonathan Wakely Jan. 5, 2018, 4:50 p.m. UTC | #12 
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.
Michele Pezzutti Jan. 5, 2018, 10:46 p.m. UTC | #13 
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.
Paolo Carlini Jan. 6, 2018, 9:23 a.m. UTC | #14 
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..
Michele Pezzutti Jan. 8, 2018, 7:08 p.m. UTC | #15 
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..
>
>
Ed Smith-Rowland May 6, 2018, 1:34 a.m. UTC | #16 
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
Michele Pezzutti May 6, 2018, 10:41 a.m. UTC | #17 
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
Ed Smith-Rowland May 6, 2018, 3:56 p.m. UTC | #18 
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
Michele Pezzutti May 20, 2018, 8:49 p.m. UTC | #19 
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 mbox series

Patch

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))