uclibc-ng/libm/float/ynf.c

/*							ynf.c
 *
 *	Bessel function of second kind of integer order
 *
 *
 *
 * SYNOPSIS:
 *
 * float x, y, ynf();
 * int n;
 *
 * y = ynf( n, x );
 *
 *
 *
 * DESCRIPTION:
 *
 * Returns Bessel function of order n, where n is a
 * (possibly negative) integer.
 *
 * The function is evaluated by forward recurrence on
 * n, starting with values computed by the routines
 * y0() and y1().
 *
 * If n = 0 or 1 the routine for y0 or y1 is called
 * directly.
 *
 *
 *
 * ACCURACY:
 *
 *
 *  Absolute error, except relative when y > 1:
 *                      
 * arithmetic   domain     # trials      peak         rms
 *    IEEE      0, 30       10000       2.3e-6      3.4e-7
 *
 *
 * ERROR MESSAGES:
 *
 *   message         condition      value returned
 * yn singularity   x = 0              MAXNUMF
 * yn overflow                         MAXNUMF
 *
 * Spot checked against tables for x, n between 0 and 100.
 *
 */

/*
Cephes Math Library Release 2.2: June, 1992
Copyright 1984, 1987, 1992 by Stephen L. Moshier
Direct inquiries to 30 Frost Street, Cambridge, MA 02140
*/

#include <math.h>
extern float MAXNUMF, MAXLOGF;

float y0f(float), y1f(float), logf(float);

float ynf( int nn, float xx )
{
float x, an, anm1, anm2, r, xinv;
int k, n, sign;

x = xx;
n = nn;
if( n < 0 )
	{
	n = -n;
	if( (n & 1) == 0 )	/* -1**n */
		sign = 1;
	else
		sign = -1;
	}
else
	sign = 1;


if( n == 0 )
	return( sign * y0f(x) );
if( n == 1 )
	return( sign * y1f(x) );

/* test for overflow */
if( x <= 0.0 )
	{
	mtherr( "ynf", SING );
	return( -MAXNUMF );
	}
if( (x < 1.0) || (n > 29) )
	{
	an = (float )n;
	r = an * logf( an/x );
	if( r > MAXLOGF )
		{
		mtherr( "ynf", OVERFLOW );
		return( -MAXNUMF );
		}
	}

/* forward recurrence on n */

anm2 = y0f(x);
anm1 = y1f(x);
k = 1;
r = 2 * k;
xinv = 1.0/x;
do
	{
	an = r * anm1 * xinv  -  anm2;
	anm2 = anm1;
	anm1 = an;
	r += 2.0;
	++k;
	}
while( k < n );


return( sign * an );
}