bugGNU Octave - Bugs: bug #48316, besselj gives NaN for moderately...

 
 

bug #48316: besselj gives NaN for moderately large inputs (1e10)

Submitter:  Colin Macdonald <cbm>
Submitted:  Sun 26 Jun 2016 05:55:18 AM UTC
   
 
Category:  Octave Function Severity:  3 - Normal
Priority:  5 - Normal Item Group:  Incorrect Result
Status:  Fixed Assigned to:  None
Originator Name:  Open/Closed:  * Closed
Release:  * dev Operating System:  * Any
Fixed Release:  None Planned Release:  None
* Mandatory Fields

Add a New Comment Rich Markup
   

Tue 27 Mar 2018 12:18:55 AM UTC, comment #3: 

Fixes from GSOC2018 special function project were incorporated during OctConf 2018.  Closing report.

Rik <rik5>
Group administrator
Thu 30 Jun 2016 07:44:31 PM UTC, comment #2: 

Scipy appears to use a combination of the amos library and the cephes bessel functions (http://www.netlib.org/cephes/).

Mike Miller <mtmiller>
Group Member
Thu 30 Jun 2016 07:35:48 PM UTC, comment #1: 

Changing status to confirmed, and adding some details from the mailing list (http://thread.gmane.org/gmane.comp.gnu.octave.maintainers/46617).

The bessel functions come from the amos library (http://netlib.org/amos/) which is embedded in liboctave, so any numerical problems are inherited from that implementation.

Mike Miller <mtmiller>
Group Member
Sun 26 Jun 2016 05:55:18 AM UTC, original submission:  

Octave gives "besselj(1, 1e10)" as "NaN + NaNI".

Octave:

>> O = besselj(1, 10.^[5:12]')
O =
   0.00184675756288257 +                    0i
  -0.000725968356813763 +                    0i
  0.000236899240312058 +                    0i
  7.30639118155185e-05 +                    0i
  -5.21042264155388e-06 +                    0i
                 NaN +               NaNi
                 NaN +               NaNi
                 NaN +               NaNi


Symbolic pkg, converted to double:

>> S = double(besselj(1, sym(10).^[5:12]'))
S =
   0.00184675756288257
  -0.000725968356813763
  0.000236899240312058
  7.30639118155185e-05
  -5.21042264155388e-06
  -7.67650817568416e-06
  9.95266266795786e-07
  -7.91380268385044e-07


The relative accuracy is fine before blow-up:

>> abs(O-S)./abs(S)
ans =
  1.17416838493212e-16
  1.49345651543824e-16
  1.14415961302507e-16
                     0
  1.62565113336314e-16
                                NaN
                                NaN
                                NaN


(Ironically, I found this when I went looking for a function with good relative error behaviour even though it decays to zero!  viz., #48307)

Colin Macdonald <cbm>

 

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

Attach Files:
   
   
Comment:
   

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 rik5 (Posted a comment)
  • -email is unavailable- added by mtmiller (Posted a comment)
  • -email is unavailable- added by cbm (Submitted the item)
  • -email is unavailable- added by cbm
  •  

    There are 0 votes so far. Votes easily highlight which items people would like to see resolved in priority, independently of the priority of the item set by tracker managers.

    Only group members can vote.

     

    Follow 5 latest changes.

    Date Changed by Updated Field Previous Value => Replaced by
    2018-03-27 rik5 StatusConfirmed Fixed
        Open/ClosedOpen Closed
    2016-11-16 mtmiller Operating SystemOther Any
    2016-06-30 mtmiller StatusNone Confirmed
    2016-06-26 cbm Carbon-Copy- Added sebald

    Back to the top

    Powered by Savane 3.13-cf05.
    Corresponding source code