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

[Ed Smith-Rowland <3dw4rd at verizon dot net>]

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