search

MichaelChirico / r-bugs

A ⚠️read-only⚠️mirror of https://bugs.r-project.org/
20 stars 0 forks source link

[BUGZILLA #15554] bessel function issues #5129

Open MichaelChirico opened 4 years ago

MichaelChirico commented 4 years ago

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.

METADATA

MichaelChirico commented 4 years ago

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.

METADATA

MichaelChirico commented 4 years ago

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.

METADATA

MichaelChirico commented 4 years ago

Notes to myself (:-)

METADATA

MichaelChirico commented 4 years ago

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.

METADATA

MichaelChirico commented 4 years ago

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

METADATA

MichaelChirico commented 4 years ago

(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.

METADATA

MichaelChirico commented 4 years ago

(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

METADATA

MichaelChirico commented 4 years ago

This thread says that a fix for issue 2 was committed but running the test code I still get the same results as in the original reports (for both 1 and 2).

METADATA