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

[Ed Smith-Rowland]

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