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