| /* |
| * Mathlib : A C Library of Special Functions |
| * Copyright (C) 1998-2014 Ross Ihaka and the R Core team. |
| * |
| * This program is free software; you can redistribute it and/or modify |
| * it under the terms of the GNU General Public License as published by |
| * the Free Software Foundation; either version 2 of the License, or |
| * (at your option) any later version. |
| * |
| * This program is distributed in the hope that it will be useful, |
| * but WITHOUT ANY WARRANTY; without even the implied warranty of |
| * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the |
| * GNU General Public License for more details. |
| * |
| * You should have received a copy of the GNU General Public License |
| * along with this program; if not, a copy is available at |
| * https://www.R-project.org/Licenses/ |
| */ |
| |
| /* DESCRIPTION --> see below */ |
| |
| |
| /* From http://www.netlib.org/specfun/ribesl Fortran translated by f2c,... |
| * ------------------------------=#---- Martin Maechler, ETH Zurich |
| */ |
| #include "nmath.h" |
| #include "bessel.h" |
| |
| #ifndef MATHLIB_STANDALONE |
| #include <R_ext/Memory.h> |
| #endif |
| |
| #define min0(x, y) (((x) <= (y)) ? (x) : (y)) |
| |
| static void I_bessel(double *x, double *alpha, int *nb, |
| int *ize, double *bi, int *ncalc); |
| |
| /* .Internal(besselI(*)) : */ |
| double bessel_i(double x, double alpha, double expo) |
| { |
| int nb, ncalc, ize; |
| double na, *bi; |
| #ifndef MATHLIB_STANDALONE |
| const void *vmax; |
| #endif |
| |
| #ifdef IEEE_754 |
| /* NaNs propagated correctly */ |
| if (ISNAN(x) || ISNAN(alpha)) return x + alpha; |
| #endif |
| if (x < 0) { |
| ML_WARNING(ME_RANGE, "bessel_i"); |
| return ML_NAN; |
| } |
| ize = (int)expo; |
| na = floor(alpha); |
| if (alpha < 0) { |
| /* Using Abramowitz & Stegun 9.6.2 & 9.6.6 |
| * this may not be quite optimal (CPU and accuracy wise) */ |
| return(bessel_i(x, -alpha, expo) + |
| ((alpha == na) ? /* sin(pi * alpha) = 0 */ 0 : |
| bessel_k(x, -alpha, expo) * |
| ((ize == 1)? 2. : 2.*exp(-2.*x))/M_PI * sinpi(-alpha))); |
| } |
| nb = 1 + (int)na;/* nb-1 <= alpha < nb */ |
| alpha -= (double)(nb-1); |
| #ifdef MATHLIB_STANDALONE |
| bi = (double *) calloc(nb, sizeof(double)); |
| if (!bi) MATHLIB_ERROR("%s", _("bessel_i allocation error")); |
| #else |
| vmax = vmaxget(); |
| bi = (double *) R_alloc((size_t) nb, sizeof(double)); |
| #endif |
| I_bessel(&x, &alpha, &nb, &ize, bi, &ncalc); |
| if(ncalc != nb) {/* error input */ |
| if(ncalc < 0) |
| MATHLIB_WARNING4(_("bessel_i(%g): ncalc (=%d) != nb (=%d); alpha=%g. Arg. out of range?\n"), |
| x, ncalc, nb, alpha); |
| else |
| MATHLIB_WARNING2(_("bessel_i(%g,nu=%g): precision lost in result\n"), |
| x, alpha+(double)nb-1); |
| } |
| x = bi[nb-1]; |
| #ifdef MATHLIB_STANDALONE |
| free(bi); |
| #else |
| vmaxset(vmax); |
| #endif |
| return x; |
| } |
| |
| /* modified version of bessel_i that accepts a work array instead of |
| allocating one. */ |
| double bessel_i_ex(double x, double alpha, double expo, double *bi) |
| { |
| int nb, ncalc, ize; |
| double na; |
| |
| #ifdef IEEE_754 |
| /* NaNs propagated correctly */ |
| if (ISNAN(x) || ISNAN(alpha)) return x + alpha; |
| #endif |
| if (x < 0) { |
| ML_WARNING(ME_RANGE, "bessel_i"); |
| return ML_NAN; |
| } |
| ize = (int)expo; |
| na = floor(alpha); |
| if (alpha < 0) { |
| /* Using Abramowitz & Stegun 9.6.2 & 9.6.6 |
| * this may not be quite optimal (CPU and accuracy wise) */ |
| return(bessel_i_ex(x, -alpha, expo, bi) + |
| ((alpha == na) ? 0 : |
| bessel_k_ex(x, -alpha, expo, bi) * |
| ((ize == 1)? 2. : 2.*exp(-2.*x))/M_PI * sinpi(-alpha))); |
| } |
| nb = 1 + (int)na;/* nb-1 <= alpha < nb */ |
| alpha -= (double)(nb-1); |
| I_bessel(&x, &alpha, &nb, &ize, bi, &ncalc); |
| if(ncalc != nb) {/* error input */ |
| if(ncalc < 0) |
| MATHLIB_WARNING4(_("bessel_i(%g): ncalc (=%d) != nb (=%d); alpha=%g. Arg. out of range?\n"), |
| x, ncalc, nb, alpha); |
| else |
| MATHLIB_WARNING2(_("bessel_i(%g,nu=%g): precision lost in result\n"), |
| x, alpha+(double)nb-1); |
| } |
| x = bi[nb-1]; |
| return x; |
| } |
| |
| static void I_bessel(double *x, double *alpha, int *nb, |
| int *ize, double *bi, int *ncalc) |
| { |
| /* ------------------------------------------------------------------- |
| |
| This routine calculates Bessel functions I_(N+ALPHA) (X) |
| for non-negative argument X, and non-negative order N+ALPHA, |
| with or without exponential scaling. |
| |
| |
| Explanation of variables in the calling sequence |
| |
| X - Non-negative argument for which |
| I's or exponentially scaled I's (I*EXP(-X)) |
| are to be calculated. If I's are to be calculated, |
| X must be less than exparg_BESS (IZE=1) or xlrg_BESS_IJ (IZE=2), |
| (see bessel.h). |
| ALPHA - Fractional part of order for which |
| I's or exponentially scaled I's (I*EXP(-X)) are |
| to be calculated. 0 <= ALPHA < 1.0. |
| NB - Number of functions to be calculated, NB > 0. |
| The first function calculated is of order ALPHA, and the |
| last is of order (NB - 1 + ALPHA). |
| IZE - Type. IZE = 1 if unscaled I's are to be calculated, |
| = 2 if exponentially scaled I's are to be calculated. |
| BI - Output vector of length NB. If the routine |
| terminates normally (NCALC=NB), the vector BI contains the |
| functions I(ALPHA,X) through I(NB-1+ALPHA,X), or the |
| corresponding exponentially scaled functions. |
| NCALC - Output variable indicating possible errors. |
| Before using the vector BI, the user should check that |
| NCALC=NB, i.e., all orders have been calculated to |
| the desired accuracy. See error returns below. |
| |
| |
| ******************************************************************* |
| ******************************************************************* |
| |
| Error returns |
| |
| In case of an error, NCALC != NB, and not all I's are |
| calculated to the desired accuracy. |
| |
| NCALC < 0: An argument is out of range. For example, |
| NB <= 0, IZE is not 1 or 2, or IZE=1 and ABS(X) >= EXPARG_BESS. |
| In this case, the BI-vector is not calculated, and NCALC is |
| set to MIN0(NB,0)-1 so that NCALC != NB. |
| |
| NB > NCALC > 0: Not all requested function values could |
| be calculated accurately. This usually occurs because NB is |
| much larger than ABS(X). In this case, BI[N] is calculated |
| to the desired accuracy for N <= NCALC, but precision |
| is lost for NCALC < N <= NB. If BI[N] does not vanish |
| for N > NCALC (because it is too small to be represented), |
| and BI[N]/BI[NCALC] = 10**(-K), then only the first NSIG-K |
| significant figures of BI[N] can be trusted. |
| |
| |
| Intrinsic functions required are: |
| |
| DBLE, EXP, gamma_cody, INT, MAX, MIN, REAL, SQRT |
| |
| |
| Acknowledgement |
| |
| This program is based on a program written by David J. |
| Sookne (2) that computes values of the Bessel functions J or |
| I of float argument and long order. Modifications include |
| the restriction of the computation to the I Bessel function |
| of non-negative float argument, the extension of the computation |
| to arbitrary positive order, the inclusion of optional |
| exponential scaling, and the elimination of most underflow. |
| An earlier version was published in (3). |
| |
| References: "A Note on Backward Recurrence Algorithms," Olver, |
| F. W. J., and Sookne, D. J., Math. Comp. 26, 1972, |
| pp 941-947. |
| |
| "Bessel Functions of Real Argument and Integer Order," |
| Sookne, D. J., NBS Jour. of Res. B. 77B, 1973, pp |
| 125-132. |
| |
| "ALGORITHM 597, Sequence of Modified Bessel Functions |
| of the First Kind," Cody, W. J., Trans. Math. Soft., |
| 1983, pp. 242-245. |
| |
| Latest modification: May 30, 1989 |
| |
| Modified by: W. J. Cody and L. Stoltz |
| Applied Mathematics Division |
| Argonne National Laboratory |
| Argonne, IL 60439 |
| */ |
| |
| /*------------------------------------------------------------------- |
| Mathematical constants |
| -------------------------------------------------------------------*/ |
| const static double const__ = 1.585; |
| |
| /* Local variables */ |
| int nend, intx, nbmx, k, l, n, nstart; |
| double pold, test, p, em, en, empal, emp2al, halfx, |
| aa, bb, cc, psave, plast, tover, psavel, sum, nu, twonu; |
| |
| /*Parameter adjustments */ |
| --bi; |
| nu = *alpha; |
| twonu = nu + nu; |
| |
| /*------------------------------------------------------------------- |
| Check for X, NB, OR IZE out of range. |
| ------------------------------------------------------------------- */ |
| if (*nb > 0 && *x >= 0. && (0. <= nu && nu < 1.) && |
| (1 <= *ize && *ize <= 2) ) { |
| |
| *ncalc = *nb; |
| if(*ize == 1 && *x > exparg_BESS) { |
| for(k=1; k <= *nb; k++) |
| bi[k]=ML_POSINF; /* the limit *is* = Inf */ |
| return; |
| } |
| if(*ize == 2 && *x > xlrg_BESS_IJ) { |
| for(k=1; k <= *nb; k++) |
| bi[k]= 0.; /* The limit exp(-x) * I_nu(x) --> 0 : */ |
| return; |
| } |
| intx = (int) (*x);/* fine, since *x <= xlrg_BESS_IJ <<< LONG_MAX */ |
| if (*x >= rtnsig_BESS) { /* "non-small" x ( >= 1e-4 ) */ |
| /* ------------------------------------------------------------------- |
| Initialize the forward sweep, the P-sequence of Olver |
| ------------------------------------------------------------------- */ |
| nbmx = *nb - intx; |
| n = intx + 1; |
| en = (double) (n + n) + twonu; |
| plast = 1.; |
| p = en / *x; |
| /* ------------------------------------------------ |
| Calculate general significance test |
| ------------------------------------------------ */ |
| test = ensig_BESS + ensig_BESS; |
| if (intx << 1 > nsig_BESS * 5) { |
| test = sqrt(test * p); |
| } else { |
| test /= R_pow_di(const__, intx); |
| } |
| if (nbmx >= 3) { |
| /* -------------------------------------------------- |
| Calculate P-sequence until N = NB-1 |
| Check for possible overflow. |
| ------------------------------------------------ */ |
| tover = enten_BESS / ensig_BESS; |
| nstart = intx + 2; |
| nend = *nb - 1; |
| for (k = nstart; k <= nend; ++k) { |
| n = k; |
| en += 2.; |
| pold = plast; |
| plast = p; |
| p = en * plast / *x + pold; |
| if (p > tover) { |
| /* ------------------------------------------------ |
| To avoid overflow, divide P-sequence by TOVER. |
| Calculate P-sequence until ABS(P) > 1. |
| ---------------------------------------------- */ |
| tover = enten_BESS; |
| p /= tover; |
| plast /= tover; |
| psave = p; |
| psavel = plast; |
| nstart = n + 1; |
| do { |
| ++n; |
| en += 2.; |
| pold = plast; |
| plast = p; |
| p = en * plast / *x + pold; |
| } |
| while (p <= 1.); |
| |
| bb = en / *x; |
| /* ------------------------------------------------ |
| Calculate backward test, and find NCALC, |
| the highest N such that the test is passed. |
| ------------------------------------------------ */ |
| test = pold * plast / ensig_BESS; |
| test *= .5 - .5 / (bb * bb); |
| p = plast * tover; |
| --n; |
| en -= 2.; |
| nend = min0(*nb,n); |
| for (l = nstart; l <= nend; ++l) { |
| *ncalc = l; |
| pold = psavel; |
| psavel = psave; |
| psave = en * psavel / *x + pold; |
| if (psave * psavel > test) { |
| goto L90; |
| } |
| } |
| *ncalc = nend + 1; |
| L90: |
| --(*ncalc); |
| goto L120; |
| } |
| } |
| n = nend; |
| en = (double)(n + n) + twonu; |
| /*--------------------------------------------------- |
| Calculate special significance test for NBMX > 2. |
| --------------------------------------------------- */ |
| test = fmax2(test,sqrt(plast * ensig_BESS) * sqrt(p + p)); |
| } |
| /* -------------------------------------------------------- |
| Calculate P-sequence until significance test passed. |
| -------------------------------------------------------- */ |
| do { |
| ++n; |
| en += 2.; |
| pold = plast; |
| plast = p; |
| p = en * plast / *x + pold; |
| } while (p < test); |
| |
| L120: |
| /* ------------------------------------------------------------------- |
| Initialize the backward recursion and the normalization sum. |
| ------------------------------------------------------------------- */ |
| ++n; |
| en += 2.; |
| bb = 0.; |
| aa = 1. / p; |
| em = (double) n - 1.; |
| empal = em + nu; |
| emp2al = em - 1. + twonu; |
| sum = aa * empal * emp2al / em; |
| nend = n - *nb; |
| if (nend < 0) { |
| /* ----------------------------------------------------- |
| N < NB, so store BI[N] and set higher orders to 0.. |
| ----------------------------------------------------- */ |
| bi[n] = aa; |
| nend = -nend; |
| for (l = 1; l <= nend; ++l) { |
| bi[n + l] = 0.; |
| } |
| } else { |
| if (nend > 0) { |
| /* ----------------------------------------------------- |
| Recur backward via difference equation, |
| calculating (but not storing) BI[N], until N = NB. |
| --------------------------------------------------- */ |
| |
| for (l = 1; l <= nend; ++l) { |
| --n; |
| en -= 2.; |
| cc = bb; |
| bb = aa; |
| /* for x ~= 1500, sum would overflow to 'inf' here, |
| * and the final bi[] /= sum would give 0 wrongly; |
| * RE-normalize (aa, sum) here -- no need to undo */ |
| if(nend > 100 && aa > 1e200) { |
| /* multiply by 2^-900 = 1.18e-271 */ |
| cc = ldexp(cc, -900); |
| bb = ldexp(bb, -900); |
| sum = ldexp(sum,-900); |
| } |
| aa = en * bb / *x + cc; |
| em -= 1.; |
| emp2al -= 1.; |
| if (n == 1) { |
| break; |
| } |
| if (n == 2) { |
| emp2al = 1.; |
| } |
| empal -= 1.; |
| sum = (sum + aa * empal) * emp2al / em; |
| } |
| } |
| /* --------------------------------------------------- |
| Store BI[NB] |
| --------------------------------------------------- */ |
| bi[n] = aa; |
| if (*nb <= 1) { |
| sum = sum + sum + aa; |
| goto L230; |
| } |
| /* ------------------------------------------------- |
| Calculate and Store BI[NB-1] |
| ------------------------------------------------- */ |
| --n; |
| en -= 2.; |
| bi[n] = en * aa / *x + bb; |
| if (n == 1) { |
| goto L220; |
| } |
| em -= 1.; |
| if (n == 2) |
| emp2al = 1.; |
| else |
| emp2al -= 1.; |
| |
| empal -= 1.; |
| sum = (sum + bi[n] * empal) * emp2al / em; |
| } |
| nend = n - 2; |
| if (nend > 0) { |
| /* -------------------------------------------- |
| Calculate via difference equation |
| and store BI[N], until N = 2. |
| ------------------------------------------ */ |
| for (l = 1; l <= nend; ++l) { |
| --n; |
| en -= 2.; |
| bi[n] = en * bi[n + 1] / *x + bi[n + 2]; |
| em -= 1.; |
| if (n == 2) |
| emp2al = 1.; |
| else |
| emp2al -= 1.; |
| empal -= 1.; |
| sum = (sum + bi[n] * empal) * emp2al / em; |
| } |
| } |
| /* ---------------------------------------------- |
| Calculate BI[1] |
| -------------------------------------------- */ |
| bi[1] = 2. * empal * bi[2] / *x + bi[3]; |
| L220: |
| sum = sum + sum + bi[1]; |
| |
| L230: |
| /* --------------------------------------------------------- |
| Normalize. Divide all BI[N] by sum. |
| --------------------------------------------------------- */ |
| if (nu != 0.) |
| sum *= (Rf_gamma_cody(1. + nu) * pow(*x * .5, -nu)); |
| if (*ize == 1) |
| sum *= exp(-(*x)); |
| aa = enmten_BESS; |
| if (sum > 1.) |
| aa *= sum; |
| for (n = 1; n <= *nb; ++n) { |
| if (bi[n] < aa) |
| bi[n] = 0.; |
| else |
| bi[n] /= sum; |
| } |
| return; |
| } else { /* small x < 1e-4 */ |
| /* ----------------------------------------------------------- |
| Two-term ascending series for small X. |
| -----------------------------------------------------------*/ |
| aa = 1.; |
| empal = 1. + nu; |
| #ifdef IEEE_754 |
| /* No need to check for underflow */ |
| halfx = .5 * *x; |
| #else |
| if (*x > enmten_BESS) */ |
| halfx = .5 * *x; |
| else |
| halfx = 0.; |
| #endif |
| if (nu != 0.) |
| aa = pow(halfx, nu) / Rf_gamma_cody(empal); |
| if (*ize == 2) |
| aa *= exp(-(*x)); |
| bb = halfx * halfx; |
| bi[1] = aa + aa * bb / empal; |
| if (*x != 0. && bi[1] == 0.) |
| *ncalc = 0; |
| if (*nb > 1) { |
| if (*x == 0.) { |
| for (n = 2; n <= *nb; ++n) |
| bi[n] = 0.; |
| } else { |
| /* ------------------------------------------------- |
| Calculate higher-order functions. |
| ------------------------------------------------- */ |
| cc = halfx; |
| tover = (enmten_BESS + enmten_BESS) / *x; |
| if (bb != 0.) |
| tover = enmten_BESS / bb; |
| for (n = 2; n <= *nb; ++n) { |
| aa /= empal; |
| empal += 1.; |
| aa *= cc; |
| if (aa <= tover * empal) |
| bi[n] = aa = 0.; |
| else |
| bi[n] = aa + aa * bb / empal; |
| if (bi[n] == 0. && *ncalc > n) |
| *ncalc = n - 1; |
| } |
| } |
| } |
| } |
| } else { /* argument out of range */ |
| *ncalc = min0(*nb,0) - 1; |
| } |
| } |