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

[Michele Pezzutti]

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