/* iv.c * * Modified Bessel function of noninteger order * * * * SYNOPSIS: * * double v, x, y, iv(); * * y = iv( v, x ); * * * * DESCRIPTION: * * Returns modified Bessel function of order v of the * argument. If x is negative, v must be integer valued. * */ /* iv.c */ /* Modified Bessel function of noninteger order */ /* If x < 0, then v must be an integer. */ /* * Parts of the code are copyright: * * Cephes Math Library Release 2.8: June, 2000 * Copyright 1984, 1987, 1988, 2000 by Stephen L. Moshier * * And other parts: * * Copyright (c) 2006 Xiaogang Zhang * Use, modification and distribution are subject to the * Boost Software License, Version 1.0. * * Boost Software License - Version 1.0 - August 17th, 2003 * * Permission is hereby granted, free of charge, to any person or * organization obtaining a copy of the software and accompanying * documentation covered by this license (the "Software") to use, reproduce, * display, distribute, execute, and transmit the Software, and to prepare * derivative works of the Software, and to permit third-parties to whom the * Software is furnished to do so, all subject to the following: * * The copyright notices in the Software and this entire statement, * including the above license grant, this restriction and the following * disclaimer, must be included in all copies of the Software, in whole or * in part, and all derivative works of the Software, unless such copies or * derivative works are solely in the form of machine-executable object code * generated by a source language processor. * * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF * MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE, TITLE AND * NON-INFRINGEMENT. IN NO EVENT SHALL THE COPYRIGHT HOLDERS OR ANYONE * DISTRIBUTING THE SOFTWARE BE LIABLE FOR ANY DAMAGES OR OTHER LIABILITY, * WHETHER IN CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN * CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE * SOFTWARE. * * And the rest are: * * Copyright (C) 2009 Pauli Virtanen * Distributed under the same license as Scipy. * */ #include "mconf.h" #include #include extern double MACHEP; static double iv_asymptotic(double v, double x); static void ikv_asymptotic_uniform(double v, double x, double *Iv, double *Kv); static void ikv_temme(double v, double x, double *Iv, double *Kv); double iv(double v, double x) { int sign; double t, ax, res; if (npy_isnan(v) || npy_isnan(x)) { return NPY_NAN; } /* If v is a negative integer, invoke symmetry */ t = floor(v); if (v < 0.0) { if (t == v) { v = -v; /* symmetry */ t = -t; } } /* If x is negative, require v to be an integer */ sign = 1; if (x < 0.0) { if (t != v) { sf_error("iv", SF_ERROR_DOMAIN, NULL); return (NPY_NAN); } if (v != 2.0 * floor(v / 2.0)) { sign = -1; } } /* Avoid logarithm singularity */ if (x == 0.0) { if (v == 0.0) { return 1.0; } if (v < 0.0) { sf_error("iv", SF_ERROR_OVERFLOW, NULL); return NPY_INFINITY; } else return 0.0; } ax = fabs(x); if (fabs(v) > 50) { /* * Uniform asymptotic expansion for large orders. * * This appears to overflow slightly later than the Boost * implementation of Temme's method. */ ikv_asymptotic_uniform(v, ax, &res, NULL); } else { /* Otherwise: Temme's method */ ikv_temme(v, ax, &res, NULL); } res *= sign; return res; } /* * Compute Iv from (AMS5 9.7.1), asymptotic expansion for large |z| * Iv ~ exp(x)/sqrt(2 pi x) ( 1 + (4*v*v-1)/8x + (4*v*v-1)(4*v*v-9)/8x/2! + ...) */ static double iv_asymptotic(double v, double x) { double mu; double sum, term, prefactor, factor; int k; prefactor = exp(x) / sqrt(2 * NPY_PI * x); if (prefactor == NPY_INFINITY) { return prefactor; } mu = 4 * v * v; sum = 1.0; term = 1.0; k = 1; do { factor = (mu - (2 * k - 1) * (2 * k - 1)) / (8 * x) / k; if (k > 100) { /* didn't converge */ sf_error("iv(iv_asymptotic)", SF_ERROR_NO_RESULT, NULL); break; } term *= -factor; sum += term; ++k; } while (fabs(term) > MACHEP * fabs(sum)); return sum * prefactor; } /* * Uniform asymptotic expansion factors, (AMS5 9.3.9; AMS5 9.3.10) * * Computed with: * -------------------- import numpy as np t = np.poly1d([1,0]) def up1(p): return .5*t*t*(1-t*t)*p.deriv() + 1/8. * ((1-5*t*t)*p).integ() us = [np.poly1d([1])] for k in range(10): us.append(up1(us[-1])) n = us[-1].order for p in us: print "{" + ", ".join(["0"]*(n-p.order) + map(repr, p)) + "}," print "N_UFACTORS", len(us) print "N_UFACTOR_TERMS", us[-1].order + 1 * -------------------- */ #define N_UFACTORS 11 #define N_UFACTOR_TERMS 31 static const double asymptotic_ufactors[N_UFACTORS][N_UFACTOR_TERMS] = { {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1}, {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -0.20833333333333334, 0.0, 0.125, 0.0}, {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.3342013888888889, 0.0, -0.40104166666666669, 0.0, 0.0703125, 0.0, 0.0}, {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -1.0258125964506173, 0.0, 1.8464626736111112, 0.0, -0.89121093750000002, 0.0, 0.0732421875, 0.0, 0.0, 0.0}, {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 4.6695844234262474, 0.0, -11.207002616222995, 0.0, 8.78912353515625, 0.0, -2.3640869140624998, 0.0, 0.112152099609375, 0.0, 0.0, 0.0, 0.0}, {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -28.212072558200244, 0.0, 84.636217674600744, 0.0, -91.818241543240035, 0.0, 42.534998745388457, 0.0, -7.3687943594796312, 0.0, 0.22710800170898438, 0.0, 0.0, 0.0, 0.0, 0.0}, {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 212.5701300392171, 0.0, -765.25246814118157, 0.0, 1059.9904525279999, 0.0, -699.57962737613275, 0.0, 218.19051174421159, 0.0, -26.491430486951554, 0.0, 0.57250142097473145, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0}, {0, 0, 0, 0, 0, 0, 0, 0, 0, -1919.4576623184068, 0.0, 8061.7221817373083, 0.0, -13586.550006434136, 0.0, 11655.393336864536, 0.0, -5305.6469786134048, 0.0, 1200.9029132163525, 0.0, -108.09091978839464, 0.0, 1.7277275025844574, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0}, {0, 0, 0, 0, 0, 0, 20204.291330966149, 0.0, -96980.598388637503, 0.0, 192547.0012325315, 0.0, -203400.17728041555, 0.0, 122200.46498301747, 0.0, -41192.654968897557, 0.0, 7109.5143024893641, 0.0, -493.915304773088, 0.0, 6.074042001273483, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0}, {0, 0, 0, -242919.18790055133, 0.0, 1311763.6146629769, 0.0, -2998015.9185381061, 0.0, 3763271.2976564039, 0.0, -2813563.2265865342, 0.0, 1268365.2733216248, 0.0, -331645.17248456361, 0.0, 45218.768981362737, 0.0, -2499.8304818112092, 0.0, 24.380529699556064, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0}, {3284469.8530720375, 0.0, -19706819.11843222, 0.0, 50952602.492664628, 0.0, -74105148.211532637, 0.0, 66344512.274729028, 0.0, -37567176.660763353, 0.0, 13288767.166421819, 0.0, -2785618.1280864552, 0.0, 308186.40461266245, 0.0, -13886.089753717039, 0.0, 110.01714026924674, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0} }; /* * Compute Iv, Kv from (AMS5 9.7.7 + 9.7.8), asymptotic expansion for large v */ static void ikv_asymptotic_uniform(double v, double x, double *i_value, double *k_value) { double i_prefactor, k_prefactor; double t, t2, eta, z; double i_sum, k_sum, term, divisor; int k, n; int sign = 1; if (v < 0) { /* Negative v; compute I_{-v} and K_{-v} and use (AMS 9.6.2) */ sign = -1; v = -v; } z = x / v; t = 1 / sqrt(1 + z * z); t2 = t * t; eta = sqrt(1 + z * z) + log(z / (1 + 1 / t)); i_prefactor = sqrt(t / (2 * NPY_PI * v)) * exp(v * eta); i_sum = 1.0; k_prefactor = sqrt(NPY_PI * t / (2 * v)) * exp(-v * eta); k_sum = 1.0; divisor = v; for (n = 1; n < N_UFACTORS; ++n) { /* * Evaluate u_k(t) with Horner's scheme; * (using the knowledge about which coefficients are zero) */ term = 0; for (k = N_UFACTOR_TERMS - 1 - 3 * n; k < N_UFACTOR_TERMS - n; k += 2) { term *= t2; term += asymptotic_ufactors[n][k]; } for (k = 1; k < n; k += 2) { term *= t2; } if (n % 2 == 1) { term *= t; } /* Sum terms */ term /= divisor; i_sum += term; k_sum += (n % 2 == 0) ? term : -term; /* Check convergence */ if (fabs(term) < MACHEP) { break; } divisor *= v; } if (fabs(term) > 1e-3 * fabs(i_sum)) { /* Didn't converge */ sf_error("ikv_asymptotic_uniform", SF_ERROR_NO_RESULT, NULL); } if (fabs(term) > MACHEP * fabs(i_sum)) { /* Some precision lost */ sf_error("ikv_asymptotic_uniform", SF_ERROR_LOSS, NULL); } if (k_value != NULL) { /* symmetric in v */ *k_value = k_prefactor * k_sum; } if (i_value != NULL) { if (sign == 1) { *i_value = i_prefactor * i_sum; } else { /* (AMS 9.6.2) */ *i_value = (i_prefactor * i_sum + (2 / NPY_PI) * sin(NPY_PI * v) * k_prefactor * k_sum); } } } /* * The following code originates from the Boost C++ library, * from file `boost/math/special_functions/detail/bessel_ik.hpp`, * converted from C++ to C. */ #ifdef DEBUG #define BOOST_ASSERT(a) assert(a) #else #define BOOST_ASSERT(a) #endif /* * Modified Bessel functions of the first and second kind of fractional order * * Calculate K(v, x) and K(v+1, x) by method analogous to * Temme, Journal of Computational Physics, vol 21, 343 (1976) */ static int temme_ik_series(double v, double x, double *K, double *K1) { double f, h, p, q, coef, sum, sum1, tolerance; double a, b, c, d, sigma, gamma1, gamma2; unsigned long k; double gp; double gm; /* * |x| <= 2, Temme series converge rapidly * |x| > 2, the larger the |x|, the slower the convergence */ BOOST_ASSERT(fabs(x) <= 2); BOOST_ASSERT(fabs(v) <= 0.5f); gp = gamma(v + 1) - 1; gm = gamma(-v + 1) - 1; a = log(x / 2); b = exp(v * a); sigma = -a * v; c = fabs(v) < MACHEP ? 1 : sin(NPY_PI * v) / (v * NPY_PI); d = fabs(sigma) < MACHEP ? 1 : sinh(sigma) / sigma; gamma1 = fabs(v) < MACHEP ? -NPY_EULER : (0.5f / v) * (gp - gm) * c; gamma2 = (2 + gp + gm) * c / 2; /* initial values */ p = (gp + 1) / (2 * b); q = (1 + gm) * b / 2; f = (cosh(sigma) * gamma1 + d * (-a) * gamma2) / c; h = p; coef = 1; sum = coef * f; sum1 = coef * h; /* series summation */ tolerance = MACHEP; for (k = 1; k < MAXITER; k++) { f = (k * f + p + q) / (k * k - v * v); p /= k - v; q /= k + v; h = p - k * f; coef *= x * x / (4 * k); sum += coef * f; sum1 += coef * h; if (fabs(coef * f) < fabs(sum) * tolerance) { break; } } if (k == MAXITER) { sf_error("ikv_temme(temme_ik_series)", SF_ERROR_NO_RESULT, NULL); } *K = sum; *K1 = 2 * sum1 / x; return 0; } /* Evaluate continued fraction fv = I_(v+1) / I_v, derived from * Abramowitz and Stegun, Handbook of Mathematical Functions, 1972, 9.1.73 */ static int CF1_ik(double v, double x, double *fv) { double C, D, f, a, b, delta, tiny, tolerance; unsigned long k; /* * |x| <= |v|, CF1_ik converges rapidly * |x| > |v|, CF1_ik needs O(|x|) iterations to converge */ /* * modified Lentz's method, see * Lentz, Applied Optics, vol 15, 668 (1976) */ tolerance = 2 * MACHEP; tiny = 1 / sqrt(DBL_MAX); C = f = tiny; /* b0 = 0, replace with tiny */ D = 0; for (k = 1; k < MAXITER; k++) { a = 1; b = 2 * (v + k) / x; C = b + a / C; D = b + a * D; if (C == 0) { C = tiny; } if (D == 0) { D = tiny; } D = 1 / D; delta = C * D; f *= delta; if (fabs(delta - 1) <= tolerance) { break; } } if (k == MAXITER) { sf_error("ikv_temme(CF1_ik)", SF_ERROR_NO_RESULT, NULL); } *fv = f; return 0; } /* * Calculate K(v, x) and K(v+1, x) by evaluating continued fraction * z1 / z0 = U(v+1.5, 2v+1, 2x) / U(v+0.5, 2v+1, 2x), see * Thompson and Barnett, Computer Physics Communications, vol 47, 245 (1987) */ static int CF2_ik(double v, double x, double *Kv, double *Kv1) { double S, C, Q, D, f, a, b, q, delta, tolerance, current, prev; unsigned long k; /* * |x| >= |v|, CF2_ik converges rapidly * |x| -> 0, CF2_ik fails to converge */ BOOST_ASSERT(fabs(x) > 1); /* * Steed's algorithm, see Thompson and Barnett, * Journal of Computational Physics, vol 64, 490 (1986) */ tolerance = MACHEP; a = v * v - 0.25f; b = 2 * (x + 1); /* b1 */ D = 1 / b; /* D1 = 1 / b1 */ f = delta = D; /* f1 = delta1 = D1, coincidence */ prev = 0; /* q0 */ current = 1; /* q1 */ Q = C = -a; /* Q1 = C1 because q1 = 1 */ S = 1 + Q * delta; /* S1 */ for (k = 2; k < MAXITER; k++) { /* starting from 2 */ /* continued fraction f = z1 / z0 */ a -= 2 * (k - 1); b += 2; D = 1 / (b + a * D); delta *= b * D - 1; f += delta; /* series summation S = 1 + \sum_{n=1}^{\infty} C_n * z_n / z_0 */ q = (prev - (b - 2) * current) / a; prev = current; current = q; /* forward recurrence for q */ C *= -a / k; Q += C * q; S += Q * delta; /* S converges slower than f */ if (fabs(Q * delta) < fabs(S) * tolerance) { break; } } if (k == MAXITER) { sf_error("ikv_temme(CF2_ik)", SF_ERROR_NO_RESULT, NULL); } *Kv = sqrt(NPY_PI / (2 * x)) * exp(-x) / S; *Kv1 = *Kv * (0.5f + v + x + (v * v - 0.25f) * f) / x; return 0; } /* Flags for what to compute */ enum { need_i = 0x1, need_k = 0x2 }; /* * Compute I(v, x) and K(v, x) simultaneously by Temme's method, see * Temme, Journal of Computational Physics, vol 19, 324 (1975) */ static void ikv_temme(double v, double x, double *Iv_p, double *Kv_p) { /* Kv1 = K_(v+1), fv = I_(v+1) / I_v */ /* Ku1 = K_(u+1), fu = I_(u+1) / I_u */ double u, Iv, Kv, Kv1, Ku, Ku1, fv; double W, current, prev, next; int reflect = 0; unsigned n, k; int kind; kind = 0; if (Iv_p != NULL) { kind |= need_i; } if (Kv_p != NULL) { kind |= need_k; } if (v < 0) { reflect = 1; v = -v; /* v is non-negative from here */ kind |= need_k; } n = round(v); u = v - n; /* -1/2 <= u < 1/2 */ if (x < 0) { if (Iv_p != NULL) *Iv_p = NPY_NAN; if (Kv_p != NULL) *Kv_p = NPY_NAN; sf_error("ikv_temme", SF_ERROR_DOMAIN, NULL); return; } if (x == 0) { Iv = (v == 0) ? 1 : 0; if (kind & need_k) { sf_error("ikv_temme", SF_ERROR_OVERFLOW, NULL); Kv = NPY_INFINITY; } else { Kv = NPY_NAN; /* any value will do */ } if (reflect && (kind & need_i)) { double z = (u + n % 2); Iv = sin((double)NPY_PI * z) == 0 ? Iv : NPY_INFINITY; if (Iv == NPY_INFINITY || Iv == -NPY_INFINITY) { sf_error("ikv_temme", SF_ERROR_OVERFLOW, NULL); } } if (Iv_p != NULL) { *Iv_p = Iv; } if (Kv_p != NULL) { *Kv_p = Kv; } return; } /* x is positive until reflection */ W = 1 / x; /* Wronskian */ if (x <= 2) { /* x in (0, 2] */ temme_ik_series(u, x, &Ku, &Ku1); /* Temme series */ } else { /* x in (2, \infty) */ CF2_ik(u, x, &Ku, &Ku1); /* continued fraction CF2_ik */ } prev = Ku; current = Ku1; for (k = 1; k <= n; k++) { /* forward recurrence for K */ next = 2 * (u + k) * current / x + prev; prev = current; current = next; } Kv = prev; Kv1 = current; if (kind & need_i) { double lim = (4 * v * v + 10) / (8 * x); lim *= lim; lim *= lim; lim /= 24; if ((lim < MACHEP * 10) && (x > 100)) { /* * x is huge compared to v, CF1 may be very slow * to converge so use asymptotic expansion for large * x case instead. Note that the asymptotic expansion * isn't very accurate - so it's deliberately very hard * to get here - probably we're going to overflow: */ Iv = iv_asymptotic(v, x); } else { CF1_ik(v, x, &fv); /* continued fraction CF1_ik */ Iv = W / (Kv * fv + Kv1); /* Wronskian relation */ } } else { Iv = NPY_NAN; /* any value will do */ } if (reflect) { double z = (u + n % 2); if (Iv_p != NULL) { *Iv_p = Iv + (2 / NPY_PI) * sin(NPY_PI * z) * Kv; /* reflection formula */ } if (Kv_p != NULL) { *Kv_p = Kv; } } else { if (Iv_p != NULL) { *Iv_p = Iv; } if (Kv_p != NULL) { *Kv_p = Kv; } } return; }