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

[Jonathan Wakely]

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.