TMB Documentation  v1.9.11
bessel_j.cpp
1 /*
2  * Mathlib : A C Library of Special Functions
3  * Copyright (C) 1998-2015 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/rjbesl Fortran translated by f2c,...
24  * ------------------------------=#---- Martin Maechler, ETH Zurich
25  * Additional code for nu == alpha < 0 MM
26  */
27 
28 #include "bessel.h"
29 
30 #ifndef MATHLIB_STANDALONE
31 #include <R_ext/Memory.h>
32 #endif
33 
34 #define min0(x, y) (((x) <= (y)) ? (x) : (y))
35 
36 template<class Float>
37 static void J_bessel(Float *x, Float *alpha, int *nb,
38  Float *b, int *ncalc);
39 
40 // unused now from R
41 template<class Float>
42 Float bessel_j(Float x, Float alpha)
43 {
44  int nb, ncalc;
45  Float na, *bj;
46 #ifndef MATHLIB_STANDALONE
47  const void *vmax;
48 #endif
49 
50 #ifdef IEEE_754
51  /* NaNs propagated correctly */
52  if (ISNAN(x) || ISNAN(alpha)) return x + alpha;
53 #endif
54  if (x < 0) {
55  ML_ERROR(ME_RANGE, "bessel_j");
56  return ML_NAN;
57  }
58  na = floor(alpha);
59  if (alpha < 0) {
60  /* Using Abramowitz & Stegun 9.1.2
61  * this may not be quite optimal (CPU and accuracy wise) */
62  return(((alpha - na == 0.5) ? 0 : bessel_j<Float>(x, -alpha) * cospi(alpha)) +
63  ((alpha == na ) ? 0 : bessel_y<Float>(x, -alpha) * sinpi(alpha)));
64  }
65  else if (alpha > 1e7) {
66  MATHLIB_WARNING(_("besselJ(x, nu): nu=%g too large for bessel_j() algorithm"),
67  alpha);
68  return ML_NAN;
69  }
70  nb = 1 + (int)trunc(na); /* nb-1 <= alpha < nb */
71  alpha -= (double)(nb-1); // ==> alpha' in [0, 1)
72 #ifdef MATHLIB_STANDALONE
73  bj = (Float *) calloc(nb, sizeof(Float));
74  if (!bj) MATHLIB_ERROR("%s", _("bessel_j allocation error"));
75 #else
76  vmax = vmaxget();
77  bj = (Float *) R_alloc((size_t) nb, sizeof(Float));
78 #endif
79  J_bessel(&x, &alpha, &nb, bj, &ncalc);
80  if(ncalc != nb) {/* error input */
81  if(ncalc < 0)
82  MATHLIB_WARNING4(_("bessel_j(%g): ncalc (=%d) != nb (=%d); alpha=%g. Arg. out of range?\n"),
83  x, ncalc, nb, alpha);
84  else
85  MATHLIB_WARNING2(_("bessel_j(%g,nu=%g): precision lost in result\n"),
86  x, alpha+(double)nb-1);
87  }
88  x = bj[nb-1];
89 #ifdef MATHLIB_STANDALONE
90  free(bj);
91 #else
92  vmaxset(vmax);
93 #endif
94  return x;
95 }
96 
97 /* Called from R: modified version of bessel_j(), accepting a work array
98  * instead of allocating one. */
99 template<class Float>
100 Float bessel_j_ex(Float x, Float alpha, Float *bj)
101 {
102  int nb, ncalc;
103  double na;
104 
105 #ifdef IEEE_754
106  /* NaNs propagated correctly */
107  if (ISNAN(x) || ISNAN(alpha)) return x + alpha;
108 #endif
109  if (x < 0) {
110  ML_ERROR(ME_RANGE, "bessel_j");
111  return ML_NAN;
112  }
113  na = floor(alpha);
114  if (alpha < 0) {
115  /* Using Abramowitz & Stegun 9.1.2
116  * this may not be quite optimal (CPU and accuracy wise) */
117  return(((alpha - na == 0.5) ? 0 : bessel_j_ex(x, -alpha, bj) * cospi(alpha)) +
118  ((alpha == na ) ? 0 : bessel_y_ex(x, -alpha, bj) * sinpi(alpha)));
119  }
120  else if (alpha > 1e7) {
121  MATHLIB_WARNING(_("besselJ(x, nu): nu=%g too large for bessel_j() algorithm"),
122  alpha);
123  return ML_NAN;
124  }
125  nb = 1 + (int)trunc(na); /* nb-1 <= alpha < nb */
126  alpha -= (double)(nb-1); // ==> alpha' in [0, 1)
127  J_bessel(&x, &alpha, &nb, bj, &ncalc);
128  if(ncalc != nb) {/* error input */
129  if(ncalc < 0)
130  MATHLIB_WARNING4(_("bessel_j(%g): ncalc (=%d) != nb (=%d); alpha=%g. Arg. out of range?\n"),
131  x, ncalc, nb, alpha);
132  else
133  MATHLIB_WARNING2(_("bessel_j(%g,nu=%g): precision lost in result\n"),
134  x, alpha+(double)nb-1);
135  }
136  x = bj[nb-1];
137  return x;
138 }
139 
140 template<class Float>
141 static void J_bessel(Float *x, Float *alpha, int *nb,
142  Float *b, int *ncalc)
143 {
144 /*
145  Calculates Bessel functions J_{n+alpha} (x)
146  for non-negative argument x, and non-negative order n+alpha, n = 0,1,..,nb-1.
147 
148  Explanation of variables in the calling sequence.
149 
150  X - Non-negative argument for which J's are to be calculated.
151  ALPHA - Fractional part of order for which
152  J's are to be calculated. 0 <= ALPHA < 1.
153  NB - Number of functions to be calculated, NB >= 1.
154  The first function calculated is of order ALPHA, and the
155  last is of order (NB - 1 + ALPHA).
156  B - Output vector of length NB. If RJBESL
157  terminates normally (NCALC=NB), the vector B contains the
158  functions J/ALPHA/(X) through J/NB-1+ALPHA/(X).
159  NCALC - Output variable indicating possible errors.
160  Before using the vector B, the user should check that
161  NCALC=NB, i.e., all orders have been calculated to
162  the desired accuracy. See the following
163 
164  ****************************************************************
165 
166  Error return codes
167 
168  In case of an error, NCALC != NB, and not all J's are
169  calculated to the desired accuracy.
170 
171  NCALC < 0: An argument is out of range. For example,
172  NBES <= 0, ALPHA < 0 or > 1, or X is too large.
173  In this case, b[1] is set to zero, the remainder of the
174  B-vector is not calculated, and NCALC is set to
175  MIN(NB,0)-1 so that NCALC != NB.
176 
177  NB > NCALC > 0: Not all requested function values could
178  be calculated accurately. This usually occurs because NB is
179  much larger than ABS(X). In this case, b[N] is calculated
180  to the desired accuracy for N <= NCALC, but precision
181  is lost for NCALC < N <= NB. If b[N] does not vanish
182  for N > NCALC (because it is too small to be represented),
183  and b[N]/b[NCALC] = 10^(-K), then only the first NSIG - K
184  significant figures of b[N] can be trusted.
185 
186 
187  Acknowledgement
188 
189  This program is based on a program written by David J. Sookne
190  (2) that computes values of the Bessel functions J or I of float
191  argument and long order. Modifications include the restriction
192  of the computation to the J Bessel function of non-negative float
193  argument, the extension of the computation to arbitrary positive
194  order, and the elimination of most underflow.
195 
196  References:
197 
198  Olver, F.W.J., and Sookne, D.J. (1972)
199  "A Note on Backward Recurrence Algorithms";
200  Math. Comp. 26, 941-947.
201 
202  Sookne, D.J. (1973)
203  "Bessel Functions of Real Argument and Integer Order";
204  NBS Jour. of Res. B. 77B, 125-132.
205 
206  Latest modification: March 19, 1990
207 
208  Author: W. J. Cody
209  Applied Mathematics Division
210  Argonne National Laboratory
211  Argonne, IL 60439
212  *******************************************************************
213  */
214 
215 /* ---------------------------------------------------------------------
216  Mathematical constants
217 
218  PI2 = 2 / PI
219  TWOPI1 = first few significant digits of 2 * PI
220  TWOPI2 = (2*PI - TWOPI1) to working precision, i.e.,
221  TWOPI1 + TWOPI2 = 2 * PI to extra precision.
222  --------------------------------------------------------------------- */
223  const static double pi2 = .636619772367581343075535;
224  const static double twopi1 = 6.28125;
225  const static double twopi2 = .001935307179586476925286767;
226 
227 /*---------------------------------------------------------------------
228  * Factorial(N)
229  *--------------------------------------------------------------------- */
230  const static double fact[25] = { 1.,1.,2.,6.,24.,120.,720.,5040.,40320.,
231  362880.,3628800.,39916800.,479001600.,6227020800.,87178291200.,
232  1.307674368e12,2.0922789888e13,3.55687428096e14,6.402373705728e15,
233  1.21645100408832e17,2.43290200817664e18,5.109094217170944e19,
234  1.12400072777760768e21,2.585201673888497664e22,
235  6.2044840173323943936e23 };
236 
237  /* Local variables */
238  int nend, intx, nbmx, i, j, k, l, m, n, nstart;
239 
240  Float nu, twonu, capp, capq, pold, vcos, test, vsin;
241  Float p, s, t, z, alpem, halfx, aa, bb, cc, psave, plast;
242  Float tover, t1, alp2em, em, en, xc, xk, xm, psavel, gnu, xin, sum;
243 
244 
245  /* Parameter adjustment */
246  --b;
247 
248  nu = *alpha;
249  twonu = nu + nu;
250 
251  /*-------------------------------------------------------------------
252  Check for out of range arguments.
253  -------------------------------------------------------------------*/
254  if (*nb > 0 && *x >= 0. && 0. <= nu && nu < 1.) {
255 
256  *ncalc = *nb;
257  if(*x > xlrg_BESS_IJ) {
258  ML_ERROR(ME_RANGE, "J_bessel");
259  /* indeed, the limit is 0,
260  * but the cutoff happens too early */
261  for(i=1; i <= *nb; i++)
262  b[i] = 0.; /*was ML_POSINF (really nonsense) */
263  return;
264  }
265  intx = (int) trunc(*x);
266  /* Initialize result array to zero. */
267  for (i = 1; i <= *nb; ++i)
268  b[i] = 0.;
269 
270  /*===================================================================
271  Branch into 3 cases :
272  1) use 2-term ascending series for small X
273  2) use asymptotic form for large X when NB is not too large
274  3) use recursion otherwise
275  ===================================================================*/
276 
277  if (*x < rtnsig_BESS) {
278  /* ---------------------------------------------------------------
279  Two-term ascending series for small X.
280  --------------------------------------------------------------- */
281  alpem = 1. + nu;
282 
283  halfx = (*x > enmten_BESS) ? .5 * *x : 0.;
284  aa = (nu != 0.) ? pow(halfx, nu) / (nu * Rf_gamma_cody(nu)) : 1.;
285  bb = (*x + 1. > 1.)? -halfx * halfx : 0.;
286  b[1] = aa + aa * bb / alpem;
287  if (*x != 0. && b[1] == 0.)
288  *ncalc = 0;
289 
290  if (*nb != 1) {
291  if (*x <= 0.) {
292  for (n = 2; n <= *nb; ++n)
293  b[n] = 0.;
294  }
295  else {
296  /* ----------------------------------------------
297  Calculate higher order functions.
298  ---------------------------------------------- */
299  if (bb == 0.)
300  tover = (enmten_BESS + enmten_BESS) / *x;
301  else
302  tover = enmten_BESS / bb;
303  cc = halfx;
304  for (n = 2; n <= *nb; ++n) {
305  aa /= alpem;
306  alpem += 1.;
307  aa *= cc;
308  if (aa <= tover * alpem)
309  aa = 0.;
310 
311  b[n] = aa + aa * bb / alpem;
312  if (b[n] == 0. && *ncalc > n)
313  *ncalc = n - 1;
314  }
315  }
316  }
317  } else if (*x > 25. && *nb <= intx + 1) {
318  /* ------------------------------------------------------------
319  Asymptotic series for X > 25 (and not too large nb)
320  ------------------------------------------------------------ */
321  xc = sqrt(pi2 / *x);
322  xin = 1 / (64 * *x * *x);
323  if (*x >= 130.) m = 4;
324  else if (*x >= 35.) m = 8;
325  else m = 11;
326  xm = 4. * (double) m;
327  /* ------------------------------------------------
328  Argument reduction for SIN and COS routines.
329  ------------------------------------------------ */
330  t = trunc(*x / (twopi1 + twopi2) + .5);
331  z = (*x - t * twopi1) - t * twopi2 - (nu + .5) / pi2;
332  vsin = sin(z);
333  vcos = cos(z);
334  gnu = twonu;
335  for (i = 1; i <= 2; ++i) {
336  s = (xm - 1. - gnu) * (xm - 1. + gnu) * xin * .5;
337  t = (gnu - (xm - 3.)) * (gnu + (xm - 3.));
338  t1= (gnu - (xm + 1.)) * (gnu + (xm + 1.));
339  k = m + m;
340  capp = s * t / fact[k];
341  capq = s * t1/ fact[k + 1];
342  xk = xm;
343  for (; k >= 4; k -= 2) {/* k + 2(j-2) == 2m */
344  xk -= 4.;
345  s = (xk - 1. - gnu) * (xk - 1. + gnu);
346  t1 = t;
347  t = (gnu - (xk - 3.)) * (gnu + (xk - 3.));
348  capp = (capp + 1. / fact[k - 2]) * s * t * xin;
349  capq = (capq + 1. / fact[k - 1]) * s * t1 * xin;
350 
351  }
352  capp += 1.;
353  capq = (capq + 1.) * (gnu * gnu - 1.) * (.125 / *x);
354  b[i] = xc * (capp * vcos - capq * vsin);
355  if (*nb == 1)
356  return;
357 
358  /* vsin <--> vcos */ t = vsin; vsin = -vcos; vcos = t;
359  gnu += 2.;
360  }
361  /* -----------------------------------------------
362  If NB > 2, compute J(X,ORDER+I) for I = 2, NB-1
363  ----------------------------------------------- */
364  if (*nb > 2)
365  for (gnu = twonu + 2., j = 3; j <= *nb; j++, gnu += 2.)
366  b[j] = gnu * b[j - 1] / *x - b[j - 2];
367  }
368  else {
369  /* rtnsig_BESS <= x && ( x <= 25 || intx+1 < *nb ) :
370  --------------------------------------------------------
371  Use recurrence to generate results.
372  First initialize the calculation of P*S.
373  -------------------------------------------------------- */
374  nbmx = *nb - intx;
375  n = intx + 1;
376  en = (double)(n + n) + twonu;
377  plast = 1.;
378  p = en / *x;
379  /* ---------------------------------------------------
380  Calculate general significance test.
381  --------------------------------------------------- */
382  test = ensig_BESS + ensig_BESS;
383  if (nbmx >= 3) {
384  /* ------------------------------------------------------------
385  Calculate P*S until N = NB-1. Check for possible overflow.
386  ---------------------------------------------------------- */
387  tover = enten_BESS / ensig_BESS;
388  nstart = intx + 2;
389  nend = *nb - 1;
390  en = (double) (nstart + nstart) - 2. + twonu;
391  for (k = nstart; k <= nend; ++k) {
392  n = k;
393  en += 2.;
394  pold = plast;
395  plast = p;
396  p = en * plast / *x - pold;
397  if (p > tover) {
398  /* -------------------------------------------
399  To avoid overflow, divide P*S by TOVER.
400  Calculate P*S until ABS(P) > 1.
401  -------------------------------------------*/
402  tover = enten_BESS;
403  p /= tover;
404  plast /= tover;
405  psave = p;
406  psavel = plast;
407  nstart = n + 1;
408  do {
409  ++n;
410  en += 2.;
411  pold = plast;
412  plast = p;
413  p = en * plast / *x - pold;
414  } while (p <= 1.);
415 
416  bb = en / *x;
417  /* -----------------------------------------------
418  Calculate backward test and find NCALC,
419  the highest N such that the test is passed.
420  ----------------------------------------------- */
421  test = pold * plast * (.5 - .5 / (bb * bb));
422  test /= ensig_BESS;
423  p = plast * tover;
424  --n;
425  en -= 2.;
426  nend = min0(*nb,n);
427  for (l = nstart; l <= nend; ++l) {
428  pold = psavel;
429  psavel = psave;
430  psave = en * psavel / *x - pold;
431  if (psave * psavel > test) {
432  *ncalc = l - 1;
433  goto L190;
434  }
435  }
436  *ncalc = nend;
437  goto L190;
438  }
439  }
440  n = nend;
441  en = (double) (n + n) + twonu;
442  /* -----------------------------------------------------
443  Calculate special significance test for NBMX > 2.
444  -----------------------------------------------------*/
445  test = fmax2(test, sqrt(plast * ensig_BESS) * sqrt(p + p));
446  }
447  /* ------------------------------------------------
448  Calculate P*S until significance test passes. */
449  do {
450  ++n;
451  en += 2.;
452  pold = plast;
453  plast = p;
454  p = en * plast / *x - pold;
455  } while (p < test);
456 
457 L190:
458  /*---------------------------------------------------------------
459  Initialize the backward recursion and the normalization sum.
460  --------------------------------------------------------------- */
461  ++n;
462  en += 2.;
463  bb = 0.;
464  aa = 1. / p;
465  m = n / 2;
466  em = (double)m;
467  m = (n << 1) - (m << 2);/* = 2 n - 4 (n/2)
468  = 0 for even, 2 for odd n */
469  if (m == 0)
470  sum = 0.;
471  else {
472  alpem = em - 1. + nu;
473  alp2em = em + em + nu;
474  sum = aa * alpem * alp2em / em;
475  }
476  nend = n - *nb;
477  /* if (nend > 0) */
478  /* --------------------------------------------------------
479  Recur backward via difference equation, calculating
480  (but not storing) b[N], until N = NB.
481  -------------------------------------------------------- */
482  for (l = 1; l <= nend; ++l) {
483  --n;
484  en -= 2.;
485  cc = bb;
486  bb = aa;
487  aa = en * bb / *x - cc;
488  m = m ? 0 : 2; /* m = 2 - m failed on gcc4-20041019 */
489  if (m != 0) {
490  em -= 1.;
491  alp2em = em + em + nu;
492  if (n == 1)
493  break;
494 
495  alpem = em - 1. + nu;
496  if (alpem == 0.)
497  alpem = 1.;
498  sum = (sum + aa * alp2em) * alpem / em;
499  }
500  }
501  /*--------------------------------------------------
502  Store b[NB].
503  --------------------------------------------------*/
504  b[n] = aa;
505  if (nend >= 0) {
506  if (*nb <= 1) {
507  if (nu + 1. == 1.)
508  alp2em = 1.;
509  else
510  alp2em = nu;
511  sum += b[1] * alp2em;
512  goto L250;
513  }
514  else {/*-- nb >= 2 : ---------------------------
515  Calculate and store b[NB-1].
516  ----------------------------------------*/
517  --n;
518  en -= 2.;
519  b[n] = en * aa / *x - bb;
520  if (n == 1)
521  goto L240;
522 
523  m = m ? 0 : 2; /* m = 2 - m failed on gcc4-20041019 */
524  if (m != 0) {
525  em -= 1.;
526  alp2em = em + em + nu;
527  alpem = em - 1. + nu;
528  if (alpem == 0.)
529  alpem = 1.;
530  sum = (sum + b[n] * alp2em) * alpem / em;
531  }
532  }
533  }
534 
535  /* if (n - 2 != 0) */
536  /* --------------------------------------------------------
537  Calculate via difference equation and store b[N],
538  until N = 2.
539  -------------------------------------------------------- */
540  for (n = n-1; n >= 2; n--) {
541  en -= 2.;
542  b[n] = en * b[n + 1] / *x - b[n + 2];
543  m = m ? 0 : 2; /* m = 2 - m failed on gcc4-20041019 */
544  if (m != 0) {
545  em -= 1.;
546  alp2em = em + em + nu;
547  alpem = em - 1. + nu;
548  if (alpem == 0.)
549  alpem = 1.;
550  sum = (sum + b[n] * alp2em) * alpem / em;
551  }
552  }
553  /* ---------------------------------------
554  Calculate b[1].
555  -----------------------------------------*/
556  b[1] = 2. * (nu + 1.) * b[2] / *x - b[3];
557 
558 L240:
559  em -= 1.;
560  alp2em = em + em + nu;
561  if (alp2em == 0.)
562  alp2em = 1.;
563  sum += b[1] * alp2em;
564 
565 L250:
566  /* ---------------------------------------------------
567  Normalize. Divide all b[N] by sum.
568  ---------------------------------------------------*/
569 /* if (nu + 1. != 1.) poor test */
570  if(fabs(nu) > 1e-15)
571  sum *= (Rf_gamma_cody(nu) * pow(.5* *x, -nu));
572 
573  aa = enmten_BESS;
574  if (sum > 1.)
575  aa *= sum;
576  for (n = 1; n <= *nb; ++n) {
577  if (fabs(b[n]) < aa)
578  b[n] = 0.;
579  else
580  b[n] /= sum;
581  }
582  }
583 
584  }
585  else {
586  /* Error return -- X, NB, or ALPHA is out of range : */
587  b[1] = 0.;
588  *ncalc = min0(*nb,0) - 1;
589  }
590 }
sum
Type sum(Vector< Type > x)
Definition: convenience.hpp:33
License: GPL v2