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

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.

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"

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

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