1 /*
2 * Mathlib : A C Library of Special Functions
3 * Copyright (C) 1998-2012 Ross Ihaka and the R Core team.
4 *
5 * This program is free software; you can redistribute it and/or modify
6 * it under the terms of the GNU General Public License as published by
7 * the Free Software Foundation; either version 2 of the License, or
8 * (at your option) any later version.
9 *
10 * This program is distributed in the hope that it will be useful,
11 * but WITHOUT ANY WARRANTY; without even the implied warranty of
12 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
13 * GNU General Public License for more details.
14 *
15 * You should have received a copy of the GNU General Public License
16 * along with this program; if not, a copy is available at
17 * http://www.r-project.org/Licenses/
18 */
19
20 /* DESCRIPTION --> see below */
21
22
23 /* From http://www.netlib.org/specfun/ribesl Fortran translated by f2c,...
24 * ------------------------------=#---- Martin Maechler, ETH Zurich
25 */
26 #include "nmath.h"
27 #include "bessel.h"
28
29 #ifndef MATHLIB_STANDALONE
30 #include <R_ext/Memory.h>
31 #endif
32
33 #define min0(x, y) (((x) <= (y)) ? (x) : (y))
34
35 static void I_bessel(double *x, double *alpha, long *nb,
36 long *ize, double *bi, long *ncalc);
37
38 /* .Internal(besselI(*)) : */
bessel_i(double x,double alpha,double expo)39 double bessel_i(double x, double alpha, double expo)
40 {
41 long nb, ncalc, ize;
42 double na, *bi;
43 #ifndef MATHLIB_STANDALONE
44 const void *vmax;
45 #endif
46
47 #ifdef IEEE_754
48 /* NaNs propagated correctly */
49 if (ISNAN(x) || ISNAN(alpha)) return x + alpha;
50 #endif
51 if (x < 0) {
52 ML_ERROR(ME_RANGE, "bessel_i");
53 return ML_NAN;
54 }
55 ize = (long)expo;
56 na = floor(alpha);
57 if (alpha < 0) {
58 /* Using Abramowitz & Stegun 9.6.2 & 9.6.6
59 * this may not be quite optimal (CPU and accuracy wise) */
60 return(bessel_i(x, -alpha, expo) +
61 ((alpha == na) ? /* sin(pi * alpha) = 0 */ 0 :
62 bessel_k(x, -alpha, expo) *
63 ((ize == 1)? 2. : 2.*exp(-2.*x))/M_PI * sin(-M_PI * alpha)));
64 }
65 nb = 1 + (long)na;/* nb-1 <= alpha < nb */
66 alpha -= (double)(nb-1);
67 #ifdef MATHLIB_STANDALONE
68 bi = (double *) calloc(nb, sizeof(double));
69 if (!bi) MATHLIB_ERROR("%s", _("bessel_i allocation error"));
70 #else
71 vmax = vmaxget();
72 bi = (double *) R_alloc((size_t) nb, sizeof(double));
73 #endif
74 I_bessel(&x, &alpha, &nb, &ize, bi, &ncalc);
75 if(ncalc != nb) {/* error input */
76 if(ncalc < 0)
77 MATHLIB_WARNING4(_("bessel_i(%g): ncalc (=%ld) != nb (=%ld); alpha=%g. Arg. out of range?\n"),
78 x, ncalc, nb, alpha);
79 else
80 MATHLIB_WARNING2(_("bessel_i(%g,nu=%g): precision lost in result\n"),
81 x, alpha+(double)nb-1);
82 }
83 x = bi[nb-1];
84 #ifdef MATHLIB_STANDALONE
85 free(bi);
86 #else
87 vmaxset(vmax);
88 #endif
89 return x;
90 }
91
92 /* modified version of bessel_i that accepts a work array instead of
93 allocating one. */
bessel_i_ex(double x,double alpha,double expo,double * bi)94 double bessel_i_ex(double x, double alpha, double expo, double *bi)
95 {
96 long nb, ncalc, ize;
97 double na;
98
99 #ifdef IEEE_754
100 /* NaNs propagated correctly */
101 if (ISNAN(x) || ISNAN(alpha)) return x + alpha;
102 #endif
103 if (x < 0) {
104 ML_ERROR(ME_RANGE, "bessel_i");
105 return ML_NAN;
106 }
107 ize = (long)expo;
108 na = floor(alpha);
109 if (alpha < 0) {
110 /* Using Abramowitz & Stegun 9.6.2 & 9.6.6
111 * this may not be quite optimal (CPU and accuracy wise) */
112 return(bessel_i_ex(x, -alpha, expo, bi) +
113 ((alpha == na) ? 0 :
114 bessel_k_ex(x, -alpha, expo, bi) *
115 ((ize == 1)? 2. : 2.*exp(-2.*x))/M_PI * sin(-M_PI * alpha)));
116 }
117 nb = 1 + (long)na;/* nb-1 <= alpha < nb */
118 alpha -= (double)(nb-1);
119 I_bessel(&x, &alpha, &nb, &ize, bi, &ncalc);
120 if(ncalc != nb) {/* error input */
121 if(ncalc < 0)
122 MATHLIB_WARNING4(_("bessel_i(%g): ncalc (=%ld) != nb (=%ld); alpha=%g. Arg. out of range?\n"),
123 x, ncalc, nb, alpha);
124 else
125 MATHLIB_WARNING2(_("bessel_i(%g,nu=%g): precision lost in result\n"),
126 x, alpha+(double)nb-1);
127 }
128 x = bi[nb-1];
129 return x;
130 }
131
I_bessel(double * x,double * alpha,long * nb,long * ize,double * bi,long * ncalc)132 static void I_bessel(double *x, double *alpha, long *nb,
133 long *ize, double *bi, long *ncalc)
134 {
135 /* -------------------------------------------------------------------
136
137 This routine calculates Bessel functions I_(N+ALPHA) (X)
138 for non-negative argument X, and non-negative order N+ALPHA,
139 with or without exponential scaling.
140
141
142 Explanation of variables in the calling sequence
143
144 X - Non-negative argument for which
145 I's or exponentially scaled I's (I*EXP(-X))
146 are to be calculated. If I's are to be calculated,
147 X must be less than exparg_BESS (IZE=1) or xlrg_BESS_IJ (IZE=2),
148 (see bessel.h).
149 ALPHA - Fractional part of order for which
150 I's or exponentially scaled I's (I*EXP(-X)) are
151 to be calculated. 0 <= ALPHA < 1.0.
152 NB - Number of functions to be calculated, NB > 0.
153 The first function calculated is of order ALPHA, and the
154 last is of order (NB - 1 + ALPHA).
155 IZE - Type. IZE = 1 if unscaled I's are to be calculated,
156 = 2 if exponentially scaled I's are to be calculated.
157 BI - Output vector of length NB. If the routine
158 terminates normally (NCALC=NB), the vector BI contains the
159 functions I(ALPHA,X) through I(NB-1+ALPHA,X), or the
160 corresponding exponentially scaled functions.
161 NCALC - Output variable indicating possible errors.
162 Before using the vector BI, the user should check that
163 NCALC=NB, i.e., all orders have been calculated to
164 the desired accuracy. See error returns below.
165
166
167 *******************************************************************
168 *******************************************************************
169
170 Error returns
171
172 In case of an error, NCALC != NB, and not all I's are
173 calculated to the desired accuracy.
174
175 NCALC < 0: An argument is out of range. For example,
176 NB <= 0, IZE is not 1 or 2, or IZE=1 and ABS(X) >= EXPARG_BESS.
177 In this case, the BI-vector is not calculated, and NCALC is
178 set to MIN0(NB,0)-1 so that NCALC != NB.
179
180 NB > NCALC > 0: Not all requested function values could
181 be calculated accurately. This usually occurs because NB is
182 much larger than ABS(X). In this case, BI[N] is calculated
183 to the desired accuracy for N <= NCALC, but precision
184 is lost for NCALC < N <= NB. If BI[N] does not vanish
185 for N > NCALC (because it is too small to be represented),
186 and BI[N]/BI[NCALC] = 10**(-K), then only the first NSIG-K
187 significant figures of BI[N] can be trusted.
188
189
190 Intrinsic functions required are:
191
192 DBLE, EXP, gamma_cody, INT, MAX, MIN, REAL, SQRT
193
194
195 Acknowledgement
196
197 This program is based on a program written by David J.
198 Sookne (2) that computes values of the Bessel functions J or
199 I of float argument and long order. Modifications include
200 the restriction of the computation to the I Bessel function
201 of non-negative float argument, the extension of the computation
202 to arbitrary positive order, the inclusion of optional
203 exponential scaling, and the elimination of most underflow.
204 An earlier version was published in (3).
205
206 References: "A Note on Backward Recurrence Algorithms," Olver,
207 F. W. J., and Sookne, D. J., Math. Comp. 26, 1972,
208 pp 941-947.
209
210 "Bessel Functions of Real Argument and Integer Order,"
211 Sookne, D. J., NBS Jour. of Res. B. 77B, 1973, pp
212 125-132.
213
214 "ALGORITHM 597, Sequence of Modified Bessel Functions
215 of the First Kind," Cody, W. J., Trans. Math. Soft.,
216 1983, pp. 242-245.
217
218 Latest modification: May 30, 1989
219
220 Modified by: W. J. Cody and L. Stoltz
221 Applied Mathematics Division
222 Argonne National Laboratory
223 Argonne, IL 60439
224 */
225
226 /*-------------------------------------------------------------------
227 Mathematical constants
228 -------------------------------------------------------------------*/
229 const static double const__ = 1.585;
230
231 /* Local variables */
232 long nend, intx, nbmx, k, l, n, nstart;
233 double pold, test, p, em, en, empal, emp2al, halfx,
234 aa, bb, cc, psave, plast, tover, psavel, sum, nu, twonu;
235
236 /*Parameter adjustments */
237 --bi;
238 nu = *alpha;
239 twonu = nu + nu;
240
241 /*-------------------------------------------------------------------
242 Check for X, NB, OR IZE out of range.
243 ------------------------------------------------------------------- */
244 if (*nb > 0 && *x >= 0. && (0. <= nu && nu < 1.) &&
245 (1 <= *ize && *ize <= 2) ) {
246
247 *ncalc = *nb;
248 if(*ize == 1 && *x > exparg_BESS) {
249 for(k=1; k <= *nb; k++)
250 bi[k]=ML_POSINF; /* the limit *is* = Inf */
251 return;
252 }
253 if(*ize == 2 && *x > xlrg_BESS_IJ) {
254 for(k=1; k <= *nb; k++)
255 bi[k]= 0.; /* The limit exp(-x) * I_nu(x) --> 0 : */
256 return;
257 }
258 intx = (long) (*x);/* fine, since *x <= xlrg_BESS_IJ <<< LONG_MAX */
259 if (*x >= rtnsig_BESS) { /* "non-small" x ( >= 1e-4 ) */
260 /* -------------------------------------------------------------------
261 Initialize the forward sweep, the P-sequence of Olver
262 ------------------------------------------------------------------- */
263 nbmx = *nb - intx;
264 n = intx + 1;
265 en = (double) (n + n) + twonu;
266 plast = 1.;
267 p = en / *x;
268 /* ------------------------------------------------
269 Calculate general significance test
270 ------------------------------------------------ */
271 test = ensig_BESS + ensig_BESS;
272 if (intx << 1 > nsig_BESS * 5) {
273 test = sqrt(test * p);
274 } else {
275 test /= pow(const__, (double)intx);
276 }
277 if (nbmx >= 3) {
278 /* --------------------------------------------------
279 Calculate P-sequence until N = NB-1
280 Check for possible overflow.
281 ------------------------------------------------ */
282 tover = enten_BESS / ensig_BESS;
283 nstart = intx + 2;
284 nend = *nb - 1;
285 for (k = nstart; k <= nend; ++k) {
286 n = k;
287 en += 2.;
288 pold = plast;
289 plast = p;
290 p = en * plast / *x + pold;
291 if (p > tover) {
292 /* ------------------------------------------------
293 To avoid overflow, divide P-sequence by TOVER.
294 Calculate P-sequence until ABS(P) > 1.
295 ---------------------------------------------- */
296 tover = enten_BESS;
297 p /= tover;
298 plast /= tover;
299 psave = p;
300 psavel = plast;
301 nstart = n + 1;
302 do {
303 ++n;
304 en += 2.;
305 pold = plast;
306 plast = p;
307 p = en * plast / *x + pold;
308 }
309 while (p <= 1.);
310
311 bb = en / *x;
312 /* ------------------------------------------------
313 Calculate backward test, and find NCALC,
314 the highest N such that the test is passed.
315 ------------------------------------------------ */
316 test = pold * plast / ensig_BESS;
317 test *= .5 - .5 / (bb * bb);
318 p = plast * tover;
319 --n;
320 en -= 2.;
321 nend = min0(*nb,n);
322 for (l = nstart; l <= nend; ++l) {
323 *ncalc = l;
324 pold = psavel;
325 psavel = psave;
326 psave = en * psavel / *x + pold;
327 if (psave * psavel > test) {
328 goto L90;
329 }
330 }
331 *ncalc = nend + 1;
332 L90:
333 --(*ncalc);
334 goto L120;
335 }
336 }
337 n = nend;
338 en = (double)(n + n) + twonu;
339 /*---------------------------------------------------
340 Calculate special significance test for NBMX > 2.
341 --------------------------------------------------- */
342 test = fmax2(test,sqrt(plast * ensig_BESS) * sqrt(p + p));
343 }
344 /* --------------------------------------------------------
345 Calculate P-sequence until significance test passed.
346 -------------------------------------------------------- */
347 do {
348 ++n;
349 en += 2.;
350 pold = plast;
351 plast = p;
352 p = en * plast / *x + pold;
353 } while (p < test);
354
355 L120:
356 /* -------------------------------------------------------------------
357 Initialize the backward recursion and the normalization sum.
358 ------------------------------------------------------------------- */
359 ++n;
360 en += 2.;
361 bb = 0.;
362 aa = 1. / p;
363 em = (double) n - 1.;
364 empal = em + nu;
365 emp2al = em - 1. + twonu;
366 sum = aa * empal * emp2al / em;
367 nend = n - *nb;
368 if (nend < 0) {
369 /* -----------------------------------------------------
370 N < NB, so store BI[N] and set higher orders to 0..
371 ----------------------------------------------------- */
372 bi[n] = aa;
373 nend = -nend;
374 for (l = 1; l <= nend; ++l) {
375 bi[n + l] = 0.;
376 }
377 } else {
378 if (nend > 0) {
379 /* -----------------------------------------------------
380 Recur backward via difference equation,
381 calculating (but not storing) BI[N], until N = NB.
382 --------------------------------------------------- */
383
384 for (l = 1; l <= nend; ++l) {
385 --n;
386 en -= 2.;
387 cc = bb;
388 bb = aa;
389 /* for x ~= 1500, sum would overflow to 'inf' here,
390 * and the final bi[] /= sum would give 0 wrongly;
391 * RE-normalize (aa, sum) here -- no need to undo */
392 if(nend > 100 && aa > 1e200) {
393 /* multiply by 2^-900 = 1.18e-271 */
394 cc = ldexp(cc, -900);
395 bb = ldexp(bb, -900);
396 sum = ldexp(sum,-900);
397 }
398 aa = en * bb / *x + cc;
399 em -= 1.;
400 emp2al -= 1.;
401 if (n == 1) {
402 break;
403 }
404 if (n == 2) {
405 emp2al = 1.;
406 }
407 empal -= 1.;
408 sum = (sum + aa * empal) * emp2al / em;
409 }
410 }
411 /* ---------------------------------------------------
412 Store BI[NB]
413 --------------------------------------------------- */
414 bi[n] = aa;
415 if (*nb <= 1) {
416 sum = sum + sum + aa;
417 goto L230;
418 }
419 /* -------------------------------------------------
420 Calculate and Store BI[NB-1]
421 ------------------------------------------------- */
422 --n;
423 en -= 2.;
424 bi[n] = en * aa / *x + bb;
425 if (n == 1) {
426 goto L220;
427 }
428 em -= 1.;
429 if (n == 2)
430 emp2al = 1.;
431 else
432 emp2al -= 1.;
433
434 empal -= 1.;
435 sum = (sum + bi[n] * empal) * emp2al / em;
436 }
437 nend = n - 2;
438 if (nend > 0) {
439 /* --------------------------------------------
440 Calculate via difference equation
441 and store BI[N], until N = 2.
442 ------------------------------------------ */
443 for (l = 1; l <= nend; ++l) {
444 --n;
445 en -= 2.;
446 bi[n] = en * bi[n + 1] / *x + bi[n + 2];
447 em -= 1.;
448 if (n == 2)
449 emp2al = 1.;
450 else
451 emp2al -= 1.;
452 empal -= 1.;
453 sum = (sum + bi[n] * empal) * emp2al / em;
454 }
455 }
456 /* ----------------------------------------------
457 Calculate BI[1]
458 -------------------------------------------- */
459 bi[1] = 2. * empal * bi[2] / *x + bi[3];
460 L220:
461 sum = sum + sum + bi[1];
462
463 L230:
464 /* ---------------------------------------------------------
465 Normalize. Divide all BI[N] by sum.
466 --------------------------------------------------------- */
467 if (nu != 0.)
468 sum *= (gamma_cody(1. + nu) * pow(*x * .5, -nu));
469 if (*ize == 1)
470 sum *= exp(-(*x));
471 aa = enmten_BESS;
472 if (sum > 1.)
473 aa *= sum;
474 for (n = 1; n <= *nb; ++n) {
475 if (bi[n] < aa)
476 bi[n] = 0.;
477 else
478 bi[n] /= sum;
479 }
480 return;
481 } else { /* small x < 1e-4 */
482 /* -----------------------------------------------------------
483 Two-term ascending series for small X.
484 -----------------------------------------------------------*/
485 aa = 1.;
486 empal = 1. + nu;
487 #ifdef IEEE_754
488 /* No need to check for underflow */
489 halfx = .5 * *x;
490 #else
491 if (*x > enmten_BESS) */
492 halfx = .5 * *x;
493 else
494 halfx = 0.;
495 #endif
496 if (nu != 0.)
497 aa = pow(halfx, nu) / gamma_cody(empal);
498 if (*ize == 2)
499 aa *= exp(-(*x));
500 bb = halfx * halfx;
501 bi[1] = aa + aa * bb / empal;
502 if (*x != 0. && bi[1] == 0.)
503 *ncalc = 0;
504 if (*nb > 1) {
505 if (*x == 0.) {
506 for (n = 2; n <= *nb; ++n)
507 bi[n] = 0.;
508 } else {
509 /* -------------------------------------------------
510 Calculate higher-order functions.
511 ------------------------------------------------- */
512 cc = halfx;
513 tover = (enmten_BESS + enmten_BESS) / *x;
514 if (bb != 0.)
515 tover = enmten_BESS / bb;
516 for (n = 2; n <= *nb; ++n) {
517 aa /= empal;
518 empal += 1.;
519 aa *= cc;
520 if (aa <= tover * empal)
521 bi[n] = aa = 0.;
522 else
523 bi[n] = aa + aa * bb / empal;
524 if (bi[n] == 0. && *ncalc > n)
525 *ncalc = n - 1;
526 }
527 }
528 }
529 }
530 } else { /* argument out of range */
531 *ncalc = min0(*nb,0) - 1;
532 }
533 }
534