/*							qjn.c	*/
/* jn.c		1 Dec 83
 * Bessel function of order n
 */

#define DEBUG 0
#define ERRCK 0
#include <stdio.h>
#include "qhead.h"

extern QELT qhalf[], qone[], qtwo[], qpi[];

static QELT t[NQ];
static QELT u[NQ];
QELT hankzz[NQ];
static QELT k[NQ];
static QELT ans[NQ];
QELT hankc[NQ];
QELT hanks[NQ];
static QELT j[NQ];
static QELT m[NQ];
QELT hankpp[NQ];
QELT hankqq[NQ];
static QELT rans[NQ];
static QELT ru[NQ];

#ifdef ANSIPROT
int qhank( QELT *, QELT *, QELT * );
static int qrecur( QELT *, QELT *, QELT *, QELT * );
#else
double floor(), sqrt();
static int qrecur();
int mtherr(), qhank(), qgamma();
#endif

#define fabs(x) ((x) < 0 ? -(x) : (x))
#if DEBUG
static QELT yh[NQ];
#endif

int qjn( nn, xx, y )
QELT nn[], xx[], y[];
{
QELT n[NQ], x[NQ];
double dx, dn, an, du;
int i, sign;
QELT bt;

bt = 0;
qmov( nn, n );
qmov( xx, x );
qtoe( n, (unsigned short *) &dn );
qtoe( x, (unsigned short *) &dx );
sign = 1;
an = fabs( dn );
if( an == floor(an) )
	{
	i = an - 16384.0 * floor( an/16384.0 );
	if( n[0] != 0 )
		{
		if( i & 1 )
			sign = -sign;
		n[0] = 0;
		}
	if( x[0] != 0 )
		{
		if( i & 1 )
			sign = -sign;
		x[0] = 0;
		}
	}
else
	{
	if( x[0] != 0 )
		{
		mtherr( "qjv", DOMAIN );
		qclear( y );
		goto done;
		}
	}

dx = fabs(dx);


#if DEBUG
	qhank( n, x, yh);
#endif

#if NBITS > 336
if( dx > 81.0 )
#else
#if NBITS > 144
if( dx > 100.0 )
#else
if( dx > 34.0 )
#endif
#endif
	{
/*	if( dx > 0.95*an )*/
	if( dx > 1.4*an )
		{
		qhank( n, x, y );
		goto done;
		}
	if( dx > 0.7 * an )
		{
		du = 3.6 * sqrt(dx);
/*		du = 0.8 * dx;*/
		etoq( (unsigned short *) &du, ru );
		ru[0] = 0;
		qfloor( n, m );
		qsub( m, n, m );
		qfloor( ru, ru );
		qadd( ru, m, ru );
		if( n[0] == 0 )
			{
			qrecur( n, x, ru, rans );
			}
		else
			{
			qmov( ru, m );
			qmov( n, ru );
			qrecur( m, x, ru, rans );
			qmov( m, ru );
			}
		if( rans[1] == 0 )
			{
			qclear( y );
			goto done;
			}
		qhank( ru, x, y );
		if( n[0] == 0 )
			qdiv( rans, y, y );
		else
			qmul( rans, y, y );
		goto done;
		}
	}

qmul( x, x, hankzz );	/*z = -x * x / 4.0;		*/
/* x = 0, n = 0 is special case	*/

if( hankzz[1] < 3 )
	{
	if( n[1] < 3 )
		qmov( qone, y );
	else
		qclear( y );
	goto done;
	}

hankzz[1] -= 2;
hankzz[0] = -1;

qmov( n, ans );
qadd( qone, ans, ans );
if( (n[1] < 3) || (qcmp( qone, n ) == 0) )
	qmov( qone, u );
else
	{
	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( u[1] > (QELT) (ans[1] - NBITS)
      || u[1] > (QELT) (qone[1] - NBITS) )
	{
	qadd( n, k, t );
	qmul( t, k, t );	/*u *= z / (k * (n+k));*/
	qdiv( t, hankzz, t );
	qmul( t, u, u );
#if 0
/* #if DEBUG */
	qtoe( u, (unsigned short *) &du );
	fprintf(stderr, "pseries term: %.5e\n", du );
#endif
/* remember largest term summed */
#if ERRCK
	if( u[1] > bt )
		bt = u[1];
	if( ans[1] > bt )
		bt = ans[1];
#endif
	qadd( u, ans, ans );	/*ans += u;*/
	qadd( qone, k, k );	/*k += 1.0;*/
	}

/* estimate cancellation error */
#if ERRCK
i = bt - ans[1];
if( i > NBITS/2
    || DEBUG )
	fprintf(stderr, "qjn pseries: %d bits cancellation\n", i );
#endif

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

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

#if DEBUG
qsub( y, yh, yh );
qtoe( yh, (unsigned short *) &du );
printf( "qjn - qhank = %.5e\n", du );
#endif
done:
if( sign < 0 )
	y[0] = ~y[0];
return(0);
}


/* Hankel's asymptotic expansion
 * for large x.
 * AMS55 #9.2.5.
 */

int qhank( n, x, y )
QELT n[], x[], y[];
{
QELT bt;
int flag, sign, nsum;

#if ERRCK || DEBUG
int    i;
double dconv;
#endif

bt = 0;
nsum = 0;
qmul( n, n, m );	/* m = 4.0*n*n;*/
m[1] += 2;
qmov( qone, j );	/* j = 1.0;*/
qmov( x, hankzz );	/* z = 8.0 * x;*/
hankzz[1] += 3;
qmov( qone, k );	/* k = 1.0;*/
qmov( qone, hankc );	/* hankc = 1.0;*/
qsub( qone, m, u );	/* u = (m - 1.0)/z;*/
qdiv( hankzz, u, u );
qmov( u, hanks );	/* hanks = u;*/
sign = 1;
qmov( qone, ans );	/* conv = 1.0;*/
flag = 0;
qmov( qone, t );	/* t = 1.0;*/

while( t[1] > (qone[1] - NBITS)
    || u[1] > (qone[1] - NBITS) )
	{
	qadd( qtwo, k, k );	/* k += 2.0;*/
	qadd( qone, j, j );	/* j += 1.0;*/
	sign = -sign;
	qmul( k, k, t );	/* u *= (m - k * k)/(j * z);*/
	qsub( t, m, t );
	qdiv( j, t, t );
	qdiv( hankzz, t, t );
	qmul( t, u, u );
	if( sign < 0 )		/* hankc += sign * u;*/
		qsub( u, hankc, hankc );
	else
		qadd( u, hankc, hankc );
/* remember largest term summed */
#if ERRCK
	if( u[1] > bt )
		bt = u[1];
	if( hankc[1] > bt )
		bt = hankc[1];
#endif
/*	printf( "Hank P: %.5E %.5E", u, p ); */
	qadd( qtwo, k, k );	/* k += 2.0;*/
	qadd( qone, j, j );	/* j += 1.0;*/
	qmul( k, k, t );	/* u *= (m - k * k)/(j * z);*/
	qsub( t, m, t );
	qdiv( j, t, t );
	qdiv( hankzz, t, t );
	qmul( t, u, u );
	if( sign < 0 )		/* q += sign * u;*/
		qsub( u, hanks, hanks );
	else
		qadd( u, hanks, hanks );
/* remember largest term summed */
#if ERRCK
	if( u[1] > bt )
		bt = u[1];
	if( hanks[1] > bt )
		bt = hanks[1];
#endif
/*	printf( " Q: %.5E %.5E\n", u, q ); */
	qdiv( hankc, u, t );	/* t = fabs(u/p);*/
	t[0] = 0;
	if( qcmp( t, ans ) < 0 )	/* ( t < conv )*/
		{
		qmov( t, ans );	/* conv = t;*/
		qmov( hanks, hankqq );	/* qq = hanks; */
		qmov( hankc, hankpp );	/* pp = hankc; */
		flag = 1;
		nsum += 1;
		}
/* stop if the terms start getting larger */
	else
		{
		if( flag != 0 )
			{
			goto hank1;
			}
		}
	}

hank1:

/* estimate cancellation error */
#if ERRCK
i = bt - hankpp[1];
if(i > NBITS/2
    || DEBUG )
  {
    qtoe( n, (unsigned short *) &dconv );
    printf( "qhank(%.5e,", dconv );
    qtoe( x, (unsigned short *) &dconv );
    printf( "%.5e): ", dconv );
    printf( "%d bits cancellation after %d terms\n", i, nsum );
  }
#endif

#if DEBUG
    qtoe( ans, (unsigned short *) &dconv );
    printf( "qhank: last term / sum = %.4E\n", dconv );
#endif

qmov( n, t );	/* u = x - (0.5*n + 0.25) * PI;*/
qmul( qhalf, t, t );
qmov( qone, hanks );
hanks[1] -= 2;
qadd( hanks, t, t );
qmul( qpi, t, t );
qsub( t, x, u );

/* t = sqrt( 2.0/(PI*x) ) * ( pp * cos(u) - qq * sin(u) ); */
qmul( qpi, x, t );
qdiv( t, qtwo, t );
qsqrt( t, hankzz );
qsin( u, hanks );
qcos( u, hankc );
qmul( hankc, hankpp, k );
qmul( hanks, hankqq, m );
qsub( m, k, k );
qmul( k, hankzz, y );
#if 0
/* #if DEBUG */
qtoe( y, (unsigned short *) &dconv );
printf( "%.4e\n", dconv );
#endif
return 0;
}



/* Reduce the order by backward recurrence.
 * AMS55 #9.1.27 and 9.1.73.
 */

static QELT pkm2[NQ];
static QELT pkm1[NQ];
static QELT pk[NQ];
static QELT pkp1[NQ];
static QELT qkm2[NQ];
static QELT qkm1[NQ];
static QELT qk[NQ];
static QELT xk[NQ];
static QELT yk[NQ];
static QELT r[NQ];
static QELT kf[NQ];
static QELT k1[NQ];

static int qrecur( n, x, newn, ans )
QELT n[], x[], newn[], ans[];
{
int nflag, ctr;
#if DEBUG
double da, db;
#endif

nflag = n[0];

/* continued fraction for Jn(x)/Jn-1(x)  */

/* fstart: */

#if DEBUG
qtoe( n, (unsigned short *) &da );
qtoe( newn, (unsigned short *) &db );
printf( "qn = %.6e, qnewn = %.6e, qcfrac = ", da, db );
#endif

qclear( pkm2 );
qmov( qone, qkm2 );
qmov( x, pkm1 );
qadd( n, n, qkm1 );
qmul( x, x, xk );   /* xk = -x * x; */
qneg( xk );
qmov( qkm1, yk );
qmov( qone, ans );
ctr = 0;
do
	{
	qadd( qtwo, yk, yk );
/*	pk = pkm1 * yk +  pkm2 * xk; */
	qmul( pkm1, yk, pk );
	qmul( pkm2, xk, t );
	qadd( t, pk, pk );
/*	qk = qkm1 * yk +  qkm2 * xk;*/
	qmul( qkm1, yk, qk );
	qmul( qkm2, xk, t );
	qadd( t, qk, qk );
	qmov( pkm1, pkm2 );
	qmov( pk, pkm1 );
	qmov( qkm1, qkm2 );
	qmov( qk, qkm1 );
	if( qk[1] != 0 )
		qdiv( qk, pk, r );
	if( r[1] != 0 )
		{
/*		t = fabs( (ans - r)/r ); */
		qsub( r, ans, t );
		qdiv( r, t, t );
		t[0] = 0;
		qmov( r, ans );
		}
	else
		qmov( qone, t );

	if( ++ctr > 1000 )
		{
		printf( "qrecur: continued fraction did not converge\n" );
		goto done;
		}
/*
	if( t < MACHEP )
		goto done;
*/

	if( pk[1] > (qone[1]+NBITS) )
		{
		pkm2[1] -= NBITS;
		pkm1[1] -= NBITS;
		qkm2[1] -= NBITS;
		qkm1[1] -= NBITS;
		}
	}
while( t[1] > (qone[1] - NBITS) );

done:

#if DEBUG
qtoe( ans, (unsigned short *) &da );
printf( "%.6e\n", da );
#endif

/* Change n to n-1 if n < 0 and the continued fraction is small
 */
/*
if( nflag )
	{
	if( fabs(ans) < 0.125 )
		{
		nflag = -1;
		*n = *n - 1.0;
		goto fstart;
		}
	}
*/

qmov( newn, kf );

/* backward recurrence
 *              2k
 *  J   (x)  =  --- J (x)  -  J   (x)
 *   k-1         x   k         k+1
 */

qmov( qone, pk );
qdiv( ans, qone, pkm1 );
qsub( qone, n, k1 );
qadd( k1, k1, r );
do
	{
/*	pkm2 = (pkm1 * r  -  pk * x) / x;*/
	qmul( pkm1, r, pkm2 );
	qmul( pk, x, t );
	qsub( t, pkm2, pkm2 );
	qdiv( x, pkm2, pkm2 );
	qmov( pk, pkp1 );
	qmov( pkm1, pk );
	qmov( pkm2, pkm1 );
	qsub( qtwo, r, r );
/*
	t = fabs(pkp1) + fabs(pk);
	if( (k > (kf + 2.5)) && (fabs(pkm1) < 0.25*t) )
		{
		k1 -= 1.0;
		t = x*x;
		pkm2 = ( (r*(r+2.0)-t)*pk - r*x*pkp1 )/t;
		pkp1 = pk;
		pk = pkm1;
		pkm1 = pkm2;
		r -= 2.0;
		}
*/
	qsub( qone, k1, k1 );
	qadd( qhalf, kf, t );
	}
while( qcmp( k1, t ) > 0 );

/* Take the larger of the last two iterates
 * on the theory that it may have less cancellation error.
 */
/*
if( cancel )
	{
	if( (kf >= 0.0) && (fabs(pk) > fabs(pkm1)) )
		{
		k += 1.0;
		pkm2 = pk;
		}
	}
*/
qmov( k1, newn );
#if DEBUG
qtoe( k1, (unsigned short *) &da );
qtoe( pkm2, (unsigned short *) &db );
printf( "qnewn %.6e qrans %.6e\n", da, db );
#endif
qmov( pkm2, ans );
return(0);
}