src/nmath/bessel_k.c - R - Git at Google

 /*
  *  Mathlib : A C Library of Special Functions
  *  Copyright (C) 1998-2014 Ross Ihaka and the R Core team.
  *  Copyright (C) 2002-3    The R Foundation
  *
  *  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/rkbesl	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))
 #define max0(x, y) (((x) <= (y)) ? (y) : (x))

 static void K_bessel(double *x, double *alpha, int *nb,
 		     int *ize, double *bk, int *ncalc);

 double bessel_k(double x, double alpha, double expo)
 {
     int nb, ncalc, ize;
     double *bk;
 #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_ERROR(ME_RANGE, "bessel_k");
 	return ML_NAN;
     }
     ize = (int)expo;
     if(alpha < 0)
 	alpha = -alpha;
     nb = 1+ (int)floor(alpha);/* nb-1 <= |alpha| < nb */
     alpha -= (double)(nb-1);
 #ifdef MATHLIB_STANDALONE
     bk = (double *) calloc(nb, sizeof(double));
     if (!bk) MATHLIB_ERROR("%s", _("bessel_k allocation error"));
 #else
     vmax = vmaxget();
     bk = (double *) R_alloc((size_t) nb, sizeof(double));
 #endif
     K_bessel(&x, &alpha, &nb, &ize, bk, &ncalc);
     if(ncalc != nb) {/* error input */
       if(ncalc < 0)
 	MATHLIB_WARNING4(_("bessel_k(%g): ncalc (=%d) != nb (=%d); alpha=%g. Arg. out of range?\n"),
 			 x, ncalc, nb, alpha);
       else
 	MATHLIB_WARNING2(_("bessel_k(%g,nu=%g): precision lost in result\n"),
 			 x, alpha+(double)nb-1);
     }
     x = bk[nb-1];
 #ifdef MATHLIB_STANDALONE
     free(bk);
 #else
     vmaxset(vmax);
 #endif
     return x;
 }

 /* modified version of bessel_k that accepts a work array instead of
    allocating one. */
 double bessel_k_ex(double x, double alpha, double expo, double *bk)
 {
     int nb, ncalc, ize;

 #ifdef IEEE_754
     /* NaNs propagated correctly */
     if (ISNAN(x) || ISNAN(alpha)) return x + alpha;
 #endif
     if (x < 0) {
 	ML_ERROR(ME_RANGE, "bessel_k");
 	return ML_NAN;
     }
     ize = (int)expo;
     if(alpha < 0)
 	alpha = -alpha;
     nb = 1+ (int)floor(alpha);/* nb-1 <= |alpha| < nb */
     alpha -= (double)(nb-1);
     K_bessel(&x, &alpha, &nb, &ize, bk, &ncalc);
     if(ncalc != nb) {/* error input */
       if(ncalc < 0)
 	MATHLIB_WARNING4(_("bessel_k(%g): ncalc (=%d) != nb (=%d); alpha=%g. Arg. out of range?\n"),
 			 x, ncalc, nb, alpha);
       else
 	MATHLIB_WARNING2(_("bessel_k(%g,nu=%g): precision lost in result\n"),
 			 x, alpha+(double)nb-1);
     }
     x = bk[nb-1];
     return x;
 }

 static void K_bessel(double *x, double *alpha, int *nb,
 		     int *ize, double *bk, int *ncalc)
 {
 /*-------------------------------------------------------------------

   This routine calculates modified Bessel functions
   of the third kind, K_(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
 	 K's or exponentially scaled K's (K*EXP(X))
 	 are to be calculated.	If K's are to be calculated,
 	 X must not be greater than XMAX_BESS_K.
  ALPHA - Fractional part of order for which
 	 K's or exponentially scaled K's (K*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 K's are to be calculated,
 		    = 2 if exponentially scaled K's are to be calculated.
  BK    - Output vector of length NB.	If the
 	 routine terminates normally (NCALC=NB), the vector BK
 	 contains the functions K(ALPHA,X), ... , K(NB-1+ALPHA,X),
 	 or the corresponding exponentially scaled functions.
 	 If (0 < NCALC < NB), BK(I) contains correct function
 	 values for I <= NCALC, and contains the ratios
 	 K(ALPHA+I-1,X)/K(ALPHA+I-2,X) for the rest of the array.
  NCALC - Output variable indicating possible errors.
 	 Before using the vector BK, 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 K's are
   calculated to the desired accuracy.

   NCALC < -1:  An argument is out of range. For example,
 	NB <= 0, IZE is not 1 or 2, or IZE=1 and ABS(X) >= XMAX_BESS_K.
 	In this case, the B-vector is not calculated,
 	and NCALC is set to MIN0(NB,0)-2	 so that NCALC != NB.
   NCALC = -1:  Either  K(ALPHA,X) >= XINF  or
 	K(ALPHA+NB-1,X)/K(ALPHA+NB-2,X) >= XINF.	 In this case,
 	the B-vector is not calculated.	Note that again
 	NCALC != NB.

   0 < NCALC < NB: Not all requested function values could
 	be calculated accurately.  BK(I) contains correct function
 	values for I <= NCALC, and contains the ratios
 	K(ALPHA+I-1,X)/K(ALPHA+I-2,X) for the rest of the array.


  Intrinsic functions required are:

      ABS, AINT, EXP, INT, LOG, MAX, MIN, SINH, SQRT


  Acknowledgement

 	This program is based on a program written by J. B. Campbell
 	(2) that computes values of the Bessel functions K of float
 	argument and float order.  Modifications include the addition
 	of non-scaled functions, parameterization of machine
 	dependencies, and the use of more accurate approximations
 	for SINH and SIN.

  References: "On Temme's Algorithm for the Modified Bessel
 	      Functions of the Third Kind," Campbell, J. B.,
 	      TOMS 6(4), Dec. 1980, pp. 581-586.

 	     "A FORTRAN IV Subroutine for the Modified Bessel
 	      Functions of the Third Kind of Real Order and Real
 	      Argument," Campbell, J. B., Report NRC/ERB-925,
 	      National Research Council, Canada.

   Latest modification: May 30, 1989

   Modified by: W. J. Cody and L. Stoltz
 	       Applied Mathematics Division
 	       Argonne National Laboratory
 	       Argonne, IL  60439

  -------------------------------------------------------------------
 */
     /*---------------------------------------------------------------------
      * Mathematical constants
      *	A = LOG(2) - Euler's constant
      *	D = SQRT(2/PI)
      ---------------------------------------------------------------------*/
     const static double a = .11593151565841244881;

     /*---------------------------------------------------------------------
       P, Q - Approximation for LOG(GAMMA(1+ALPHA))/ALPHA + Euler's constant
       Coefficients converted from hex to decimal and modified
       by W. J. Cody, 2/26/82 */
     const static double p[8] = { .805629875690432845,20.4045500205365151,
 	    157.705605106676174,536.671116469207504,900.382759291288778,
 	    730.923886650660393,229.299301509425145,.822467033424113231 };
     const static double q[7] = { 29.4601986247850434,277.577868510221208,
 	    1206.70325591027438,2762.91444159791519,3443.74050506564618,
 	    2210.63190113378647,572.267338359892221 };
     /* R, S - Approximation for (1-ALPHA*PI/SIN(ALPHA*PI))/(2.D0*ALPHA) */
     const static double r[5] = { -.48672575865218401848,13.079485869097804016,
 	    -101.96490580880537526,347.65409106507813131,
 	    3.495898124521934782e-4 };
     const static double s[4] = { -25.579105509976461286,212.57260432226544008,
 	    -610.69018684944109624,422.69668805777760407 };
     /* T    - Approximation for SINH(Y)/Y */
     const static double t[6] = { 1.6125990452916363814e-10,
 	    2.5051878502858255354e-8,2.7557319615147964774e-6,
 	    1.9841269840928373686e-4,.0083333333333334751799,
 	    .16666666666666666446 };
     /*---------------------------------------------------------------------*/
     const static double estm[6] = { 52.0583,5.7607,2.7782,14.4303,185.3004, 9.3715 };
     const static double estf[7] = { 41.8341,7.1075,6.4306,42.511,1.35633,84.5096,20.};

     /* Local variables */
     int iend, i, j, k, m, ii, mplus1;
     double x2by4, twox, c, blpha, ratio, wminf;
     double d1, d2, d3, f0, f1, f2, p0, q0, t1, t2, twonu;
     double dm, ex, bk1, bk2, nu;

     ii = 0; /* -Wall */

     ex = *x;
     nu = *alpha;
     *ncalc = min0(*nb,0) - 2;
     if (*nb > 0 && (0. <= nu && nu < 1.) && (1 <= *ize && *ize <= 2)) {
 	if(ex <= 0 || (*ize == 1 && ex > xmax_BESS_K)) {
 	    if(ex <= 0) {
 		if(ex < 0) ML_ERROR(ME_RANGE, "K_bessel");
 		for(i=0; i < *nb; i++)
 		    bk[i] = ML_POSINF;
 	    } else /* would only have underflow */
 		for(i=0; i < *nb; i++)
 		    bk[i] = 0.;
 	    *ncalc = *nb;
 	    return;
 	}
 	k = 0;
 	if (nu < sqxmin_BESS_K) {
 	    nu = 0.;
 	} else if (nu > .5) {
 	    k = 1;
 	    nu -= 1.;
 	}
 	twonu = nu + nu;
 	iend = *nb + k - 1;
 	c = nu * nu;
 	d3 = -c;
 	if (ex <= 1.) {
 	    /* ------------------------------------------------------------
 	       Calculation of P0 = GAMMA(1+ALPHA) * (2/X)**ALPHA
 			      Q0 = GAMMA(1-ALPHA) * (X/2)**ALPHA
 	       ------------------------------------------------------------ */
 	    d1 = 0.; d2 = p[0];
 	    t1 = 1.; t2 = q[0];
 	    for (i = 2; i <= 7; i += 2) {
 		d1 = c * d1 + p[i - 1];
 		d2 = c * d2 + p[i];
 		t1 = c * t1 + q[i - 1];
 		t2 = c * t2 + q[i];
 	    }
 	    d1 = nu * d1;
 	    t1 = nu * t1;
 	    f1 = log(ex);
 	    f0 = a + nu * (p[7] - nu * (d1 + d2) / (t1 + t2)) - f1;
 	    q0 = exp(-nu * (a - nu * (p[7] + nu * (d1-d2) / (t1-t2)) - f1));
 	    f1 = nu * f0;
 	    p0 = exp(f1);
 	    /* -----------------------------------------------------------
 	       Calculation of F0 =
 	       ----------------------------------------------------------- */
 	    d1 = r[4];
 	    t1 = 1.;
 	    for (i = 0; i < 4; ++i) {
 		d1 = c * d1 + r[i];
 		t1 = c * t1 + s[i];
 	    }
 	    /* d2 := sinh(f1)/ nu = sinh(f1)/(f1/f0)
 	     *	   = f0 * sinh(f1)/f1 */
 	    if (fabs(f1) <= .5) {
 		f1 *= f1;
 		d2 = 0.;
 		for (i = 0; i < 6; ++i) {
 		    d2 = f1 * d2 + t[i];
 		}
 		d2 = f0 + f0 * f1 * d2;
 	    } else {
 		d2 = sinh(f1) / nu;
 	    }
 	    f0 = d2 - nu * d1 / (t1 * p0);
 	    if (ex <= 1e-10) {
 		/* ---------------------------------------------------------
 		   X <= 1.0E-10
 		   Calculation of K(ALPHA,X) and X*K(ALPHA+1,X)/K(ALPHA,X)
 		   --------------------------------------------------------- */
 		bk[0] = f0 + ex * f0;
 		if (*ize == 1) {
 		    bk[0] -= ex * bk[0];
 		}
 		ratio = p0 / f0;
 		c = ex * DBL_MAX;
 		if (k != 0) {
 		    /* ---------------------------------------------------
 		       Calculation of K(ALPHA,X)
 		       and  X*K(ALPHA+1,X)/K(ALPHA,X),	ALPHA >= 1/2
 		       --------------------------------------------------- */
 		    *ncalc = -1;
 		    if (bk[0] >= c / ratio) {
 			return;
 		    }
 		    bk[0] = ratio * bk[0] / ex;
 		    twonu += 2.;
 		    ratio = twonu;
 		}
 		*ncalc = 1;
 		if (*nb == 1)
 		    return;

 		/* -----------------------------------------------------
 		   Calculate  K(ALPHA+L,X)/K(ALPHA+L-1,X),
 		   L = 1, 2, ... , NB-1
 		   ----------------------------------------------------- */
 		*ncalc = -1;
 		for (i = 1; i < *nb; ++i) {
 		    if (ratio >= c)
 			return;

 		    bk[i] = ratio / ex;
 		    twonu += 2.;
 		    ratio = twonu;
 		}
 		*ncalc = 1;
 		goto L420;
 	    } else {
 		/* ------------------------------------------------------
 		   10^-10 < X <= 1.0
 		   ------------------------------------------------------ */
 		c = 1.;
 		x2by4 = ex * ex / 4.;
 		p0 = .5 * p0;
 		q0 = .5 * q0;
 		d1 = -1.;
 		d2 = 0.;
 		bk1 = 0.;
 		bk2 = 0.;
 		f1 = f0;
 		f2 = p0;
 		do {
 		    d1 += 2.;
 		    d2 += 1.;
 		    d3 = d1 + d3;
 		    c = x2by4 * c / d2;
 		    f0 = (d2 * f0 + p0 + q0) / d3;
 		    p0 /= d2 - nu;
 		    q0 /= d2 + nu;
 		    t1 = c * f0;
 		    t2 = c * (p0 - d2 * f0);
 		    bk1 += t1;
 		    bk2 += t2;
 		} while (fabs(t1 / (f1 + bk1)) > DBL_EPSILON ||
 			 fabs(t2 / (f2 + bk2)) > DBL_EPSILON);
 		bk1 = f1 + bk1;
 		bk2 = 2. * (f2 + bk2) / ex;
 		if (*ize == 2) {
 		    d1 = exp(ex);
 		    bk1 *= d1;
 		    bk2 *= d1;
 		}
 		wminf = estf[0] * ex + estf[1];
 	    }
 	} else if (DBL_EPSILON * ex > 1.) {
 	    /* -------------------------------------------------
 	       X > 1./EPS
 	       ------------------------------------------------- */
 	    *ncalc = *nb;
 	    bk1 = 1. / (M_SQRT_2dPI * sqrt(ex));
 	    for (i = 0; i < *nb; ++i)
 		bk[i] = bk1;
 	    return;

 	} else {
 	    /* -------------------------------------------------------
 	       X > 1.0
 	       ------------------------------------------------------- */
 	    twox = ex + ex;
 	    blpha = 0.;
 	    ratio = 0.;
 	    if (ex <= 4.) {
 		/* ----------------------------------------------------------
 		   Calculation of K(ALPHA+1,X)/K(ALPHA,X),  1.0 <= X <= 4.0
 		   ----------------------------------------------------------*/
 		d2 = trunc(estm[0] / ex + estm[1]);
 		m = (int) d2;
 		d1 = d2 + d2;
 		d2 -= .5;
 		d2 *= d2;
 		for (i = 2; i <= m; ++i) {
 		    d1 -= 2.;
 		    d2 -= d1;
 		    ratio = (d3 + d2) / (twox + d1 - ratio);
 		}
 		/* -----------------------------------------------------------
 		   Calculation of I(|ALPHA|,X) and I(|ALPHA|+1,X) by backward
 		   recurrence and K(ALPHA,X) from the wronskian
 		   -----------------------------------------------------------*/
 		d2 = trunc(estm[2] * ex + estm[3]);
 		m = (int) d2;
 		c = fabs(nu);
 		d3 = c + c;
 		d1 = d3 - 1.;
 		f1 = DBL_MIN;
 		f0 = (2. * (c + d2) / ex + .5 * ex / (c + d2 + 1.)) * DBL_MIN;
 		for (i = 3; i <= m; ++i) {
 		    d2 -= 1.;
 		    f2 = (d3 + d2 + d2) * f0;
 		    blpha = (1. + d1 / d2) * (f2 + blpha);
 		    f2 = f2 / ex + f1;
 		    f1 = f0;
 		    f0 = f2;
 		}
 		f1 = (d3 + 2.) * f0 / ex + f1;
 		d1 = 0.;
 		t1 = 1.;
 		for (i = 1; i <= 7; ++i) {
 		    d1 = c * d1 + p[i - 1];
 		    t1 = c * t1 + q[i - 1];
 		}
 		p0 = exp(c * (a + c * (p[7] - c * d1 / t1) - log(ex))) / ex;
 		f2 = (c + .5 - ratio) * f1 / ex;
 		bk1 = p0 + (d3 * f0 - f2 + f0 + blpha) / (f2 + f1 + f0) * p0;
 		if (*ize == 1) {
 		    bk1 *= exp(-ex);
 		}
 		wminf = estf[2] * ex + estf[3];
 	    } else {
 		/* ---------------------------------------------------------
 		   Calculation of K(ALPHA,X) and K(ALPHA+1,X)/K(ALPHA,X), by
 		   backward recurrence, for  X > 4.0
 		   ----------------------------------------------------------*/
 		dm = trunc(estm[4] / ex + estm[5]);
 		m = (int) dm;
 		d2 = dm - .5;
 		d2 *= d2;
 		d1 = dm + dm;
 		for (i = 2; i <= m; ++i) {
 		    dm -= 1.;
 		    d1 -= 2.;
 		    d2 -= d1;
 		    ratio = (d3 + d2) / (twox + d1 - ratio);
 		    blpha = (ratio + ratio * blpha) / dm;
 		}
 		bk1 = 1. / ((M_SQRT_2dPI + M_SQRT_2dPI * blpha) * sqrt(ex));
 		if (*ize == 1)
 		    bk1 *= exp(-ex);
 		wminf = estf[4] * (ex - fabs(ex - estf[6])) + estf[5];
 	    }
 	    /* ---------------------------------------------------------
 	       Calculation of K(ALPHA+1,X)
 	       from K(ALPHA,X) and  K(ALPHA+1,X)/K(ALPHA,X)
 	       --------------------------------------------------------- */
 	    bk2 = bk1 + bk1 * (nu + .5 - ratio) / ex;
 	}
 	/*--------------------------------------------------------------------
 	  Calculation of 'NCALC', K(ALPHA+I,X),	I  =  0, 1, ... , NCALC-1,
 	  &	  K(ALPHA+I,X)/K(ALPHA+I-1,X),	I = NCALC, NCALC+1, ... , NB-1
 	  -------------------------------------------------------------------*/
 	*ncalc = *nb;
 	bk[0] = bk1;
 	if (iend == 0)
 	    return;

 	j = 1 - k;
 	if (j >= 0)
 	    bk[j] = bk2;

 	if (iend == 1)
 	    return;

 	m = min0((int) (wminf - nu),iend);
 	for (i = 2; i <= m; ++i) {
 	    t1 = bk1;
 	    bk1 = bk2;
 	    twonu += 2.;
 	    if (ex < 1.) {
 		if (bk1 >= DBL_MAX / twonu * ex)
 		    break;
 	    } else {
 		if (bk1 / ex >= DBL_MAX / twonu)
 		    break;
 	    }
 	    bk2 = twonu / ex * bk1 + t1;
 	    ii = i;
 	    ++j;
 	    if (j >= 0) {
 		bk[j] = bk2;
 	    }
 	}

 	m = ii;
 	if (m == iend) {
 	    return;
 	}
 	ratio = bk2 / bk1;
 	mplus1 = m + 1;
 	*ncalc = -1;
 	for (i = mplus1; i <= iend; ++i) {
 	    twonu += 2.;
 	    ratio = twonu / ex + 1./ratio;
 	    ++j;
 	    if (j >= 1) {
 		bk[j] = ratio;
 	    } else {
 		if (bk2 >= DBL_MAX / ratio)
 		    return;

 		bk2 *= ratio;
 	    }
 	}
 	*ncalc = max0(1, mplus1 - k);
 	if (*ncalc == 1)
 	    bk[0] = bk2;
 	if (*nb == 1)
 	    return;

 L420:
 	for (i = *ncalc; i < *nb; ++i) { /* i == *ncalc */
 #ifndef IEEE_754
 	    if (bk[i-1] >= DBL_MAX / bk[i])
 		return;
 #endif
 	    bk[i] *= bk[i-1];
 	    (*ncalc)++;
 	}
     }
 }
	/*
	* Mathlib : A C Library of Special Functions
	* Copyright (C) 1998-2014 Ross Ihaka and the R Core team.
	* Copyright (C) 2002-3 The R Foundation
	*
	* 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/rkbesl 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))
	#define max0(x, y) (((x) <= (y)) ? (y) : (x))

	static void K_bessel(double x, double alpha, int *nb,
	int ize, double bk, int *ncalc);

	double bessel_k(double x, double alpha, double expo)
	{
	int nb, ncalc, ize;
	double *bk;
	#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_ERROR(ME_RANGE, "bessel_k");
	return ML_NAN;
	}
	ize = (int)expo;
	if(alpha < 0)
	alpha = -alpha;
	nb = 1+ (int)floor(alpha);/* nb-1 <= \|alpha\| < nb */
	alpha -= (double)(nb-1);
	#ifdef MATHLIB_STANDALONE
	bk = (double *) calloc(nb, sizeof(double));
	if (!bk) MATHLIB_ERROR("%s", _("bessel_k allocation error"));
	#else
	vmax = vmaxget();
	bk = (double *) R_alloc((size_t) nb, sizeof(double));
	#endif
	K_bessel(&x, &alpha, &nb, &ize, bk, &ncalc);
	if(ncalc != nb) {/* error input */
	if(ncalc < 0)
	MATHLIB_WARNING4(_("bessel_k(%g): ncalc (=%d) != nb (=%d); alpha=%g. Arg. out of range?\n"),
	x, ncalc, nb, alpha);
	else
	MATHLIB_WARNING2(_("bessel_k(%g,nu=%g): precision lost in result\n"),
	x, alpha+(double)nb-1);
	}
	x = bk[nb-1];
	#ifdef MATHLIB_STANDALONE
	free(bk);
	#else
	vmaxset(vmax);
	#endif
	return x;
	}

	/* modified version of bessel_k that accepts a work array instead of
	allocating one. */
	double bessel_k_ex(double x, double alpha, double expo, double *bk)
	{
	int nb, ncalc, ize;

	#ifdef IEEE_754
	/* NaNs propagated correctly */
	if (ISNAN(x) \|\| ISNAN(alpha)) return x + alpha;
	#endif
	if (x < 0) {
	ML_ERROR(ME_RANGE, "bessel_k");
	return ML_NAN;
	}
	ize = (int)expo;
	if(alpha < 0)
	alpha = -alpha;
	nb = 1+ (int)floor(alpha);/* nb-1 <= \|alpha\| < nb */
	alpha -= (double)(nb-1);
	K_bessel(&x, &alpha, &nb, &ize, bk, &ncalc);
	if(ncalc != nb) {/* error input */
	if(ncalc < 0)
	MATHLIB_WARNING4(_("bessel_k(%g): ncalc (=%d) != nb (=%d); alpha=%g. Arg. out of range?\n"),
	x, ncalc, nb, alpha);
	else
	MATHLIB_WARNING2(_("bessel_k(%g,nu=%g): precision lost in result\n"),
	x, alpha+(double)nb-1);
	}
	x = bk[nb-1];
	return x;
	}

	static void K_bessel(double x, double alpha, int *nb,
	int ize, double bk, int *ncalc)
	{
	/*-------------------------------------------------------------------

	This routine calculates modified Bessel functions
	of the third kind, K_(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
	K's or exponentially scaled K's (K*EXP(X))
	are to be calculated. If K's are to be calculated,
	X must not be greater than XMAX_BESS_K.
	ALPHA - Fractional part of order for which
	K's or exponentially scaled K's (K*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 K's are to be calculated,
	= 2 if exponentially scaled K's are to be calculated.
	BK - Output vector of length NB. If the
	routine terminates normally (NCALC=NB), the vector BK
	contains the functions K(ALPHA,X), ... , K(NB-1+ALPHA,X),
	or the corresponding exponentially scaled functions.
	If (0 < NCALC < NB), BK(I) contains correct function
	values for I <= NCALC, and contains the ratios
	K(ALPHA+I-1,X)/K(ALPHA+I-2,X) for the rest of the array.
	NCALC - Output variable indicating possible errors.
	Before using the vector BK, 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 K's are
	calculated to the desired accuracy.

	NCALC < -1: An argument is out of range. For example,
	NB <= 0, IZE is not 1 or 2, or IZE=1 and ABS(X) >= XMAX_BESS_K.
	In this case, the B-vector is not calculated,
	and NCALC is set to MIN0(NB,0)-2 so that NCALC != NB.
	NCALC = -1: Either K(ALPHA,X) >= XINF or
	K(ALPHA+NB-1,X)/K(ALPHA+NB-2,X) >= XINF. In this case,
	the B-vector is not calculated. Note that again
	NCALC != NB.

	0 < NCALC < NB: Not all requested function values could
	be calculated accurately. BK(I) contains correct function
	values for I <= NCALC, and contains the ratios
	K(ALPHA+I-1,X)/K(ALPHA+I-2,X) for the rest of the array.


	Intrinsic functions required are:

	ABS, AINT, EXP, INT, LOG, MAX, MIN, SINH, SQRT


	Acknowledgement

	This program is based on a program written by J. B. Campbell
	(2) that computes values of the Bessel functions K of float
	argument and float order. Modifications include the addition
	of non-scaled functions, parameterization of machine
	dependencies, and the use of more accurate approximations
	for SINH and SIN.

	References: "On Temme's Algorithm for the Modified Bessel
	Functions of the Third Kind," Campbell, J. B.,
	TOMS 6(4), Dec. 1980, pp. 581-586.

	"A FORTRAN IV Subroutine for the Modified Bessel
	Functions of the Third Kind of Real Order and Real
	Argument," Campbell, J. B., Report NRC/ERB-925,
	National Research Council, Canada.

	Latest modification: May 30, 1989

	Modified by: W. J. Cody and L. Stoltz
	Applied Mathematics Division
	Argonne National Laboratory
	Argonne, IL 60439

	-------------------------------------------------------------------
	*/
	/*---------------------------------------------------------------------
	* Mathematical constants
	* A = LOG(2) - Euler's constant
	* D = SQRT(2/PI)
	---------------------------------------------------------------------*/
	const static double a = .11593151565841244881;

	/*---------------------------------------------------------------------
	P, Q - Approximation for LOG(GAMMA(1+ALPHA))/ALPHA + Euler's constant
	Coefficients converted from hex to decimal and modified
	by W. J. Cody, 2/26/82 */
	const static double p[8] = { .805629875690432845,20.4045500205365151,
	157.705605106676174,536.671116469207504,900.382759291288778,
	730.923886650660393,229.299301509425145,.822467033424113231 };
	const static double q[7] = { 29.4601986247850434,277.577868510221208,
	1206.70325591027438,2762.91444159791519,3443.74050506564618,
	2210.63190113378647,572.267338359892221 };
	/* R, S - Approximation for (1-ALPHAPI/SIN(ALPHAPI))/(2.D0ALPHA) /
	const static double r[5] = { -.48672575865218401848,13.079485869097804016,
	-101.96490580880537526,347.65409106507813131,
	3.495898124521934782e-4 };
	const static double s[4] = { -25.579105509976461286,212.57260432226544008,
	-610.69018684944109624,422.69668805777760407 };
	/* T - Approximation for SINH(Y)/Y */
	const static double t[6] = { 1.6125990452916363814e-10,
	2.5051878502858255354e-8,2.7557319615147964774e-6,
	1.9841269840928373686e-4,.0083333333333334751799,
	.16666666666666666446 };
	/---------------------------------------------------------------------/
	const static double estm[6] = { 52.0583,5.7607,2.7782,14.4303,185.3004, 9.3715 };
	const static double estf[7] = { 41.8341,7.1075,6.4306,42.511,1.35633,84.5096,20.};

	/* Local variables */
	int iend, i, j, k, m, ii, mplus1;
	double x2by4, twox, c, blpha, ratio, wminf;
	double d1, d2, d3, f0, f1, f2, p0, q0, t1, t2, twonu;
	double dm, ex, bk1, bk2, nu;

	ii = 0; /* -Wall */

	ex = *x;
	nu = *alpha;
	ncalc = min0(nb,0) - 2;
	if (nb > 0 && (0. <= nu && nu < 1.) && (1 <= ize && *ize <= 2)) {
	if(ex <= 0 \|\| (*ize == 1 && ex > xmax_BESS_K)) {
	if(ex <= 0) {
	if(ex < 0) ML_ERROR(ME_RANGE, "K_bessel");
	for(i=0; i < *nb; i++)
	bk[i] = ML_POSINF;
	} else /* would only have underflow */
	for(i=0; i < *nb; i++)
	bk[i] = 0.;
	ncalc = nb;
	return;
	}
	k = 0;
	if (nu < sqxmin_BESS_K) {
	nu = 0.;
	} else if (nu > .5) {
	k = 1;
	nu -= 1.;
	}
	twonu = nu + nu;
	iend = *nb + k - 1;
	c = nu * nu;
	d3 = -c;
	if (ex <= 1.) {
	/* ------------------------------------------------------------
	Calculation of P0 = GAMMA(1+ALPHA) * (2/X)**ALPHA
	Q0 = GAMMA(1-ALPHA) * (X/2)**ALPHA
	------------------------------------------------------------ */
	d1 = 0.; d2 = p[0];
	t1 = 1.; t2 = q[0];
	for (i = 2; i <= 7; i += 2) {
	d1 = c * d1 + p[i - 1];
	d2 = c * d2 + p[i];
	t1 = c * t1 + q[i - 1];
	t2 = c * t2 + q[i];
	}
	d1 = nu * d1;
	t1 = nu * t1;
	f1 = log(ex);
	f0 = a + nu * (p[7] - nu * (d1 + d2) / (t1 + t2)) - f1;
	q0 = exp(-nu * (a - nu * (p[7] + nu * (d1-d2) / (t1-t2)) - f1));
	f1 = nu * f0;
	p0 = exp(f1);
	/* -----------------------------------------------------------
	Calculation of F0 =
	----------------------------------------------------------- */
	d1 = r[4];
	t1 = 1.;
	for (i = 0; i < 4; ++i) {
	d1 = c * d1 + r[i];
	t1 = c * t1 + s[i];
	}
	/* d2 := sinh(f1)/ nu = sinh(f1)/(f1/f0)
	* = f0 * sinh(f1)/f1 */
	if (fabs(f1) <= .5) {
	f1 *= f1;
	d2 = 0.;
	for (i = 0; i < 6; ++i) {
	d2 = f1 * d2 + t[i];
	}
	d2 = f0 + f0 * f1 * d2;
	} else {
	d2 = sinh(f1) / nu;
	}
	f0 = d2 - nu * d1 / (t1 * p0);
	if (ex <= 1e-10) {
	/* ---------------------------------------------------------
	X <= 1.0E-10
	Calculation of K(ALPHA,X) and X*K(ALPHA+1,X)/K(ALPHA,X)
	--------------------------------------------------------- */
	bk[0] = f0 + ex * f0;
	if (*ize == 1) {
	bk[0] -= ex * bk[0];
	}
	ratio = p0 / f0;
	c = ex * DBL_MAX;
	if (k != 0) {
	/* ---------------------------------------------------
	Calculation of K(ALPHA,X)
	and X*K(ALPHA+1,X)/K(ALPHA,X), ALPHA >= 1/2
	--------------------------------------------------- */
	*ncalc = -1;
	if (bk[0] >= c / ratio) {
	return;
	}
	bk[0] = ratio * bk[0] / ex;
	twonu += 2.;
	ratio = twonu;
	}
	*ncalc = 1;
	if (*nb == 1)
	return;

	/* -----------------------------------------------------
	Calculate K(ALPHA+L,X)/K(ALPHA+L-1,X),
	L = 1, 2, ... , NB-1
	----------------------------------------------------- */
	*ncalc = -1;
	for (i = 1; i < *nb; ++i) {
	if (ratio >= c)
	return;

	bk[i] = ratio / ex;
	twonu += 2.;
	ratio = twonu;
	}
	*ncalc = 1;
	goto L420;
	} else {
	/* ------------------------------------------------------
	10^-10 < X <= 1.0
	------------------------------------------------------ */
	c = 1.;
	x2by4 = ex * ex / 4.;
	p0 = .5 * p0;
	q0 = .5 * q0;
	d1 = -1.;
	d2 = 0.;
	bk1 = 0.;
	bk2 = 0.;
	f1 = f0;
	f2 = p0;
	do {
	d1 += 2.;
	d2 += 1.;
	d3 = d1 + d3;
	c = x2by4 * c / d2;
	f0 = (d2 * f0 + p0 + q0) / d3;
	p0 /= d2 - nu;
	q0 /= d2 + nu;
	t1 = c * f0;
	t2 = c * (p0 - d2 * f0);
	bk1 += t1;
	bk2 += t2;
	} while (fabs(t1 / (f1 + bk1)) > DBL_EPSILON \|\|
	fabs(t2 / (f2 + bk2)) > DBL_EPSILON);
	bk1 = f1 + bk1;
	bk2 = 2. * (f2 + bk2) / ex;
	if (*ize == 2) {
	d1 = exp(ex);
	bk1 *= d1;
	bk2 *= d1;
	}
	wminf = estf[0] * ex + estf[1];
	}
	} else if (DBL_EPSILON * ex > 1.) {
	/* -------------------------------------------------
	X > 1./EPS
	------------------------------------------------- */
	ncalc = nb;
	bk1 = 1. / (M_SQRT_2dPI * sqrt(ex));
	for (i = 0; i < *nb; ++i)
	bk[i] = bk1;
	return;

	} else {
	/* -------------------------------------------------------
	X > 1.0
	------------------------------------------------------- */
	twox = ex + ex;
	blpha = 0.;
	ratio = 0.;
	if (ex <= 4.) {
	/* ----------------------------------------------------------
	Calculation of K(ALPHA+1,X)/K(ALPHA,X), 1.0 <= X <= 4.0
	----------------------------------------------------------*/
	d2 = trunc(estm[0] / ex + estm[1]);
	m = (int) d2;
	d1 = d2 + d2;
	d2 -= .5;
	d2 *= d2;
	for (i = 2; i <= m; ++i) {
	d1 -= 2.;
	d2 -= d1;
	ratio = (d3 + d2) / (twox + d1 - ratio);
	}
	/* -----------------------------------------------------------
	Calculation of I(\|ALPHA\|,X) and I(\|ALPHA\|+1,X) by backward
	recurrence and K(ALPHA,X) from the wronskian
	-----------------------------------------------------------*/
	d2 = trunc(estm[2] * ex + estm[3]);
	m = (int) d2;
	c = fabs(nu);
	d3 = c + c;
	d1 = d3 - 1.;
	f1 = DBL_MIN;
	f0 = (2. * (c + d2) / ex + .5 * ex / (c + d2 + 1.)) * DBL_MIN;
	for (i = 3; i <= m; ++i) {
	d2 -= 1.;
	f2 = (d3 + d2 + d2) * f0;
	blpha = (1. + d1 / d2) * (f2 + blpha);
	f2 = f2 / ex + f1;
	f1 = f0;
	f0 = f2;
	}
	f1 = (d3 + 2.) * f0 / ex + f1;
	d1 = 0.;
	t1 = 1.;
	for (i = 1; i <= 7; ++i) {
	d1 = c * d1 + p[i - 1];
	t1 = c * t1 + q[i - 1];
	}
	p0 = exp(c * (a + c * (p[7] - c * d1 / t1) - log(ex))) / ex;
	f2 = (c + .5 - ratio) * f1 / ex;
	bk1 = p0 + (d3 * f0 - f2 + f0 + blpha) / (f2 + f1 + f0) * p0;
	if (*ize == 1) {
	bk1 *= exp(-ex);
	}
	wminf = estf[2] * ex + estf[3];
	} else {
	/* ---------------------------------------------------------
	Calculation of K(ALPHA,X) and K(ALPHA+1,X)/K(ALPHA,X), by
	backward recurrence, for X > 4.0
	----------------------------------------------------------*/
	dm = trunc(estm[4] / ex + estm[5]);
	m = (int) dm;
	d2 = dm - .5;
	d2 *= d2;
	d1 = dm + dm;
	for (i = 2; i <= m; ++i) {
	dm -= 1.;
	d1 -= 2.;
	d2 -= d1;
	ratio = (d3 + d2) / (twox + d1 - ratio);
	blpha = (ratio + ratio * blpha) / dm;
	}
	bk1 = 1. / ((M_SQRT_2dPI + M_SQRT_2dPI * blpha) * sqrt(ex));
	if (*ize == 1)
	bk1 *= exp(-ex);
	wminf = estf[4] * (ex - fabs(ex - estf[6])) + estf[5];
	}
	/* ---------------------------------------------------------
	Calculation of K(ALPHA+1,X)
	from K(ALPHA,X) and K(ALPHA+1,X)/K(ALPHA,X)
	--------------------------------------------------------- */
	bk2 = bk1 + bk1 * (nu + .5 - ratio) / ex;
	}
	/*--------------------------------------------------------------------
	Calculation of 'NCALC', K(ALPHA+I,X), I = 0, 1, ... , NCALC-1,
	& K(ALPHA+I,X)/K(ALPHA+I-1,X), I = NCALC, NCALC+1, ... , NB-1
	-------------------------------------------------------------------*/
	ncalc = nb;
	bk[0] = bk1;
	if (iend == 0)
	return;

	j = 1 - k;
	if (j >= 0)
	bk[j] = bk2;

	if (iend == 1)
	return;

	m = min0((int) (wminf - nu),iend);
	for (i = 2; i <= m; ++i) {
	t1 = bk1;
	bk1 = bk2;
	twonu += 2.;
	if (ex < 1.) {
	if (bk1 >= DBL_MAX / twonu * ex)
	break;
	} else {
	if (bk1 / ex >= DBL_MAX / twonu)
	break;
	}
	bk2 = twonu / ex * bk1 + t1;
	ii = i;
	++j;
	if (j >= 0) {
	bk[j] = bk2;
	}
	}

	m = ii;
	if (m == iend) {
	return;
	}
	ratio = bk2 / bk1;
	mplus1 = m + 1;
	*ncalc = -1;
	for (i = mplus1; i <= iend; ++i) {
	twonu += 2.;
	ratio = twonu / ex + 1./ratio;
	++j;
	if (j >= 1) {
	bk[j] = ratio;
	} else {
	if (bk2 >= DBL_MAX / ratio)
	return;

	bk2 *= ratio;
	}
	}
	*ncalc = max0(1, mplus1 - k);
	if (*ncalc == 1)
	bk[0] = bk2;
	if (*nb == 1)
	return;

	L420:
	for (i = ncalc; i < nb; ++i) { /* i == ncalc /
	#ifndef IEEE_754
	if (bk[i-1] >= DBL_MAX / bk[i])
	return;
	#endif
	bk[i] *= bk[i-1];
	(*ncalc)++;
	}
	}
	}