Re: [Bug-gsl] [bug #51240] Bessel functions for nu < 0

[Konrad]

Dear Gerard and Patrick,

thanks for your insightful answers, it is indeed quite tempting to take
the standard trigonometric functions for granted. If I find the time
maybe as a first step separate functions sin_pi and cos_pi with proper
argument reduction could be added and once this is implemented we can
come back to this issue.

Cheers,
Konrad



On 06/16/2017 07:30 AM, Patrick Alken wrote:
> Thanks Gerard, its good to know you've thought about this already. I've
> added your comments to the bug tracker for this issue.
> 
> On 06/15/2017 09:45 PM, Gerard Jungman wrote:
>> On 06/15/2017 03:59 AM, Patrick Alken wrote:
>>>
>>> recently I needed some of the Bessel functions Y_nu, but noticed that
>>> the GSL-implementation raises a domain error for nu < 0. I'm not sure
>>> why this is the case, maybe there are good reasons I have overlooked,
>>> but to me it seems we could use the rotation between Bessel functions of
>>> first and second kind to cover all values of nu. 
>>
>>
>> There are three issues. One is mathematical, one is a related
>> observation involving the mirror implementation for J, and the
>> last is a smallish but important design issue related to the
>> form of the proposed change.
>>
>> First the math issue. Suppose nu is exactly integral. The naive
>> evaluation of sin(M_PI*nu) and cos(M_PI*nu) could be very bad,
>> since you really want these to be exactly 0.0 and 1.0, to
>> avoid strange admixtures of possibly singular values.
>> The behavior near integral nu is also suspect, since
>> you must validate the behavior of the trig functions
>> near zeros.
>>
>> Therefore, you really need a set of functions which are sometimes
>> called "sin_pi(x)" and "cos_pi(x)" which perform the argument
>> reduction of x exactly. The built-in argument reduction of M_PI*x
>> is not acceptable.
>>
>> Now for the observation involving J, which is directly related
>> to this point. It may be that admixtures of the near-singular
>> solution are ok when evaluating Y, since Y dominates any
>> solution, but that the problem is much more difficult in
>> the case of J.
>>
>> The difference may be large enough that you can imagine
>> easily implementing Y but punting on J. Then you would
>> be left with two functions which are mathematically
>> similar but for which the exported interfaces have
>> different preconditions. This is confusing enough
>> that it becomes an argument for leaving them both out.
>>
>> Now for the smallish design issue.
>>
>> If these problems can be solved in a coherent way, then
>> it would be fine to provide Y and J for negative nu by
>> rotation. But the code needs to be refactored. The proposed
>> change involves the function calling itself, with flipped argument,
>> which is bad levelization. Rather, we need separate functions
>> which handle positive nu, which are called by wrappers which
>> handle both positive and negative nu.
>>
>> The main task of the wrappers would be to manage the
>> delegation to properly argument-reduced sin_pi and
>> cos_pi functions when nu < 0, so their job is not
>> entirely trivial after all.
>>
>> Whether or not we continue to export the nu>0 interfaces
>> is a separate question. I don't think I have any special
>> feelings about that.
>>
>>
>> I do not remember much about what I was thinking 20 years ago.
>> But I can guess that the main problem was this business of
>> proper argument reduction. Maybe it is a red herring after
>> all, but this would need to be demonstrated.
>>
>> As a general comment, it is always necessary to be very careful
>> about the quality of the underlying trig functions in this kind
>> of situation. Sadly, quality will vary from platform to platform.
>>
>> -- 
>> G. Jungman
> 
>