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 authentification level.

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

Submitted by:  Rhys Ulerich <rhysu>
Submitted on:  Tue 28 Aug 2012 02:15:22 PM UTC  
 
Category: Accuracy problemSeverity: 3 - Normal
Operating System: Status: None
Assigned to: NoneOpen/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>
Project 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>
Project Administrator
Tue 28 Aug 2012 02:15:22 PM UTC, original submission:

From: -unavailable- Jakub Benda
To: -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>
Project Administrator

 

No files currently attached

 

Depends on the following items: None found

Items that depend on this one: None found

 

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

    Please enter the title of George Orwell's famous dystopian book (it's a date):

     

     

    Follow 2 latest changes.

    Date Changed By Updated Field Previous Value => Replaced By
    Tue 28 Aug 2012 02:15:22 PM UTCrhysuCarbon-Copy-=>Added -unavailable-
      Carbon-Copy-=>Added -unavailable-

    Back to the top


    Powered by Savane 3.1-cleanup