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

[Gerard Jungman]

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