Bug 15529 - Bessel functions less accurate than they should be
Summary: Bessel functions less accurate than they should be
Status: CLOSED FIXED
Alias: None
Product: R
Classification: Unclassified
Component: Accuracy (show other bugs)
Version: R 3.0.1
Hardware: x86_64/x64/amd64 (64-bit) All
: P5 minor
Assignee: R-core
URL:
Depends on:
Blocks:
 
Reported: 2013-11-04 14:25 UTC by M Welinder
Modified: 2015-12-14 13:46 UTC (History)
1 user (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-04 14:25:46 UTC
> options(digits=16)
> besselI(2.125,-5+1/1024)
[1] 0.02679209380096113

There more correct value is

    0.0267920938009571

(From Mathematica -- note that order of arguments differ there)

Luckily it is easy to fix.  The bessel functions' implementations contain
fragments like

   sin(-M_PI * alpha)
   sin(M_PI * alpha)
   cos(M_PI * alpha)

The least amount of effort is to replace "alpha" in the above expressions
by "fmod(alpha,2)".  Slightly fancier would be to create sin_pi and cos_pi
nominally computing sin(x*pi), but reducing the argument to 0.5 or less
before scaling by pi.
Comment 1 Martin Maechler 2013-11-06 13:37:39 UTC
I like the idea of solving the problem more generally, via
more generally useful  sin_pi(), cos_pi() [and tan_pi() .. I need in "my" copula package]   which will even be part of the API (Rmath.h).

I won't get to this for several days though.
Comment 2 Martin Maechler 2013-11-08 09:00:24 UTC
I've introduced cospi(), sinpi() and tanpi(),
committed to  R-devel  [svn rev 64171].

This even improves lgamma() very slightly.

{"cospi" and not "cos_pi"  as I see these in C-like language documents such as
  OpenCL or  "CUDA C"
Comment 3 M Welinder 2013-11-08 14:53:03 UTC
Your cospi looks fine, but sinpi and tanpi do not.  Specifically, for
very small negative arguments, the following line will kill accuracy:

    if(x < 0) x += 2.;

I suggest do this instead:

    if(x <= -1) x += 2.; else if (x > 1) x -= 2;

which will bring you into ]-1;+1].
Comment 4 Martin Maechler 2013-11-11 08:00:58 UTC
(In reply to comment #3)
> Your cospi looks fine, but sinpi and tanpi do not.  Specifically, for
> very small negative arguments, the following line will kill accuracy:
> 
>     if(x < 0) x += 2.;
> 
> I suggest do this instead:
> 
>     if(x <= -1) x += 2.; else if (x > 1) x -= 2;
> 
> which will bring you into ]-1;+1].

Morten, you are entirely correct... indeed something I should not have missed.
Thank you!