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

[Michele Pezzutti]

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.