bugGNU Scientific Library - Bugs: bug #37209, gsl_sf_bessel_jl_e returns NaN for...

 
 

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

bug #37209: gsl_sf_bessel_jl_e returns NaN for large inputs without setting error code

Submitter:  Rhys Ulerich <rhysu>
Submitted:  Tue 28 Aug 2012 02:15:22 PM UTC
   
 
Category:  Accuracy problem Severity:  3 - Normal
Operating System:  Status:  None
Assigned to:  None Open/Closed:  Open
Release: 

Fri 04 Oct 2013 05:39:56 PM UTC, comment #2: 

The two issues below are in fact separate bugs in the spherical bessel code.

The bug report from 28 Aug 2012 involves an instability condition in the downward recurrence relation after computing the CF1 continued fraction. This instability appears when x is small relative to l. There is a short discussion of this in Numerical Recipes, but I have been unable to find a fully acceptable solution in the literature. Many other bessel codes are susceptable to this same bug (including NR and Barnett's SBESJY.F routine).

I was able to fix this particular test case by starting with a smaller initial value to the recurrence relation, which delays (but does not completely fix) the overflow condition. The code now works for l = 364, x = 36.62, however there are probably other combinations of large l, small x which will cause an overflow in the downward recurrence.

This is partially fixed in rev 4850.

------

The bug report from 23 Sep 2013 involves a different branch of the code and the gsl_sf_bessel_IJ_taylor_e function - this is still open.

Patrick Alken <psa>
Group administrator
Mon 23 Sep 2013 03:46:15 PM UTC, comment #1: 

A related bug report from JE Campagne (campagne@lal.in2p3.fr):

Dear developpers

I am running GSL 1.16 with such architecture

Linux opera3.lal.in2p3.fr 2.6.32-220.13.1.el6.x86_64 #1 SMP Tue Apr 17
15:16:22 CDT 2012 x86_64 x86_64 x86_64 GNU/Linux

and  gcc compilor version gcc (GCC) 4.4.5 20110214

I have a problem both with
gsl_sf_bessel_jl_e(int l, double x, gsl_sf_result* r)
or
gsl_sf_bessel_jl (int l, double x)

routine that fails to compute with l>~256 and produce such message

#include <gsl/gsl_sf_bessel.h>
#include <gsl/gsl_errno.h>
#include <stdlib.h>

int main(int argc, char* argv[]){

   If xval=10 it crashes for l=244
   but if xval=1 it crashes for l=149

   double xval = 1.0;
   int lmin=100;
   int lmax=200;

   gsl_sf_result result;

   for( int l=lmin;l<=lmax;l++){
     int status = gsl_sf_bessel_jl_e(l,xval,&result);
     printf("order[%d]: %s",l,gsl_strerror(status));
     if ( status == 0 ) printf("\t result: %le",result.val);
     printf("\n");
   }
}

and the execution output
order[100]: success      result: 7.444728e-190
order[101]: success      result: 3.667442e-192
order[102]: success      result: 1.789038e-194
order[103]: success      result: 8.642896e-197
order[104]: success      result: 4.135451e-199
order[105]: success      result: 1.959973e-201
order[106]: success      result: 9.201951e-204
order[107]: success      result: 4.280069e-206
order[108]: success      result: 1.972424e-208
order[109]: success      result: 9.006687e-211
order[110]: success      result: 4.075506e-213
order[111]: success      result: 1.827618e-215
order[112]: success      result: 8.122905e-218
order[113]: success      result: 3.578441e-220
order[114]: success      result: 1.562668e-222
order[115]: success      result: 6.764920e-225
order[116]: success      result: 2.903452e-227
order[117]: success      result: 1.235534e-229
order[118]: success      result: 5.213315e-232
order[119]: success      result: 2.181341e-234
order[120]: success      result: 9.051363e-237
order[121]: success      result: 3.724903e-239
order[122]: success      result: 1.520394e-241
order[123]: success      result: 6.155540e-244
order[124]: success      result: 2.472144e-246
order[125]: success      result: 9.849335e-249
order[126]: success      result: 3.893078e-251
order[127]: success      result: 1.526721e-253
order[128]: success      result: 5.940636e-256
order[129]: success      result: 2.293716e-258
order[130]: success      result: 8.788311e-261
order[131]: success      result: 3.341611e-263
order[132]: success      result: 1.261003e-265
order[133]: success      result: 4.722924e-268
order[134]: success      result: 1.755758e-270
order[135]: success      result: 6.478899e-273
order[136]: success      result: 2.373255e-275
order[137]: success      result: 8.630131e-278
order[138]: success      result: 3.115611e-280
order[139]: success      result: 1.116721e-282
order[140]: success      result: 3.974145e-285
order[141]: success      result: 1.404309e-287
order[142]: success      result: 4.927460e-290
order[143]: success      result: 1.716906e-292
order[144]: success      result: 5.940921e-295
order[145]: success      result: 2.041577e-297
order[146]: success      result: 6.967922e-300
order[147]: success      result: 2.362034e-302
order[148]: success      result: 7.953067e-305
gsl: gamma.c:1467: ERROR: underflow
Default GSL error handler invoked.
Abort (core dumped)

Patrick Alken <psa>
Group administrator
Tue 28 Aug 2012 02:15:22 PM UTC, original submission:  

From: -email is unavailable- Jakub Benda
To:  -email is unavailable-
Date: Tue, 28 Aug 2012 01:33:59 -0500
Subject: [Bug-gsl] Spherical Bessel function failure
 
The spherical Bessel function returns NaN for some (large) inputs without
setting the error code.
Example program:

--------------------------------------------------------------------------------
#include <stdio.h>
#include <gsl/gsl_sf.h>

int main(void)
{
gsl_sf_result res;
int err = gsl_sf_bessel_jl_e(364, 36.62, &res);

printf("err = %d\n", err);
printf("res.val = %g\n", res.val);
printf("res.err = %g\n", res.err);

return 0;
}
--------------------------------------------------------------------------------

Output:

--------------------------------------------------------------------------------
err = 0
res.val = -nan
res.err = nan
--------------------------------------------------------------------------------

Jakub

Rhys Ulerich <rhysu>
Group administrator

 

(Note: upload size limit is set to 16384 kB, after insertion of the required escape characters.)

No files currently attached

 

Depends on the following items: None found

Items that depend on this one: None found

 

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

    Follow 2 latest changes.

    Date Changed by Updated Field Previous Value => Replaced by
    2012-08-28 rhysu Carbon-Copy- Added -email is unavailable-
        Carbon-Copy- Added -email is unavailable-

    Back to the top

    Powered by Savane 3.13-cf05.
    Corresponding source code