29 #ifndef MATHLIB_STANDALONE 30 #include <R_ext/Memory.h> 33 #define min0(x, y) (((x) <= (y)) ? (x) : (y)) 36 static void I_bessel(Float *x, Float *alpha,
int *nb,
37 int *ize, Float *bi,
int *ncalc);
41 Float bessel_i(Float x, Float alpha,
double expo)
45 #ifndef MATHLIB_STANDALONE 51 if (ISNAN(x) || ISNAN(alpha))
return x + alpha;
54 ML_ERROR(ME_RANGE,
"bessel_i");
62 return(bessel_i<Float>(x, -alpha, expo) +
64 bessel_k<Float>(x, -alpha, expo) *
65 ((ize == 1)? 2. : 2.*exp(-2.*x))/M_PI * sinpi(-alpha)));
67 nb = 1 + (int)trunc(na);
68 alpha -= (double)(nb-1);
69 #ifdef MATHLIB_STANDALONE 70 bi = (Float *) calloc(nb,
sizeof(Float));
71 if (!bi) MATHLIB_ERROR(
"%s", _(
"bessel_i allocation error"));
74 bi = (Float *) R_alloc((
size_t) nb,
sizeof(Float));
76 I_bessel(&x, &alpha, &nb, &ize, bi, &ncalc);
79 MATHLIB_WARNING4(_(
"bessel_i(%g): ncalc (=%d) != nb (=%d); alpha=%g. Arg. out of range?\n"),
82 MATHLIB_WARNING2(_(
"bessel_i(%g,nu=%g): precision lost in result\n"),
83 x, alpha+(
double)nb-1);
86 #ifdef MATHLIB_STANDALONE 97 Float bessel_i_ex(Float x, Float alpha,
double expo, Float *bi)
104 if (ISNAN(x) || ISNAN(alpha))
return x + alpha;
107 ML_ERROR(ME_RANGE,
"bessel_i");
115 return(bessel_i_ex(x, -alpha, expo, bi) +
117 bessel_k_ex(x, -alpha, expo, bi) *
118 ((ize == 1)? 2. : 2.*exp(-2.*x))/M_PI * sinpi(-alpha)));
120 nb = 1 + (int)trunc(na);
121 alpha -= (double)(nb-1);
122 I_bessel(&x, &alpha, &nb, &ize, bi, &ncalc);
125 MATHLIB_WARNING4(_(
"bessel_i(%g): ncalc (=%d) != nb (=%d); alpha=%g. Arg. out of range?\n"),
126 x, ncalc, nb, alpha);
128 MATHLIB_WARNING2(_(
"bessel_i(%g,nu=%g): precision lost in result\n"),
129 x, alpha+(
double)nb-1);
135 template<
class Float>
136 static void I_bessel(Float *x, Float *alpha,
int *nb,
137 int *ize, Float *bi,
int *ncalc)
233 const static double const__ = 1.585;
236 int nend, intx, nbmx, k, l, n, nstart;
237 Float pold, test, p, em, en, empal, emp2al, halfx,
238 aa, bb, cc, psave, plast, tover, psavel,
sum, nu, twonu;
248 if (*nb > 0 && *x >= 0. && (0. <= nu && nu < 1.) &&
249 (1 <= *ize && *ize <= 2) ) {
252 if(*ize == 1 && *x > exparg_BESS) {
253 for(k=1; k <= *nb; k++)
257 if(*ize == 2 && *x > xlrg_BESS_IJ) {
258 for(k=1; k <= *nb; k++)
262 intx = (int) trunc(*x);
263 if (*x >= rtnsig_BESS) {
269 en = (double) (n + n) + twonu;
275 test = ensig_BESS + ensig_BESS;
276 if (intx << 1 > nsig_BESS * 5) {
277 test = sqrt(test * p);
279 test /= R_pow_di(const__, intx);
286 tover = enten_BESS / ensig_BESS;
289 for (k = nstart; k <= nend; ++k) {
294 p = en * plast / *x + pold;
311 p = en * plast / *x + pold;
320 test = pold * plast / ensig_BESS;
321 test *= .5 - .5 / (bb * bb);
326 for (l = nstart; l <= nend; ++l) {
330 psave = en * psavel / *x + pold;
331 if (psave * psavel > test) {
342 en = (double)(n + n) + twonu;
346 test = fmax2(test,sqrt(plast * ensig_BESS) * sqrt(p + p));
356 p = en * plast / *x + pold;
367 em = (double) n - 1.;
369 emp2al = em - 1. + twonu;
370 sum = aa * empal * emp2al / em;
378 for (l = 1; l <= nend; ++l) {
388 for (l = 1; l <= nend; ++l) {
396 if(nend > 100 && aa > 1e200) {
398 cc = ldexp(cc, -900);
399 bb = ldexp(bb, -900);
400 sum = ldexp(sum,-900);
402 aa = en * bb / *x + cc;
412 sum = (sum + aa * empal) * emp2al / em;
420 sum = sum + sum + aa;
428 bi[n] = en * aa / *x + bb;
439 sum = (sum + bi[n] * empal) * emp2al / em;
447 for (l = 1; l <= nend; ++l) {
450 bi[n] = en * bi[n + 1] / *x + bi[n + 2];
457 sum = (sum + bi[n] * empal) * emp2al / em;
463 bi[1] = 2. * empal * bi[2] / *x + bi[3];
465 sum = sum + sum + bi[1];
473 sum *= (Rf_gamma_cody(1. + nu) * pow(*x * .5, -nu));
479 for (n = 1; n <= *nb; ++n) {
496 if (*x > enmten_BESS)
503 aa = pow(halfx, nu) / Rf_gamma_cody(empal);
507 bi[1] = aa + aa * bb / empal;
508 if (*x != 0. && bi[1] == 0.)
512 for (n = 2; n <= *nb; ++n)
519 tover = (enmten_BESS + enmten_BESS) / *x;
521 tover = enmten_BESS / bb;
522 for (n = 2; n <= *nb; ++n) {
526 if (aa <= tover * empal)
529 bi[n] = aa + aa * bb / empal;
530 if (bi[n] == 0. && *ncalc > n)
537 *ncalc = min0(*nb,0) - 1;
Type sum(Vector< Type > x)