# GNU Scientific Library - Bugs: bug #51240, Bessel functions for nu < 0

You are not allowed to post comments on this tracker with your current authentication level.

# bug #51240: Bessel functions for nu < 0

 Submitted by: Patrick Alken Submitted on: Thu 15 Jun 2017 09:59:56 AM UTC Category: Runtime error Severity: 3 - Normal Operating System: Status: Fixed Assigned to: None Open/Closed: Closed Release:

## Discussion

Thu 07 Sep 2017 08:42:39 AM UTC, comment #3:

Patrick Alken <psa>
Mon 26 Jun 2017 07:46:55 AM UTC, comment #2:

After Gerards suggestion I've implemented (see attached) the routines sin_pi and cos_pi. Essentially, modf and fmod do the trick. I'm a bit surprised at how simple the solution is, so please check if you see any problems, but the numerical tests all pass.

Best,

Patrick Alken <psa>
Fri 16 Jun 2017 05:29:07 AM UTC, comment #1:

From Gerard:

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

I do not remember much about what I was thinking 20 years ago.
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

Patrick Alken <psa>
Thu 15 Jun 2017 09:59:56 AM UTC, original submission:

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. We would only need to
replace (in bessel_Ynu.c) the block

if(x <= 0.0 || nu < 0.0) {
DOMAIN_ERROR(result);
}

by something like

if(x <= 0.0) {
DOMAIN_ERROR(result);
}
else if (nu < 0.0) {
int Jstatus = gsl_sf_bessel_Jnu_e(-nu, x, result);
double Jval = result->val;
double Jerr = result->err;
int Ystatus = gsl_sf_bessel_Ynu_e(-nu, x, result);
double Yval = result->val;
double Yerr = result->err;
double s = sin(M_PI*nu), c = cos(M_PI*nu);
result->val = c*Yval - s*Jval;
result->err = fabs(c*Yerr) + fabs(s*Jerr);
return (Jstatus != GSL_SUCCESS ? Jstatus : Ystatus);
}

What do you think? I've attached a possible patch for the Bessel
functions Y_nu and also J_nu. It compiles without problems, but I
haven't written tests for this case yet, as I first wanted to see why
the case nu < 0 was originally excluded.

Patrick Alken <psa>

## Attached Files

Attached Files
file #41048:  sincos_pi.patch added by psa (6KiB - text/x-patch)
file #41049:  gsl_sf_sincos_pi.h added by psa (2KiB - text/x-chdr)
file #41050:  sincos_pi.c added by psa (2KiB - text/x-csrc)
file #41051:  test_sincos_pi.c added by psa (5KiB - text/x-csrc)
file #40920:  bessel.patch added by psa (2KiB - text/x-patch)

## Dependencies

Depends on the following items: None found

Items that depend on this one: None found

Carbon-Copy List
• -email is unavailable- added by psa (Submitted the item)
•