Bug 15554 - bessel function issues
Summary: bessel function issues
Status: ASSIGNED
Alias: None
Product: R
Classification: Unclassified
Component: Accuracy (show other bugs)
Version: R 3.0.1
Hardware: Other All
: P3 major
Assignee: R-core
URL:
Depends on:
Blocks:
 
Reported: 2013-11-12 16:27 UTC by M Welinder
Modified: 2015-09-10 12:33 UTC (History)
4 users (show)

See Also:


Attachments

Note You need to log in before you can comment on or make changes to this bug.
Description M Welinder 2013-11-12 16:27:30 UTC
Issue 1:

> besselJ(1.5,1700.5)
[0]
Warning message:
In besselJ(1.5, 1700.5) : bessel_j(1.5,nu=1700.5): precision lost in result

In this area (say, x*x/alpha < 0.1), the taylor series for bessel_j
(see 9.1.10 at http://people.math.sfu.ca/~cbm/aands/page_360.htm)
converges rapidly, so there is no need for using the expensive
recurrence formulas.  Same issue with besselJ(3,181) where the
expected result is actually non-zero.



Issue 2:

> besselJ(1.5,-1700.5)
[1] Inf
Warning messages:
1: In besselJ(1.5, -1700.5) :
  bessel_j(1.5,nu=1700.5): precision lost in result
2: In besselJ(1.5, -1700.5) :
  bessel_y(1.5,nu=1700.5): precision lost in result

In this case, cospi(-1700.5) is zero so there is no need to call
bessel_j(1.5,nu=1700.5) in this case.  Computing cospi in this case
is cheap.


Issue 3:

This cast here means that the code will do really strange things for really
large alpha.  2^61-1, 2^63-1, and 2^64 come to mind.

    nb = 1 + (long)na; /* nb-1 <= alpha < nb */

I suggest putting in an explicit limit of, say, 2^24.  At that size, the code
will allocate a table of ~128MB.
Comment 1 Martin Maechler 2013-11-13 10:33:15 UTC
Thank you for the issues!
Note that they (notably 1 and 2) will have to wait for a while,
unless you can provide  (tested) patches.
Comment 2 M Welinder 2013-11-20 16:15:42 UTC
I saw you were looking into replacing the current implementation of bessel
functions a few years back.

There's new (2007) research out there:

    http://www.carma.newcastle.edu.au/jon/bessel_arc.pdf

On the up side, it handles complex arguments and order, and promises as
much accuracy as you want.  And it ought to be fairly fast.

On the down side, it looks like quite a lot of work.

I am going to ask the authors for a reference implementation.  The only
code out there implementing some of this seems to be

    http://trac.common-lisp.net/oct/browser/qd-bessel.lisp

Note, that the lisp comments talk about a typo in the paper.
Comment 3 Martin Maechler 2014-12-29 21:52:28 UTC
Notes to myself (:-)
- My own package 'Bessel' (on R-forge and CRAN) provides several asymptotic formulas, not so many for J() yet.

- Package 'gsl' - provides R functions calling the GSL (GNU Scientfic Library) routines; it uses quite a few different asymptotic formulas (each a separate function) though the separate formulas are not "exported" from C and not available from R.

- Instead of Abramowitz and Stegun ("AaS" or "AandS"), we should really start referencing its designated successor which is completely hyperlinked available on the web: The Digital Library of Mathematical Functions (DLMF), see
http://dlmf.nist.org/  notably the Bessel chapter 10, http://dlmf.nist.org/10
e.g. with the asymptotic expansions such
- for J_nu() and Y_nu(): http://dlmf.nist.org/10.17  and  10.19
- for I_nu() and K_nu(): http://dlmf.nist.org/10.40  and  10.41
Comment 4 M Welinder 2014-12-30 00:23:30 UTC
For the record, I did ask for a reference implementation, but there
isn't one.

A and S (the book) does have its good sides: it doesn't experience
DNS problems (like the http://dlmf.nist.org server seems to) and
there's something nice about a good, physical book.  Even if the
log tables are a bit pointless in this day and age.
Comment 5 Bill Dunlap 2015-02-21 00:35:42 UTC
M Welinder wrote:
   This cast here means that the code will do really strange things for really
   large alpha.  2^61-1, 2^63-1, and 2^64 come to mind.

"Really strange things" includes crashing the R session (with R-3.1.2):
> besselJ(1, 2^64) # or nu=Inf

 *** caught segfault ***
address 0xfffffffc022e73e0, cause 'memory not mapped'

Traceback:
 1: besselJ(1, 2^64)

Possible actions:
1: abort (with core dump, if enabled)
2: normal R exit
3: exit R without saving workspace
4: exit R saving workspace
Comment 6 Martin Maechler 2015-02-23 11:33:14 UTC
(In reply to Bill Dunlap from comment #5)
> M Welinder wrote:
>    This cast here means that the code will do really strange things for
> really
>    large alpha.  2^61-1, 2^63-1, and 2^64 come to mind.
> 
> "Really strange things" includes crashing the R session (with R-3.1.2):
> > besselJ(1, 2^64) # or nu=Inf
> 
>  *** caught segfault ***
> address 0xfffffffc022e73e0, cause 'memory not mapped'
> 
> Traceback:
>  1: besselJ(1, 2^64)
> 
> Possible actions:
> 1: abort (with core dump, if enabled)
> 2: normal R exit
> 3: exit R without saving workspace
> 4: exit R saving workspace

Thank you, Bill, and Morten.

I'm about to commit a patch - along Morten's suggestion for issue 3 -
to both R-patched and R-devel, i.e., the fix will be in R >= 3.1.3.

I cannot close the bug of course, as this only fixes Issue 3.
Comment 7 Martin Maechler 2015-02-24 07:45:08 UTC
(In reply to Martin Maechler from comment #6)
.....................
.....................
> I cannot close the bug of course, as this only fixes Issue 3.

I've commited a patch for  Issue 2 - to R-devel  (svn rev 67886; the log entry wrongly mentions issue 1 instead of 2),  as well.

Issue 1 is open {I've started to explore the Taylor series validity in my package 'Bessel' https://r-forge.r-project.org/R/?group_id=611