bugGNU 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 <psa>
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: 

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

Konrad's patch added to git

Patrick Alken <psa>
Project Administrator
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,
Konrad

(file #41048, file #41049, file #41050, file #41051)

Patrick Alken <psa>
Project Administrator
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
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

Patrick Alken <psa>
Project Administrator
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>
Project Administrator

 

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)

 

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)
  •  

     

     

    Follow 7 latest changes.

    Date Changed by Updated Field Previous Value => Replaced by
    2017-09-07 psa StatusNone => Fixed
        Open/ClosedOpen => Closed
    2017-06-26 psa Attached File- => Added sincos_pi.patch, #41048
        Attached File- => Added gsl_sf_sincos_pi.h, #41049
        Attached File- => Added sincos_pi.c, #41050
        Attached File- => Added test_sincos_pi.c, #41051
    2017-06-15 psa Attached File- => Added bessel.patch, #40920

    Back to the top


    Powered by Savane 3.5