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