/*							qine.c
 *
 *	Modified Bessel function I of noninteger order
 *      Exponentially scaled
 *
 *
 * SYNOPSIS:
 *
 * int qine( v, x, y );
 * QELT *v, *x, *y;
 *
 * qine( v, x, y );
 *
 *
 *
 * DESCRIPTION:
 *
 * Returns modified Bessel function of order v of the
 * argument.
 *
 * The power series is
 *
 *                inf      2   k
 *              v  -     (z /4)
 * I (z) = (z/2)   >  --------------
 *  v              -       -
 *                k=0  k! | (v+k+1)
 *
 *
 * For large x,
 *                                    2          2       2
 *             exp(z)            u - 1     (u - 1 )(u - 3 )
 * I (z)  =  ------------ { 1 - -------- + ---------------- + ...}
 *  v        sqrt(2 pi z)              1              2
 *                              1! (8z)        2! (8z)
 *
 * asymptotically, where
 * 
 *            2
 *     u = 4 v .
 * 
 *
 * x <= 0 is not supported.
 *
 * Series expansion is truncated at less than full working precision.
 *
 */

/*
Cephes Math Library Release 2.1:  November, 1988
Copyright 1984, 1987, 1988 by Stephen L. Moshier
Direct inquiries to 30 Frost Street, Cambridge, MA 02140
*/

/*							qin.c	*/
/*
 * Modified Bessel function of first kind of order n
 */

#include <stdio.h>
#include "qhead.h"
extern QELT qone[], qtwo[], qpi[];
static QELT qspi[NQ];
static QELT nkf[NQ];
static QELT nk1f[NQ];
static QELT pk[NQ];
static QELT fn[NQ];
static QELT s[NQ];
static QELT pn[NQ];
static QELT z0[NQ];
static QELT t1[NQ];
static int qsflg = 0;
int qgamma();

int qin( n, x, y )
QELT n[], x[], y[];
{
QELT t[NQ], u[NQ], z[NQ], k[NQ], ans[NQ];
long i, ln;
double dx;

if( qsflg == 0 )
	{ /* compute 1/sqrt(2 pi) */
	qsflg = 1;
	qmov( qpi, qspi );
	qspi[1] += 1;
	qdiv( qspi, qone, qspi );
	qsqrt( qspi, qspi );
	}

qtoe( x, (unsigned short *) &dx );
if( dx > 40.0 )
	goto asymp;

qmul( x, x, z );	/*z = -x * x / 4.0;		*/

if( z[1] < 3 )		/* x = 0, n = 0 is special case	*/
	{
	if( n[1] == 0 )
		qmov( qone, y );
/*
	else
		qclear( y );
*/
	return(0);
	}

if( z[1] > 1 )
	z[1] -= 2;

if( n[1] < 3 )
	qmov( qone, u );
else
	{
	qmov( n, ans );
	qadd( qone, ans, ans );
	qgamma( ans, u );
	qdiv( u, qone, u );		/*u = 1.0/gamma(n+1);*/
	}
qmov( u, ans );		/*ans = u;*/
qmov( qone, k );	/*k = 1.0;*/

while( ((int) ans[1] - (int) u[1]) < NBITS/2 )	/* 70 */
	{
	qadd( n, k, t );
	qmul( t, k, t );	/*u *= z / (k * (n+k));*/
	qdiv( t, z, t );
	qmul( t, u, u );
	if( (int) u[1] <= 0 )
		printf( "qin overflow\n" );
	qadd( u, ans, ans );	/*ans += u;*/
	qadd( qone, k, k );	/*k += 1.0;*/
	}

/* ans *= exp( n * log( x/2.0 ) );*/

if( n[1] < 3 )
	{
	qmov( ans, y );
	goto done;
	}
qmov( x, t );
t[1] -= 1;
qlog( t, u );
qmul( u, n, u );
qexp( u, t );
qmul( ans, t, y );

done:
/*
qexp( x, t );
qdiv( t, y, y );
qsqrt( x, t );
qmul( x, y, y );
*/
return 0;


/* Asymptotic expansion for In(x) */

asymp:

qtoe( n, (unsigned short *) &dx );		/* ln = n */
ln = dx;
dx = 4.0 * dx * dx;		/* pn = 4.0 * n * n; */
etoq( (unsigned short *) &dx, pn );
qmov( qone, pk );		/* pk = 1.0; */
qmov( x, z0 );			/* z0 = -8.0 * x; */
z0[1] += 3;
qneg( z0 );
qmov( qone, fn );		/* fn = 1.0; */
qmov( qone, t );		/* t = 1.0; */
qmov( t, s );			/* s = t; */
qmov( qone, nkf );		/* nkf = MAXNUM; */
nkf[1] += 16000;
i = 0;
do
	{
	qmul( pk, pk, t1 );	/* z = pn - pk * pk; */
	qsub( t1, pn, z );
	qmul( fn, z0, t1 );	/* t = t * z /(fn * z0); */
	qdiv( t1, z, t1 );
	qmul( t, t1, t );
	qmov( t, nk1f );	/* nk1f = fabs(t); */
	nk1f[0] = 0;
	qsub( nkf, nk1f, t1 );
	if( (i >= ln) && (t1[0] == 0) ) /* nk1f > nkf */
		{
		printf( "qin: i=%ld, %d\n", i, t[1]-s[1] );
		goto adone;
		}
	qmov( nk1f, nkf );	/* nkf = nk1f; */
	qadd( s, t, s );	/* s += t; */
	qadd( qone, fn, fn );	/* fn += 1.0; */
	qadd( qtwo, pk, pk );	/* pk += 2.0; */
	i += 1;
	}
while( ((int) s[1] - (int) t[1]) < NBITS/2 );	/* fabs(t/s) > MACHEP ); */

adone:

qmul( qspi, s, y );	/* ans = s exp(x) / sqrt( 2 pi x ) */

qexp( x, t );
qsqrt( x, z );
qdiv( z, t, t );
qmul( t, y, y );
return 0;
}
