TMB Documentation  v1.9.11
bessel_k.cpp
1 /*
2  * Mathlib : A C Library of Special Functions
3  * Copyright (C) 1998-2014 Ross Ihaka and the R Core team.
4  * Copyright (C) 2002-3 The R Foundation
5  *
6  * This program is free software; you can redistribute it and/or modify
7  * it under the terms of the GNU General Public License as published by
8  * the Free Software Foundation; either version 2 of the License, or
9  * (at your option) any later version.
10  *
11  * This program is distributed in the hope that it will be useful,
12  * but WITHOUT ANY WARRANTY; without even the implied warranty of
13  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14  * GNU General Public License for more details.
15  *
16  * You should have received a copy of the GNU General Public License
17  * along with this program; if not, a copy is available at
18  * https://www.R-project.org/Licenses/
19  */
20 
21 /* DESCRIPTION --> see below */
22 
23 
24 /* From http://www.netlib.org/specfun/rkbesl Fortran translated by f2c,...
25  * ------------------------------=#---- Martin Maechler, ETH Zurich
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 #define max0(x, y) (((x) <= (y)) ? (y) : (x))
36 
37 template<class Float>
38 static void K_bessel(Float *x, Float *alpha, int *nb,
39  int *ize, Float *bk, int *ncalc);
40 
41 template<class Float>
42 Float bessel_k(Float x, Float alpha, double expo)
43 {
44  int nb, ncalc, ize;
45  Float *bk;
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_k");
56  return ML_NAN;
57  }
58  ize = (int)expo;
59  if(alpha < 0)
60  alpha = -alpha;
61  nb = 1+ (int)floor(alpha);/* nb-1 <= |alpha| < nb */
62  alpha -= (double)(nb-1);
63 #ifdef MATHLIB_STANDALONE
64  bk = (Float *) calloc(nb, sizeof(Float));
65  if (!bk) MATHLIB_ERROR("%s", _("bessel_k allocation error"));
66 #else
67  vmax = vmaxget();
68  bk = (Float *) R_alloc((size_t) nb, sizeof(Float));
69 #endif
70  K_bessel(&x, &alpha, &nb, &ize, bk, &ncalc);
71  if(ncalc != nb) {/* error input */
72  if(ncalc < 0)
73  MATHLIB_WARNING4(_("bessel_k(%g): ncalc (=%d) != nb (=%d); alpha=%g. Arg. out of range?\n"),
74  x, ncalc, nb, alpha);
75  else
76  MATHLIB_WARNING2(_("bessel_k(%g,nu=%g): precision lost in result\n"),
77  x, alpha+(double)nb-1);
78  }
79  x = bk[nb-1];
80 #ifdef MATHLIB_STANDALONE
81  free(bk);
82 #else
83  vmaxset(vmax);
84 #endif
85  return x;
86 }
87 
88 /* modified version of bessel_k that accepts a work array instead of
89  allocating one. */
90 template<class Float>
91 Float bessel_k_ex(Float x, Float alpha, double expo, Float *bk)
92 {
93  int nb, ncalc, ize;
94 
95 #ifdef IEEE_754
96  /* NaNs propagated correctly */
97  if (ISNAN(x) || ISNAN(alpha)) return x + alpha;
98 #endif
99  if (x < 0) {
100  ML_ERROR(ME_RANGE, "bessel_k");
101  return ML_NAN;
102  }
103  ize = (int)expo;
104  if(alpha < 0)
105  alpha = -alpha;
106  nb = 1+ (int)floor(alpha);/* nb-1 <= |alpha| < nb */
107  alpha -= (double)(nb-1);
108  K_bessel(&x, &alpha, &nb, &ize, bk, &ncalc);
109  if(ncalc != nb) {/* error input */
110  if(ncalc < 0)
111  MATHLIB_WARNING4(_("bessel_k(%g): ncalc (=%d) != nb (=%d); alpha=%g. Arg. out of range?\n"),
112  x, ncalc, nb, alpha);
113  else
114  MATHLIB_WARNING2(_("bessel_k(%g,nu=%g): precision lost in result\n"),
115  x, alpha+(double)nb-1);
116  }
117  x = bk[nb-1];
118  return x;
119 }
120 
121 template<class Float>
122 static void K_bessel(Float *x, Float *alpha, int *nb,
123  int *ize, Float *bk, int *ncalc)
124 {
125 /*-------------------------------------------------------------------
126 
127  This routine calculates modified Bessel functions
128  of the third kind, K_(N+ALPHA) (X), for non-negative
129  argument X, and non-negative order N+ALPHA, with or without
130  exponential scaling.
131 
132  Explanation of variables in the calling sequence
133 
134  X - Non-negative argument for which
135  K's or exponentially scaled K's (K*EXP(X))
136  are to be calculated. If K's are to be calculated,
137  X must not be greater than XMAX_BESS_K.
138  ALPHA - Fractional part of order for which
139  K's or exponentially scaled K's (K*EXP(X)) are
140  to be calculated. 0 <= ALPHA < 1.0.
141  NB - Number of functions to be calculated, NB > 0.
142  The first function calculated is of order ALPHA, and the
143  last is of order (NB - 1 + ALPHA).
144  IZE - Type. IZE = 1 if unscaled K's are to be calculated,
145  = 2 if exponentially scaled K's are to be calculated.
146  BK - Output vector of length NB. If the
147  routine terminates normally (NCALC=NB), the vector BK
148  contains the functions K(ALPHA,X), ... , K(NB-1+ALPHA,X),
149  or the corresponding exponentially scaled functions.
150  If (0 < NCALC < NB), BK(I) contains correct function
151  values for I <= NCALC, and contains the ratios
152  K(ALPHA+I-1,X)/K(ALPHA+I-2,X) for the rest of the array.
153  NCALC - Output variable indicating possible errors.
154  Before using the vector BK, the user should check that
155  NCALC=NB, i.e., all orders have been calculated to
156  the desired accuracy. See error returns below.
157 
158 
159  *******************************************************************
160 
161  Error returns
162 
163  In case of an error, NCALC != NB, and not all K's are
164  calculated to the desired accuracy.
165 
166  NCALC < -1: An argument is out of range. For example,
167  NB <= 0, IZE is not 1 or 2, or IZE=1 and ABS(X) >= XMAX_BESS_K.
168  In this case, the B-vector is not calculated,
169  and NCALC is set to MIN0(NB,0)-2 so that NCALC != NB.
170  NCALC = -1: Either K(ALPHA,X) >= XINF or
171  K(ALPHA+NB-1,X)/K(ALPHA+NB-2,X) >= XINF. In this case,
172  the B-vector is not calculated. Note that again
173  NCALC != NB.
174 
175  0 < NCALC < NB: Not all requested function values could
176  be calculated accurately. BK(I) contains correct function
177  values for I <= NCALC, and contains the ratios
178  K(ALPHA+I-1,X)/K(ALPHA+I-2,X) for the rest of the array.
179 
180 
181  Intrinsic functions required are:
182 
183  ABS, AINT, EXP, INT, LOG, MAX, MIN, SINH, SQRT
184 
185 
186  Acknowledgement
187 
188  This program is based on a program written by J. B. Campbell
189  (2) that computes values of the Bessel functions K of float
190  argument and float order. Modifications include the addition
191  of non-scaled functions, parameterization of machine
192  dependencies, and the use of more accurate approximations
193  for SINH and SIN.
194 
195  References: "On Temme's Algorithm for the Modified Bessel
196  Functions of the Third Kind," Campbell, J. B.,
197  TOMS 6(4), Dec. 1980, pp. 581-586.
198 
199  "A FORTRAN IV Subroutine for the Modified Bessel
200  Functions of the Third Kind of Real Order and Real
201  Argument," Campbell, J. B., Report NRC/ERB-925,
202  National Research Council, Canada.
203 
204  Latest modification: May 30, 1989
205 
206  Modified by: W. J. Cody and L. Stoltz
207  Applied Mathematics Division
208  Argonne National Laboratory
209  Argonne, IL 60439
210 
211  -------------------------------------------------------------------
212 */
213  /*---------------------------------------------------------------------
214  * Mathematical constants
215  * A = LOG(2) - Euler's constant
216  * D = SQRT(2/PI)
217  ---------------------------------------------------------------------*/
218  const static double a = .11593151565841244881;
219 
220  /*---------------------------------------------------------------------
221  P, Q - Approximation for LOG(GAMMA(1+ALPHA))/ALPHA + Euler's constant
222  Coefficients converted from hex to decimal and modified
223  by W. J. Cody, 2/26/82 */
224  const static double p[8] = { .805629875690432845,20.4045500205365151,
225  157.705605106676174,536.671116469207504,900.382759291288778,
226  730.923886650660393,229.299301509425145,.822467033424113231 };
227  const static double q[7] = { 29.4601986247850434,277.577868510221208,
228  1206.70325591027438,2762.91444159791519,3443.74050506564618,
229  2210.63190113378647,572.267338359892221 };
230  /* R, S - Approximation for (1-ALPHA*PI/SIN(ALPHA*PI))/(2.D0*ALPHA) */
231  const static double r[5] = { -.48672575865218401848,13.079485869097804016,
232  -101.96490580880537526,347.65409106507813131,
233  3.495898124521934782e-4 };
234  const static double s[4] = { -25.579105509976461286,212.57260432226544008,
235  -610.69018684944109624,422.69668805777760407 };
236  /* T - Approximation for SINH(Y)/Y */
237  const static double t[6] = { 1.6125990452916363814e-10,
238  2.5051878502858255354e-8,2.7557319615147964774e-6,
239  1.9841269840928373686e-4,.0083333333333334751799,
240  .16666666666666666446 };
241  /*---------------------------------------------------------------------*/
242  const static double estm[6] = { 52.0583,5.7607,2.7782,14.4303,185.3004, 9.3715 };
243  const static double estf[7] = { 41.8341,7.1075,6.4306,42.511,1.35633,84.5096,20.};
244 
245  /* Local variables */
246  int iend, i, j, k, m, ii, mplus1;
247  Float x2by4, twox, c, blpha, ratio, wminf;
248  Float d1, d2, d3, f0, f1, f2, p0, q0, t1, t2, twonu;
249  Float dm, ex, bk1, bk2, nu;
250 
251  ii = 0; /* -Wall */
252 
253  ex = *x;
254  nu = *alpha;
255  *ncalc = min0(*nb,0) - 2;
256  if (*nb > 0 && (0. <= nu && nu < 1.) && (1 <= *ize && *ize <= 2)) {
257  if(ex <= 0 || (*ize == 1 && ex > xmax_BESS_K)) {
258  if(ex <= 0) {
259  if(ex < 0) ML_ERROR(ME_RANGE, "K_bessel");
260  for(i=0; i < *nb; i++)
261  bk[i] = ML_POSINF;
262  } else /* would only have underflow */
263  for(i=0; i < *nb; i++)
264  bk[i] = 0.;
265  *ncalc = *nb;
266  return;
267  }
268  k = 0;
269  if (nu < sqxmin_BESS_K) {
270  // Kasper: Don't clear the derivatives of nu !
271  // nu = 0.;
272  } else if (nu > .5) {
273  k = 1;
274  nu -= 1.;
275  }
276  twonu = nu + nu;
277  iend = *nb + k - 1;
278  c = nu * nu;
279  d3 = -c;
280  if (ex <= 1.) {
281  /* ------------------------------------------------------------
282  Calculation of P0 = GAMMA(1+ALPHA) * (2/X)**ALPHA
283  Q0 = GAMMA(1-ALPHA) * (X/2)**ALPHA
284  ------------------------------------------------------------ */
285  d1 = 0.; d2 = p[0];
286  t1 = 1.; t2 = q[0];
287  for (i = 2; i <= 7; i += 2) {
288  d1 = c * d1 + p[i - 1];
289  d2 = c * d2 + p[i];
290  t1 = c * t1 + q[i - 1];
291  t2 = c * t2 + q[i];
292  }
293  d1 = nu * d1;
294  t1 = nu * t1;
295  f1 = log(ex);
296  f0 = a + nu * (p[7] - nu * (d1 + d2) / (t1 + t2)) - f1;
297  q0 = exp(-nu * (a - nu * (p[7] + nu * (d1-d2) / (t1-t2)) - f1));
298  f1 = nu * f0;
299  p0 = exp(f1);
300  /* -----------------------------------------------------------
301  Calculation of F0 =
302  ----------------------------------------------------------- */
303  d1 = r[4];
304  t1 = 1.;
305  for (i = 0; i < 4; ++i) {
306  d1 = c * d1 + r[i];
307  t1 = c * t1 + s[i];
308  }
309  /* d2 := sinh(f1)/ nu = sinh(f1)/(f1/f0)
310  * = f0 * sinh(f1)/f1 */
311  if (fabs(f1) <= .5) {
312  f1 *= f1;
313  d2 = 0.;
314  for (i = 0; i < 6; ++i) {
315  d2 = f1 * d2 + t[i];
316  }
317  d2 = f0 + f0 * f1 * d2;
318  } else {
319  d2 = sinh(f1) / nu;
320  }
321  f0 = d2 - nu * d1 / (t1 * p0);
322  if (ex <= 1e-10) {
323  /* ---------------------------------------------------------
324  X <= 1.0E-10
325  Calculation of K(ALPHA,X) and X*K(ALPHA+1,X)/K(ALPHA,X)
326  --------------------------------------------------------- */
327  bk[0] = f0 + ex * f0;
328  if (*ize == 1) {
329  bk[0] -= ex * bk[0];
330  }
331  ratio = p0 / f0;
332  c = ex * DBL_MAX;
333  if (k != 0) {
334  /* ---------------------------------------------------
335  Calculation of K(ALPHA,X)
336  and X*K(ALPHA+1,X)/K(ALPHA,X), ALPHA >= 1/2
337  --------------------------------------------------- */
338  *ncalc = -1;
339  if (bk[0] >= c / ratio) {
340  return;
341  }
342  bk[0] = ratio * bk[0] / ex;
343  twonu += 2.;
344  ratio = twonu;
345  }
346  *ncalc = 1;
347  if (*nb == 1)
348  return;
349 
350  /* -----------------------------------------------------
351  Calculate K(ALPHA+L,X)/K(ALPHA+L-1,X),
352  L = 1, 2, ... , NB-1
353  ----------------------------------------------------- */
354  *ncalc = -1;
355  for (i = 1; i < *nb; ++i) {
356  if (ratio >= c)
357  return;
358 
359  bk[i] = ratio / ex;
360  twonu += 2.;
361  ratio = twonu;
362  }
363  *ncalc = 1;
364  goto L420;
365  } else {
366  /* ------------------------------------------------------
367  10^-10 < X <= 1.0
368  ------------------------------------------------------ */
369  c = 1.;
370  x2by4 = ex * ex / 4.;
371  p0 = .5 * p0;
372  q0 = .5 * q0;
373  d1 = -1.;
374  d2 = 0.;
375  bk1 = 0.;
376  bk2 = 0.;
377  f1 = f0;
378  f2 = p0;
379  do {
380  d1 += 2.;
381  d2 += 1.;
382  d3 = d1 + d3;
383  c = x2by4 * c / d2;
384  f0 = (d2 * f0 + p0 + q0) / d3;
385  p0 /= d2 - nu;
386  q0 /= d2 + nu;
387  t1 = c * f0;
388  t2 = c * (p0 - d2 * f0);
389  bk1 += t1;
390  bk2 += t2;
391  } while (fabs(t1 / (f1 + bk1)) > DBL_EPSILON ||
392  fabs(t2 / (f2 + bk2)) > DBL_EPSILON);
393  bk1 = f1 + bk1;
394  bk2 = 2. * (f2 + bk2) / ex;
395  if (*ize == 2) {
396  d1 = exp(ex);
397  bk1 *= d1;
398  bk2 *= d1;
399  }
400  wminf = estf[0] * ex + estf[1];
401  }
402  } else if (DBL_EPSILON * ex > 1.) {
403  /* -------------------------------------------------
404  X > 1./EPS
405  ------------------------------------------------- */
406  *ncalc = *nb;
407  bk1 = 1. / (M_SQRT_2dPI * sqrt(ex));
408  for (i = 0; i < *nb; ++i)
409  bk[i] = bk1;
410  return;
411 
412  } else {
413  /* -------------------------------------------------------
414  X > 1.0
415  ------------------------------------------------------- */
416  twox = ex + ex;
417  blpha = 0.;
418  ratio = 0.;
419  if (ex <= 4.) {
420  /* ----------------------------------------------------------
421  Calculation of K(ALPHA+1,X)/K(ALPHA,X), 1.0 <= X <= 4.0
422  ----------------------------------------------------------*/
423  d2 = trunc(estm[0] / ex + estm[1]);
424  m = (int) trunc(d2);
425  d1 = d2 + d2;
426  d2 -= .5;
427  d2 *= d2;
428  for (i = 2; i <= m; ++i) {
429  d1 -= 2.;
430  d2 -= d1;
431  ratio = (d3 + d2) / (twox + d1 - ratio);
432  }
433  /* -----------------------------------------------------------
434  Calculation of I(|ALPHA|,X) and I(|ALPHA|+1,X) by backward
435  recurrence and K(ALPHA,X) from the wronskian
436  -----------------------------------------------------------*/
437  d2 = trunc(estm[2] * ex + estm[3]);
438  m = (int) trunc(d2);
439  // Kasper: No need for fabs here !
440  // c = fabs(nu); |
441  c = nu; // <------------+
442  d3 = c + c;
443  d1 = d3 - 1.;
444  f1 = DBL_MIN;
445  f0 = (2. * (c + d2) / ex + .5 * ex / (c + d2 + 1.)) * DBL_MIN;
446  for (i = 3; i <= m; ++i) {
447  d2 -= 1.;
448  f2 = (d3 + d2 + d2) * f0;
449  blpha = (1. + d1 / d2) * (f2 + blpha);
450  f2 = f2 / ex + f1;
451  f1 = f0;
452  f0 = f2;
453  }
454  f1 = (d3 + 2.) * f0 / ex + f1;
455  d1 = 0.;
456  t1 = 1.;
457  for (i = 1; i <= 7; ++i) {
458  d1 = c * d1 + p[i - 1];
459  t1 = c * t1 + q[i - 1];
460  }
461  p0 = exp(c * (a + c * (p[7] - c * d1 / t1) - log(ex))) / ex;
462  f2 = (c + .5 - ratio) * f1 / ex;
463  bk1 = p0 + (d3 * f0 - f2 + f0 + blpha) / (f2 + f1 + f0) * p0;
464  if (*ize == 1) {
465  bk1 *= exp(-ex);
466  }
467  wminf = estf[2] * ex + estf[3];
468  } else {
469  /* ---------------------------------------------------------
470  Calculation of K(ALPHA,X) and K(ALPHA+1,X)/K(ALPHA,X), by
471  backward recurrence, for X > 4.0
472  ----------------------------------------------------------*/
473  dm = trunc(estm[4] / ex + estm[5]);
474  m = (int) trunc(dm);
475  d2 = dm - .5;
476  d2 *= d2;
477  d1 = dm + dm;
478  for (i = 2; i <= m; ++i) {
479  dm -= 1.;
480  d1 -= 2.;
481  d2 -= d1;
482  ratio = (d3 + d2) / (twox + d1 - ratio);
483  blpha = (ratio + ratio * blpha) / dm;
484  }
485  bk1 = 1. / ((M_SQRT_2dPI + M_SQRT_2dPI * blpha) * sqrt(ex));
486  if (*ize == 1)
487  bk1 *= exp(-ex);
488  wminf = estf[4] * (ex - fabs(ex - estf[6])) + estf[5];
489  }
490  /* ---------------------------------------------------------
491  Calculation of K(ALPHA+1,X)
492  from K(ALPHA,X) and K(ALPHA+1,X)/K(ALPHA,X)
493  --------------------------------------------------------- */
494  bk2 = bk1 + bk1 * (nu + .5 - ratio) / ex;
495  }
496  /*--------------------------------------------------------------------
497  Calculation of 'NCALC', K(ALPHA+I,X), I = 0, 1, ... , NCALC-1,
498  & K(ALPHA+I,X)/K(ALPHA+I-1,X), I = NCALC, NCALC+1, ... , NB-1
499  -------------------------------------------------------------------*/
500  *ncalc = *nb;
501  bk[0] = bk1;
502  if (iend == 0)
503  return;
504 
505  j = 1 - k;
506  if (j >= 0)
507  bk[j] = bk2;
508 
509  if (iend == 1)
510  return;
511 
512  m = min0((int) trunc(wminf - nu),iend);
513  for (i = 2; i <= m; ++i) {
514  t1 = bk1;
515  bk1 = bk2;
516  twonu += 2.;
517  if (ex < 1.) {
518  if (bk1 >= DBL_MAX / twonu * ex)
519  break;
520  } else {
521  if (bk1 / ex >= DBL_MAX / twonu)
522  break;
523  }
524  bk2 = twonu / ex * bk1 + t1;
525  ii = i;
526  ++j;
527  if (j >= 0) {
528  bk[j] = bk2;
529  }
530  }
531 
532  m = ii;
533  if (m == iend) {
534  return;
535  }
536  ratio = bk2 / bk1;
537  mplus1 = m + 1;
538  *ncalc = -1;
539  for (i = mplus1; i <= iend; ++i) {
540  twonu += 2.;
541  ratio = twonu / ex + 1./ratio;
542  ++j;
543  if (j >= 1) {
544  bk[j] = ratio;
545  } else {
546  if (bk2 >= DBL_MAX / ratio)
547  return;
548 
549  bk2 *= ratio;
550  }
551  }
552  *ncalc = max0(1, mplus1 - k);
553  if (*ncalc == 1)
554  bk[0] = bk2;
555  if (*nb == 1)
556  return;
557 
558 L420:
559  for (i = *ncalc; i < *nb; ++i) { /* i == *ncalc */
560 #ifndef IEEE_754
561  if (bk[i-1] >= DBL_MAX / bk[i])
562  return;
563 #endif
564  bk[i] *= bk[i-1];
565  (*ncalc)++;
566  }
567  }
568 }
License: GPL v2