/*							qerfc.c	*/
/* error function check routine */

/* loop counts adjusted for full convergence to 9 word mantissa */

#include "qhead.h"

extern QELT oneopi[], qsqrt2[], qone[], qtwo[];

/* exp(x**2) * erfc(x) */

int qerfc( x, y )
QELT *x, *y;
{
QELT xx[NQ], z[NQ], a[NQ], b[NQ];
long i, j;
int sign;

qmov( x, xx );
qmul( xx, xx, z );	/* square of y				*/

sign = xx[0];

if( xx[1]< (QELT) (EXPONE + 2) ) /* xx < 4.0 */
{
i = 54; /*27;*/
if( xx[1] < EXPONE ) /* xx < 1.0 */
	i = 20; /*10;*/
j = 4 * i + 1;
ltoq( &j, a);		/*  2 * i  +  1			*/
qmov( a, y );

/* continued fraction expansion */
while( j > 1 )
	{
	qsub( qone, a, a );
	qmuli( a, z, b );
	qdiv( y, b, y );
	qsub( qone, a, a );
	qadd( a, y, y );

	qsub( qone, a, a );
	qmuli( a, z, b );
	qdiv( y, b, y );
	qsub( qone, a, a );
	qsub( y, a, y );
	j -= 4;
	}
qdiv( y, xx, y );

qneg(z);		/* exp( -xx**2 )	*/
qexp( z, a );
qmul( a, y, y );
qmul( oneopi, y, y );
qsub( y, qone, y );
/*qdiv( a, y, y );*/
}


/* alternate continued fraction expansion for large x */
else
{
j = 172;
if( xx[1] > (QELT) (EXPONE + 1) )  /* xx > 2.0 */
	j = 120; /*76;*/
ltoq( &j, a );
if( sign != 0 )
	xx[0] = 0;
qmov( xx, y );
while( j > 0 )
	{
	qdiv( y, a, y );
	*(y+1) -= 1;
	qadd( xx, y, y );
	qsub( qone, a, a );
	j -= 1;
	}
qdiv( y, qone, y );

qneg(z);		/* exp( -xx**2 )	*/
qexp( z, a );
qmul( a, y, y );
qmul( oneopi, y, y );
y[1] -= 1;
if( sign != 0 )
	qsub( y, qtwo, y );
}
return 0;
}