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

 /*
  *  Mathlib : A C Library of Special Functions
  *  Copyright (C) 1998-2015 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/rybesl	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 Y_bessel(double *x, double *alpha, int *nb,
 		     double *by, int *ncalc);

 // unused now from R
 double bessel_y(double x, double alpha)
 {
     int nb, ncalc;
     double na, *by;
 #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_y");
 	return ML_NAN;
     }
     na = floor(alpha);
     if (alpha < 0) {
 	/* Using Abramowitz & Stegun  9.1.2
 	 * this may not be quite optimal (CPU and accuracy wise) */
 	return(((alpha - na == 0.5) ? 0 : bessel_y(x, -alpha) * cospi(alpha)) -
 	       ((alpha      == na ) ? 0 : bessel_j(x, -alpha) * sinpi(alpha)));
     }
     else if (alpha > 1e7) {
 	MATHLIB_WARNING(_("besselY(x, nu): nu=%g too large for bessel_y() algorithm"),
 			alpha);
 	return ML_NAN;
     }
     nb = 1+ (int)na;/* nb-1 <= alpha < nb */
     alpha -= (double)(nb-1);
 #ifdef MATHLIB_STANDALONE
     by = (double *) calloc(nb, sizeof(double));
     if (!by) MATHLIB_ERROR("%s", _("bessel_y allocation error"));
 #else
     vmax = vmaxget();
     by = (double *) R_alloc((size_t) nb, sizeof(double));
 #endif
     Y_bessel(&x, &alpha, &nb, by, &ncalc);
     if(ncalc != nb) {/* error input */
 	if(ncalc == -1) {
 #ifdef MATHLIB_STANDALONE
 	    free(by);
 #else
 	    vmaxset(vmax);
 #endif
 	    return ML_POSINF;
 	}
 	else if(ncalc < -1)
 	    MATHLIB_WARNING4(_("bessel_y(%g): ncalc (=%d) != nb (=%d); alpha=%g. Arg. out of range?\n"),
 			     x, ncalc, nb, alpha);
 	else /* ncalc >= 0 */
 	    MATHLIB_WARNING2(_("bessel_y(%g,nu=%g): precision lost in result\n"),
 			     x, alpha+(double)nb-1);
     }
     x = by[nb-1];
 #ifdef MATHLIB_STANDALONE
     free(by);
 #else
     vmaxset(vmax);
 #endif
     return x;
 }

 /* Called from R: modified version of bessel_y(), accepting a work array
  * instead of allocating one. */
 double bessel_y_ex(double x, double alpha, double *by)
 {
     int nb, ncalc;
     double na;

 #ifdef IEEE_754
     /* NaNs propagated correctly */
     if (ISNAN(x) || ISNAN(alpha)) return x + alpha;
 #endif
     if (x < 0) {
 	ML_ERROR(ME_RANGE, "bessel_y");
 	return ML_NAN;
     }
     na = floor(alpha);
     if (alpha < 0) {
 	/* Using Abramowitz & Stegun  9.1.2
 	 * this may not be quite optimal (CPU and accuracy wise) */
 	return(((alpha - na == 0.5) ? 0 : bessel_y_ex(x, -alpha, by) * cospi(alpha)) -
 	       ((alpha      == na ) ? 0 : bessel_j_ex(x, -alpha, by) * sinpi(alpha)));
     }
     else if (alpha > 1e7) {
 	MATHLIB_WARNING(_("besselY(x, nu): nu=%g too large for bessel_y() algorithm"),
 			alpha);
 	return ML_NAN;
     }
     nb = 1+ (int)na;/* nb-1 <= alpha < nb */
     alpha -= (double)(nb-1);
     Y_bessel(&x, &alpha, &nb, by, &ncalc);
     if(ncalc != nb) {/* error input */
 	if(ncalc == -1)
 	    return ML_POSINF;
 	else if(ncalc < -1)
 	    MATHLIB_WARNING4(_("bessel_y(%g): ncalc (=%d) != nb (=%d); alpha=%g. Arg. out of range?\n"),
 			     x, ncalc, nb, alpha);
 	else /* ncalc >= 0 */
 	    MATHLIB_WARNING2(_("bessel_y(%g,nu=%g): precision lost in result\n"),
 			     x, alpha+(double)nb-1);
     }
     x = by[nb-1];
     return x;
 }

 static void Y_bessel(double *x, double *alpha, int *nb,
 		     double *by, int *ncalc)
 {
 /* ----------------------------------------------------------------------

  This routine calculates Bessel functions Y_(N+ALPHA) (X)
 v for non-negative argument X, and non-negative order N+ALPHA.


  Explanation of variables in the calling sequence

  X     - Non-negative argument for which
 	 Y's are to be calculated.
  ALPHA - Fractional part of order for which
 	 Y's 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).
  BY    - Output vector of length NB.	If the
 	 routine terminates normally (NCALC=NB), the vector BY
 	 contains the functions Y(ALPHA,X), ... , Y(NB-1+ALPHA,X),
 	 If (0 < NCALC < NB), BY(I) contains correct function
 	 values for I <= NCALC, and contains the ratios
 	 Y(ALPHA+I-1,X)/Y(ALPHA+I-2,X) for the rest of the array.
  NCALC - Output variable indicating possible errors.
 	 Before using the vector BY, 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 Y'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.  In this case, BY[0] = 0.0, the remainder of the
 	BY-vector is not calculated, and NCALC is set to
 	MIN0(NB,0)-2  so that NCALC != NB.
   NCALC = -1:  Y(ALPHA,X) >= XINF.  The requested function
 	values are set to 0.0.
   1 < NCALC < NB: Not all requested function values could
 	be calculated accurately.  BY(I) contains correct function
 	values for I <= NCALC, and and the remaining NB-NCALC
 	array elements contain 0.0.


  Intrinsic functions required are:

      DBLE, EXP, INT, MAX, MIN, REAL, SQRT


  Acknowledgement

 	This program draws heavily on Temme's Algol program for Y(a,x)
 	and Y(a+1,x) and on Campbell's programs for Y_nu(x).	Temme's
 	scheme is used for  x < THRESH, and Campbell's scheme is used
 	in the asymptotic region.  Segments of code from both sources
 	have been translated into Fortran 77, merged, and heavily modified.
 	Modifications include parameterization of machine dependencies,
 	use of a new approximation for ln(gamma(x)), and built-in
 	protection against over/underflow.

  References: "Bessel functions J_nu(x) and Y_nu(x) of float
 	      order and float argument," Campbell, J. B.,
 	      Comp. Phy. Comm. 18, 1979, pp. 133-142.

 	     "On the numerical evaluation of the ordinary
 	      Bessel function of the second kind," Temme,
 	      N. M., J. Comput. Phys. 21, 1976, pp. 343-350.

   Latest modification: March 19, 1990

   Modified by: W. J. Cody
 	       Applied Mathematics Division
 	       Argonne National Laboratory
 	       Argonne, IL  60439
  ----------------------------------------------------------------------*/


 /* ----------------------------------------------------------------------
   Mathematical constants
     FIVPI = 5*PI
     PIM5 = 5*PI - 15
  ----------------------------------------------------------------------*/
     const static double fivpi = 15.707963267948966192;
     const static double pim5	=   .70796326794896619231;

     /*----------------------------------------------------------------------
       Coefficients for Chebyshev polynomial expansion of
       1/gamma(1-x), abs(x) <= .5
       ----------------------------------------------------------------------*/
     const static double ch[21] = { -6.7735241822398840964e-24,
 	    -6.1455180116049879894e-23,2.9017595056104745456e-21,
 	    1.3639417919073099464e-19,2.3826220476859635824e-18,
 	    -9.0642907957550702534e-18,-1.4943667065169001769e-15,
 	    -3.3919078305362211264e-14,-1.7023776642512729175e-13,
 	    9.1609750938768647911e-12,2.4230957900482704055e-10,
 	    1.7451364971382984243e-9,-3.3126119768180852711e-8,
 	    -8.6592079961391259661e-7,-4.9717367041957398581e-6,
 	    7.6309597585908126618e-5,.0012719271366545622927,
 	    .0017063050710955562222,-.07685284084478667369,
 	    -.28387654227602353814,.92187029365045265648 };

     /* Local variables */
     int i, k, na;

     double alfa, div, ddiv, even, gamma, term, cosmu, sinmu,
 	b, c, d, e, f, g, h, p, q, r, s, d1, d2, q0, pa,pa1, qa,qa1,
 	en, en1, nu, ex,  ya,ya1, twobyx, den, odd, aye, dmu, x2, xna;

     en1 = ya = ya1 = 0;		/* -Wall */

     ex = *x;
     nu = *alpha;
     if (*nb > 0 && 0. <= nu && nu < 1.) {
 	if(ex < DBL_MIN || ex > xlrg_BESS_Y) {
 	    /* Warning is not really appropriate, give
 	     * proper limit:
 	     * ML_ERROR(ME_RANGE, "Y_bessel"); */
 	    *ncalc = *nb;
 	    if(ex > xlrg_BESS_Y)  by[0]= 0.; /*was ML_POSINF */
 	    else if(ex < DBL_MIN) by[0]=ML_NEGINF;
 	    for(i=0; i < *nb; i++)
 		by[i] = by[0];
 	    return;
 	}
 	xna = trunc(nu + .5);
 	na = (int) xna;
 	if (na == 1) {/* <==>  .5 <= *alpha < 1	 <==>  -5. <= nu < 0 */
 	    nu -= xna;
 	}
 	if (nu == -.5) {
 	    p = M_SQRT_2dPI / sqrt(ex);
 	    ya = p * sin(ex);
 	    ya1 = -p * cos(ex);
 	} else if (ex < 3.) {
 	    /* -------------------------------------------------------------
 	       Use Temme's scheme for small X
 	       ------------------------------------------------------------- */
 	    b = ex * .5;
 	    d = -log(b);
 	    f = nu * d;
 	    e = pow(b, -nu);
 	    if (fabs(nu) < M_eps_sinc)
 		c = M_1_PI;
 	    else
 		c = nu / sinpi(nu);

 	    /* ------------------------------------------------------------
 	       Computation of sinh(f)/f
 	       ------------------------------------------------------------ */
 	    if (fabs(f) < 1.) {
 		x2 = f * f;
 		en = 19.;
 		s = 1.;
 		for (i = 1; i <= 9; ++i) {
 		    s = s * x2 / en / (en - 1.) + 1.;
 		    en -= 2.;
 		}
 	    } else {
 		s = (e - 1. / e) * .5 / f;
 	    }
 	    /* --------------------------------------------------------
 	       Computation of 1/gamma(1-a) using Chebyshev polynomials */
 	    x2 = nu * nu * 8.;
 	    aye = ch[0];
 	    even = 0.;
 	    alfa = ch[1];
 	    odd = 0.;
 	    for (i = 3; i <= 19; i += 2) {
 		even = -(aye + aye + even);
 		aye = -even * x2 - aye + ch[i - 1];
 		odd = -(alfa + alfa + odd);
 		alfa = -odd * x2 - alfa + ch[i];
 	    }
 	    even = (even * .5 + aye) * x2 - aye + ch[20];
 	    odd = (odd + alfa) * 2.;
 	    gamma = odd * nu + even;
 	    /* End of computation of 1/gamma(1-a)
 	       ----------------------------------------------------------- */
 	    g = e * gamma;
 	    e = (e + 1. / e) * .5;
 	    f = 2. * c * (odd * e + even * s * d);
 	    e = nu * nu;
 	    p = g * c;
 	    q = M_1_PI / g;
 	    c = nu * M_PI_2;
 	    if (fabs(c) < M_eps_sinc)
 		r = 1.;
 	    else
 		r = sinpi(nu/2) / c;

 	    r = M_PI * c * r * r;
 	    c = 1.;
 	    d = -b * b;
 	    h = 0.;
 	    ya = f + r * q;
 	    ya1 = p;
 	    en = 1.;

 	    while (fabs(g / (1. + fabs(ya))) +
 		   fabs(h / (1. + fabs(ya1))) > DBL_EPSILON) {
 		f = (f * en + p + q) / (en * en - e);
 		c *= (d / en);
 		p /= en - nu;
 		q /= en + nu;
 		g = c * (f + r * q);
 		h = c * p - en * g;
 		ya += g;
 		ya1+= h;
 		en += 1.;
 	    }
 	    ya = -ya;
 	    ya1 = -ya1 / b;
 	} else if (ex < thresh_BESS_Y) {
 	    /* --------------------------------------------------------------
 	       Use Temme's scheme for moderate X :  3 <= x < 16
 	       -------------------------------------------------------------- */
 	    c = (.5 - nu) * (.5 + nu);
 	    b = ex + ex;
 	    e = ex * M_1_PI * cospi(nu) / DBL_EPSILON;
 	    e *= e;
 	    p = 1.;
 	    q = -ex;
 	    r = 1. + ex * ex;
 	    s = r;
 	    en = 2.;
 	    while (r * en * en < e) {
 		en1 = en + 1.;
 		d = (en - 1. + c / en) / s;
 		p = (en + en - p * d) / en1;
 		q = (-b + q * d) / en1;
 		s = p * p + q * q;
 		r *= s;
 		en = en1;
 	    }
 	    f = p / s;
 	    p = f;
 	    g = -q / s;
 	    q = g;
 L220:
 	    en -= 1.;
 	    if (en > 0.) {
 		r = en1 * (2. - p) - 2.;
 		s = b + en1 * q;
 		d = (en - 1. + c / en) / (r * r + s * s);
 		p = d * r;
 		q = d * s;
 		e = f + 1.;
 		f = p * e - g * q;
 		g = q * e + p * g;
 		en1 = en;
 		goto L220;
 	    }
 	    f = 1. + f;
 	    d = f * f + g * g;
 	    pa = f / d;
 	    qa = -g / d;
 	    d = nu + .5 - p;
 	    q += ex;
 	    pa1 = (pa * q - qa * d) / ex;
 	    qa1 = (qa * q + pa * d) / ex;
 	    b = ex - M_PI_2 * (nu + .5);
 	    c = cos(b);
 	    s = sin(b);
 	    d = M_SQRT_2dPI / sqrt(ex);
 	    ya = d * (pa * s + qa * c);
 	    ya1 = d * (qa1 * s - pa1 * c);
 	} else { /* x > thresh_BESS_Y */
 	    /* ----------------------------------------------------------
 	       Use Campbell's asymptotic scheme.
 	       ---------------------------------------------------------- */
 	    na = 0;
 	    d1 = trunc(ex / fivpi);
 	    i = (int) d1;
 	    dmu = ex - 15. * d1 - d1 * pim5 - (*alpha + .5) * M_PI_2;
 	    if (i - (i / 2 << 1) == 0) {
 		cosmu = cos(dmu);
 		sinmu = sin(dmu);
 	    } else {
 		cosmu = -cos(dmu);
 		sinmu = -sin(dmu);
 	    }
 	    ddiv = 8. * ex;
 	    dmu = *alpha;
 	    den = sqrt(ex);
 	    for (k = 1; k <= 2; ++k) {
 		p = cosmu;
 		cosmu = sinmu;
 		sinmu = -p;
 		d1 = (2. * dmu - 1.) * (2. * dmu + 1.);
 		d2 = 0.;
 		div = ddiv;
 		p = 0.;
 		q = 0.;
 		q0 = d1 / div;
 		term = q0;
 		for (i = 2; i <= 20; ++i) {
 		    d2 += 8.;
 		    d1 -= d2;
 		    div += ddiv;
 		    term = -term * d1 / div;
 		    p += term;
 		    d2 += 8.;
 		    d1 -= d2;
 		    div += ddiv;
 		    term *= (d1 / div);
 		    q += term;
 		    if (fabs(term) <= DBL_EPSILON) {
 			break;
 		    }
 		}
 		p += 1.;
 		q += q0;
 		if (k == 1)
 		    ya = M_SQRT_2dPI * (p * cosmu - q * sinmu) / den;
 		else
 		    ya1 = M_SQRT_2dPI * (p * cosmu - q * sinmu) / den;
 		dmu += 1.;
 	    }
 	}
 	if (na == 1) {
 	    h = 2. * (nu + 1.) / ex;
 	    if (h > 1.) {
 		if (fabs(ya1) > DBL_MAX / h) {
 		    h = 0.;
 		    ya = 0.;
 		}
 	    }
 	    h = h * ya1 - ya;
 	    ya = ya1;
 	    ya1 = h;
 	}

 	/* ---------------------------------------------------------------
 	   Now have first one or two Y's
 	   --------------------------------------------------------------- */
 	by[0] = ya;
 	*ncalc = 1;
 	if(*nb > 1) {
 	    by[1] = ya1;
 	    if (ya1 != 0.) {
 		aye = 1. + *alpha;
 		twobyx = 2. / ex;
 		*ncalc = 2;
 		for (i = 2; i < *nb; ++i) {
 		    if (twobyx < 1.) {
 			if (fabs(by[i - 1]) * twobyx >= DBL_MAX / aye)
 			    goto L450;
 		    } else {
 			if (fabs(by[i - 1]) >= DBL_MAX / aye / twobyx)
 			    goto L450;
 		    }
 		    by[i] = twobyx * aye * by[i - 1] - by[i - 2];
 		    aye += 1.;
 		    ++(*ncalc);
 		}
 	    }
 	}
 L450:
 	for (i = *ncalc; i < *nb; ++i)
 	    by[i] = ML_NEGINF;/* was 0 */

     } else {
 	by[0] = 0.;
 	*ncalc = min0(*nb,0) - 1;
     }
 }
	/*
	* Mathlib : A C Library of Special Functions
	* Copyright (C) 1998-2015 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/rybesl 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 Y_bessel(double x, double alpha, int *nb,
	double by, int ncalc);

	// unused now from R
	double bessel_y(double x, double alpha)
	{
	int nb, ncalc;
	double na, *by;
	#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_y");
	return ML_NAN;
	}
	na = floor(alpha);
	if (alpha < 0) {
	/* Using Abramowitz & Stegun 9.1.2
	* this may not be quite optimal (CPU and accuracy wise) */
	return(((alpha - na == 0.5) ? 0 : bessel_y(x, -alpha) * cospi(alpha)) -
	((alpha == na ) ? 0 : bessel_j(x, -alpha) * sinpi(alpha)));
	}
	else if (alpha > 1e7) {
	MATHLIB_WARNING(_("besselY(x, nu): nu=%g too large for bessel_y() algorithm"),
	alpha);
	return ML_NAN;
	}
	nb = 1+ (int)na;/* nb-1 <= alpha < nb */
	alpha -= (double)(nb-1);
	#ifdef MATHLIB_STANDALONE
	by = (double *) calloc(nb, sizeof(double));
	if (!by) MATHLIB_ERROR("%s", _("bessel_y allocation error"));
	#else
	vmax = vmaxget();
	by = (double *) R_alloc((size_t) nb, sizeof(double));
	#endif
	Y_bessel(&x, &alpha, &nb, by, &ncalc);
	if(ncalc != nb) {/* error input */
	if(ncalc == -1) {
	#ifdef MATHLIB_STANDALONE
	free(by);
	#else
	vmaxset(vmax);
	#endif
	return ML_POSINF;
	}
	else if(ncalc < -1)
	MATHLIB_WARNING4(_("bessel_y(%g): ncalc (=%d) != nb (=%d); alpha=%g. Arg. out of range?\n"),
	x, ncalc, nb, alpha);
	else /* ncalc >= 0 */
	MATHLIB_WARNING2(_("bessel_y(%g,nu=%g): precision lost in result\n"),
	x, alpha+(double)nb-1);
	}
	x = by[nb-1];
	#ifdef MATHLIB_STANDALONE
	free(by);
	#else
	vmaxset(vmax);
	#endif
	return x;
	}

	/* Called from R: modified version of bessel_y(), accepting a work array
	* instead of allocating one. */
	double bessel_y_ex(double x, double alpha, double *by)
	{
	int nb, ncalc;
	double na;

	#ifdef IEEE_754
	/* NaNs propagated correctly */
	if (ISNAN(x) \|\| ISNAN(alpha)) return x + alpha;
	#endif
	if (x < 0) {
	ML_ERROR(ME_RANGE, "bessel_y");
	return ML_NAN;
	}
	na = floor(alpha);
	if (alpha < 0) {
	/* Using Abramowitz & Stegun 9.1.2
	* this may not be quite optimal (CPU and accuracy wise) */
	return(((alpha - na == 0.5) ? 0 : bessel_y_ex(x, -alpha, by) * cospi(alpha)) -
	((alpha == na ) ? 0 : bessel_j_ex(x, -alpha, by) * sinpi(alpha)));
	}
	else if (alpha > 1e7) {
	MATHLIB_WARNING(_("besselY(x, nu): nu=%g too large for bessel_y() algorithm"),
	alpha);
	return ML_NAN;
	}
	nb = 1+ (int)na;/* nb-1 <= alpha < nb */
	alpha -= (double)(nb-1);
	Y_bessel(&x, &alpha, &nb, by, &ncalc);
	if(ncalc != nb) {/* error input */
	if(ncalc == -1)
	return ML_POSINF;
	else if(ncalc < -1)
	MATHLIB_WARNING4(_("bessel_y(%g): ncalc (=%d) != nb (=%d); alpha=%g. Arg. out of range?\n"),
	x, ncalc, nb, alpha);
	else /* ncalc >= 0 */
	MATHLIB_WARNING2(_("bessel_y(%g,nu=%g): precision lost in result\n"),
	x, alpha+(double)nb-1);
	}
	x = by[nb-1];
	return x;
	}

	static void Y_bessel(double x, double alpha, int *nb,
	double by, int ncalc)
	{
	/* ----------------------------------------------------------------------

	This routine calculates Bessel functions Y_(N+ALPHA) (X)
	v for non-negative argument X, and non-negative order N+ALPHA.


	Explanation of variables in the calling sequence

	X - Non-negative argument for which
	Y's are to be calculated.
	ALPHA - Fractional part of order for which
	Y's 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).
	BY - Output vector of length NB. If the
	routine terminates normally (NCALC=NB), the vector BY
	contains the functions Y(ALPHA,X), ... , Y(NB-1+ALPHA,X),
	If (0 < NCALC < NB), BY(I) contains correct function
	values for I <= NCALC, and contains the ratios
	Y(ALPHA+I-1,X)/Y(ALPHA+I-2,X) for the rest of the array.
	NCALC - Output variable indicating possible errors.
	Before using the vector BY, 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 Y'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. In this case, BY[0] = 0.0, the remainder of the
	BY-vector is not calculated, and NCALC is set to
	MIN0(NB,0)-2 so that NCALC != NB.
	NCALC = -1: Y(ALPHA,X) >= XINF. The requested function
	values are set to 0.0.
	1 < NCALC < NB: Not all requested function values could
	be calculated accurately. BY(I) contains correct function
	values for I <= NCALC, and and the remaining NB-NCALC
	array elements contain 0.0.


	Intrinsic functions required are:

	DBLE, EXP, INT, MAX, MIN, REAL, SQRT


	Acknowledgement

	This program draws heavily on Temme's Algol program for Y(a,x)
	and Y(a+1,x) and on Campbell's programs for Y_nu(x). Temme's
	scheme is used for x < THRESH, and Campbell's scheme is used
	in the asymptotic region. Segments of code from both sources
	have been translated into Fortran 77, merged, and heavily modified.
	Modifications include parameterization of machine dependencies,
	use of a new approximation for ln(gamma(x)), and built-in
	protection against over/underflow.

	References: "Bessel functions J_nu(x) and Y_nu(x) of float
	order and float argument," Campbell, J. B.,
	Comp. Phy. Comm. 18, 1979, pp. 133-142.

	"On the numerical evaluation of the ordinary
	Bessel function of the second kind," Temme,
	N. M., J. Comput. Phys. 21, 1976, pp. 343-350.

	Latest modification: March 19, 1990

	Modified by: W. J. Cody
	Applied Mathematics Division
	Argonne National Laboratory
	Argonne, IL 60439
	----------------------------------------------------------------------*/


	/* ----------------------------------------------------------------------
	Mathematical constants
	FIVPI = 5*PI
	PIM5 = 5*PI - 15
	----------------------------------------------------------------------*/
	const static double fivpi = 15.707963267948966192;
	const static double pim5 = .70796326794896619231;

	/*----------------------------------------------------------------------
	Coefficients for Chebyshev polynomial expansion of
	1/gamma(1-x), abs(x) <= .5
	----------------------------------------------------------------------*/
	const static double ch[21] = { -6.7735241822398840964e-24,
	-6.1455180116049879894e-23,2.9017595056104745456e-21,
	1.3639417919073099464e-19,2.3826220476859635824e-18,
	-9.0642907957550702534e-18,-1.4943667065169001769e-15,
	-3.3919078305362211264e-14,-1.7023776642512729175e-13,
	9.1609750938768647911e-12,2.4230957900482704055e-10,
	1.7451364971382984243e-9,-3.3126119768180852711e-8,
	-8.6592079961391259661e-7,-4.9717367041957398581e-6,
	7.6309597585908126618e-5,.0012719271366545622927,
	.0017063050710955562222,-.07685284084478667369,
	-.28387654227602353814,.92187029365045265648 };

	/* Local variables */
	int i, k, na;

	double alfa, div, ddiv, even, gamma, term, cosmu, sinmu,
	b, c, d, e, f, g, h, p, q, r, s, d1, d2, q0, pa,pa1, qa,qa1,
	en, en1, nu, ex, ya,ya1, twobyx, den, odd, aye, dmu, x2, xna;

	en1 = ya = ya1 = 0; /* -Wall */

	ex = *x;
	nu = *alpha;
	if (*nb > 0 && 0. <= nu && nu < 1.) {
	if(ex < DBL_MIN \|\| ex > xlrg_BESS_Y) {
	/* Warning is not really appropriate, give
	* proper limit:
	* ML_ERROR(ME_RANGE, "Y_bessel"); */
	ncalc = nb;
	if(ex > xlrg_BESS_Y) by[0]= 0.; /was ML_POSINF /
	else if(ex < DBL_MIN) by[0]=ML_NEGINF;
	for(i=0; i < *nb; i++)
	by[i] = by[0];
	return;
	}
	xna = trunc(nu + .5);
	na = (int) xna;
	if (na == 1) {/* <==> .5 <= alpha < 1 <==> -5. <= nu < 0 /
	nu -= xna;
	}
	if (nu == -.5) {
	p = M_SQRT_2dPI / sqrt(ex);
	ya = p * sin(ex);
	ya1 = -p * cos(ex);
	} else if (ex < 3.) {
	/* -------------------------------------------------------------
	Use Temme's scheme for small X
	------------------------------------------------------------- */
	b = ex * .5;
	d = -log(b);
	f = nu * d;
	e = pow(b, -nu);
	if (fabs(nu) < M_eps_sinc)
	c = M_1_PI;
	else
	c = nu / sinpi(nu);

	/* ------------------------------------------------------------
	Computation of sinh(f)/f
	------------------------------------------------------------ */
	if (fabs(f) < 1.) {
	x2 = f * f;
	en = 19.;
	s = 1.;
	for (i = 1; i <= 9; ++i) {
	s = s * x2 / en / (en - 1.) + 1.;
	en -= 2.;
	}
	} else {
	s = (e - 1. / e) * .5 / f;
	}
	/* --------------------------------------------------------
	Computation of 1/gamma(1-a) using Chebyshev polynomials */
	x2 = nu * nu * 8.;
	aye = ch[0];
	even = 0.;
	alfa = ch[1];
	odd = 0.;
	for (i = 3; i <= 19; i += 2) {
	even = -(aye + aye + even);
	aye = -even * x2 - aye + ch[i - 1];
	odd = -(alfa + alfa + odd);
	alfa = -odd * x2 - alfa + ch[i];
	}
	even = (even * .5 + aye) * x2 - aye + ch[20];
	odd = (odd + alfa) * 2.;
	gamma = odd * nu + even;
	/* End of computation of 1/gamma(1-a)
	----------------------------------------------------------- */
	g = e * gamma;
	e = (e + 1. / e) * .5;
	f = 2. * c * (odd * e + even * s * d);
	e = nu * nu;
	p = g * c;
	q = M_1_PI / g;
	c = nu * M_PI_2;
	if (fabs(c) < M_eps_sinc)
	r = 1.;
	else
	r = sinpi(nu/2) / c;

	r = M_PI * c * r * r;
	c = 1.;
	d = -b * b;
	h = 0.;
	ya = f + r * q;
	ya1 = p;
	en = 1.;

	while (fabs(g / (1. + fabs(ya))) +
	fabs(h / (1. + fabs(ya1))) > DBL_EPSILON) {
	f = (f * en + p + q) / (en * en - e);
	c *= (d / en);
	p /= en - nu;
	q /= en + nu;
	g = c * (f + r * q);
	h = c * p - en * g;
	ya += g;
	ya1+= h;
	en += 1.;
	}
	ya = -ya;
	ya1 = -ya1 / b;
	} else if (ex < thresh_BESS_Y) {
	/* --------------------------------------------------------------
	Use Temme's scheme for moderate X : 3 <= x < 16
	-------------------------------------------------------------- */
	c = (.5 - nu) * (.5 + nu);
	b = ex + ex;
	e = ex * M_1_PI * cospi(nu) / DBL_EPSILON;
	e *= e;
	p = 1.;
	q = -ex;
	r = 1. + ex * ex;
	s = r;
	en = 2.;
	while (r * en * en < e) {
	en1 = en + 1.;
	d = (en - 1. + c / en) / s;
	p = (en + en - p * d) / en1;
	q = (-b + q * d) / en1;
	s = p * p + q * q;
	r *= s;
	en = en1;
	}
	f = p / s;
	p = f;
	g = -q / s;
	q = g;
	L220:
	en -= 1.;
	if (en > 0.) {
	r = en1 * (2. - p) - 2.;
	s = b + en1 * q;
	d = (en - 1. + c / en) / (r * r + s * s);
	p = d * r;
	q = d * s;
	e = f + 1.;
	f = p * e - g * q;
	g = q * e + p * g;
	en1 = en;
	goto L220;
	}
	f = 1. + f;
	d = f * f + g * g;
	pa = f / d;
	qa = -g / d;
	d = nu + .5 - p;
	q += ex;
	pa1 = (pa * q - qa * d) / ex;
	qa1 = (qa * q + pa * d) / ex;
	b = ex - M_PI_2 * (nu + .5);
	c = cos(b);
	s = sin(b);
	d = M_SQRT_2dPI / sqrt(ex);
	ya = d * (pa * s + qa * c);
	ya1 = d * (qa1 * s - pa1 * c);
	} else { /* x > thresh_BESS_Y */
	/* ----------------------------------------------------------
	Use Campbell's asymptotic scheme.
	---------------------------------------------------------- */
	na = 0;
	d1 = trunc(ex / fivpi);
	i = (int) d1;
	dmu = ex - 15. * d1 - d1 * pim5 - (alpha + .5) M_PI_2;
	if (i - (i / 2 << 1) == 0) {
	cosmu = cos(dmu);
	sinmu = sin(dmu);
	} else {
	cosmu = -cos(dmu);
	sinmu = -sin(dmu);
	}
	ddiv = 8. * ex;
	dmu = *alpha;
	den = sqrt(ex);
	for (k = 1; k <= 2; ++k) {
	p = cosmu;
	cosmu = sinmu;
	sinmu = -p;
	d1 = (2. * dmu - 1.) * (2. * dmu + 1.);
	d2 = 0.;
	div = ddiv;
	p = 0.;
	q = 0.;
	q0 = d1 / div;
	term = q0;
	for (i = 2; i <= 20; ++i) {
	d2 += 8.;
	d1 -= d2;
	div += ddiv;
	term = -term * d1 / div;
	p += term;
	d2 += 8.;
	d1 -= d2;
	div += ddiv;
	term *= (d1 / div);
	q += term;
	if (fabs(term) <= DBL_EPSILON) {
	break;
	}
	}
	p += 1.;
	q += q0;
	if (k == 1)
	ya = M_SQRT_2dPI * (p * cosmu - q * sinmu) / den;
	else
	ya1 = M_SQRT_2dPI * (p * cosmu - q * sinmu) / den;
	dmu += 1.;
	}
	}
	if (na == 1) {
	h = 2. * (nu + 1.) / ex;
	if (h > 1.) {
	if (fabs(ya1) > DBL_MAX / h) {
	h = 0.;
	ya = 0.;
	}
	}
	h = h * ya1 - ya;
	ya = ya1;
	ya1 = h;
	}

	/* ---------------------------------------------------------------
	Now have first one or two Y's
	--------------------------------------------------------------- */
	by[0] = ya;
	*ncalc = 1;
	if(*nb > 1) {
	by[1] = ya1;
	if (ya1 != 0.) {
	aye = 1. + *alpha;
	twobyx = 2. / ex;
	*ncalc = 2;
	for (i = 2; i < *nb; ++i) {
	if (twobyx < 1.) {
	if (fabs(by[i - 1]) * twobyx >= DBL_MAX / aye)
	goto L450;
	} else {
	if (fabs(by[i - 1]) >= DBL_MAX / aye / twobyx)
	goto L450;
	}
	by[i] = twobyx * aye * by[i - 1] - by[i - 2];
	aye += 1.;
	++(*ncalc);
	}
	}
	}
	L450:
	for (i = ncalc; i < nb; ++i)
	by[i] = ML_NEGINF;/* was 0 */

	} else {
	by[0] = 0.;
	ncalc = min0(nb,0) - 1;
	}
	}