/*							qfresnl
 *
 *	Fresnel integral
 *
 *
 *
 * SYNOPSIS:
 *
 * int qfresnl( x, s, c );
 * QELT *x, *s, *c;
 *
 * qfresnl( x, s, c );
 *
 *
 * DESCRIPTION:
 *
 * Evaluates the Fresnel integrals
 *
 *           x
 *           -
 *          | |
 * C(x) =   |   cos(pi/2 t**2) dt,
 *        | |
 *         -
 *          0
 *
 *           x
 *           -
 *          | |
 * S(x) =   |   sin(pi/2 t**2) dt.
 *        | |
 *         -
 *          0
 *
 *
 * The integrals are evaluated by a power series for x < 1.
 * For large x auxiliary functions f(x) and g(x) are employed
 * such that
 *
 * C(x) = 0.5 + f(x) sin( pi/2 x**2 ) - g(x) cos( pi/2 x**2 )
 * S(x) = 0.5 - f(x) cos( pi/2 x**2 ) - g(x) sin( pi/2 x**2 )
 *
 * Routine qfresfg computes f and g.
 *
 *
 * ACCURACY:
 *
 * Series expansions are truncated at less than full working precision.
 */

/*
Cephes Math Library Release 2.1:  January, 1989
Copyright 1984, 1987, 1989 by Stephen L. Moshier
*/

#include <stdio.h>
#include "qhead.h"

/* asymptotic expansion converges to 144 bits at BRKPT >= 8.2 */
#if NBITS > 144
#define BRKPT 8.2
#else
#define BRKPT 5.7
#endif
extern QELT qpi[], qone[], qtwo[];

static double dtemp = 0.0;
static QELT x2[NQ];
static QELT c[NQ];
static QELT q[NQ];
static QELT d[NQ];
static QELT t[NQ];
static QELT x4[NQ];
static QELT x3[NQ];
static QELT t1[NQ];
static QELT cc0[NQ];
static QELT ss0[NQ];
static QELT f[NQ];
static QELT g[NQ];

#ifdef ANSIPROT
int qfresfg( QELT *, QELT *, QELT * );
#else
int qfresfg();
#endif

int qfresnl( x, ss, cc )
QELT x[], ss[], cc[];
{
QELT max;
int sign, i, iconv, ic2;
long il;

qtoe( x, (unsigned short *) &dtemp );
if( dtemp > BRKPT ) /* 7.0 5.7 */
	goto asymp0;

max = 0;
qmul( x, x, x2 );
qmul( x, x2, x3 );	/* x3 = x**3		*/
qmul( x2, x2, x4 );	/* x4 = x**4		*/
sign = -1;
qmov( qone, c );	/* c=1			*/
qmov( qpi, q );
q[1] -= 1;
qmul( q, x3, x3 );	/* x3 = pi/2 * x**3	*/
qmul( q, q, q );	/* q = (pi/2)**2	*/

qmov( qone, t );
qmov( qone, cc0 );
qmov( qone, ss0 );
il = 3;
ltoq( &il, d );
qdiv( d, ss0, ss0 );	/* ss = pi/2 * x**3 / 3	*/
i = 1;
do
	{
	qadd( c, qone, c );	/* c += 1		*/
	qdiv( c, t, t );	/* t /= c		*/
	qmul( q, t, t );	/* t *= (pi/2)**2	*/
	qmul( t, x4, t );	/* t *= x**4		*/
	il = 4*i + 1;
	ltoq( &il, d );
	qmov( t, f );
	qdiv( d, f, f );
	if( f[1] > max )
		max = f[1];
	if( sign < 0 )
		qsub( f, cc0, cc0 );
	else
		qadd( cc0, f, cc0 );
	iconv = cc0[1] - f[1];
	qadd( c, qone, c );
	qdiv( c, t, t );
	qmov( t, f );
	il = 4*i + 3;
	ltoq( &il, d );
	qdiv( d, f, f );
	if( sign < 0 )
		qsub( f, ss0, ss0 );
	else
		qadd( ss0, f, ss0 );
	if( f[1] > max )
		max = f[1];
	ic2 = (int) ss0[1] - (int) f[1];
	if( ic2 < iconv )
		iconv = ic2;
	sign = -sign;
	++i;
	}
while( iconv < 144 );

max = NBITS - max;
if( ss0[1] < cc0[1] )
	max += ss0[1];
else
	max += cc0[1];
if( max < (QELT) (3*NBITS/7) )
	printf( "qfresf: prec %d bits\n", max );

qmul( x, cc0, cc );		/*	C(x)	*/
qmul( x3, ss0, ss );		/*	S(x)	*/
goto done;


asymp0:

qfresfg( x, f, g );

qmov( qpi, t );		/*	t = PIO2 * x * x;	*/
t[1] -= 1;
qmul( x, t, t );
qmul( x, t, t );
qcos( t, c );		/*	c = cos(t);		*/
qsin( t, t );		/*	t = sin(t);		*/
qmov( qone, cc );	/*	cc = 0.5  +  f * t  -  g * c;	*/
cc[1] -= 1;
qmul( f, t, x3 );
qadd( x3, cc, cc );
qmul( g, c, x3 );
qsub( x3, cc, cc );

qmov( qone, ss );	/*	ss = 0.5  -  f * c  -  g * t;	*/
ss[1] -= 1;
qmul( f, c, x3 );
qsub( x3, ss, ss );
qmul( g, t, x3 );
qsub( x3, ss, ss );

done:

if( x[0] != 0 )
	{
	cc[0] = -1;
	ss[0] = -1;
	}

return(0);
}


/*		the auxiliary functions		*/
static QELT lss[NQ] = {0};
static QELT lcc[NQ] = {0};

int qfresfg( x, f, g )
QELT x[], f[], g[];
{
int sign, max;
#define co c
#define si q

qtoe( x, (unsigned short *) &dtemp );
if( dtemp <= BRKPT )
	goto nasymp;


/*		Asymptotic power series auxiliary functions
 *		for large argument
 */

qmov( qone, c );	/* 	c = 1.0;		*/
qmov( qone, t );	/*	t = 1.0;		*/
qadd( qone, qone, t1 );	/*	t1 = 2.0;		*/
qmul( qpi, x, q );	/*	q = PIO2 * 2.0 * x;	*/
qmul( q, x, x3 );	/*	x3 = nn * x;		*/
qmul( x3, x3, x4 );	/*	x4 = x3 * x3;		*/
sign = -1;
qmov( qone, f );	/*	f = 1.0;		*/
qmov( qone, g );	/*	g = 1.0;		*/

/*
if( dtemp > 20000.0 )
		goto done5;
*/

while( ((int) qone[1] - (int) t[1]) < NBITS/2 )
{
qadd( qtwo, c, c );	/*	c += 2.0;		*/
qmul( t, c, t );	/*	t *= c / x4;		*/
qdiv( x4, t, t );
if( sign < 0 )
	qsub( t, f, f ); /*	f -= t;			*/
else
	qadd( t, f, f ); /*	f += t;			*/
qadd( qtwo, c, c );	/*	c += 2.0;		*/
qmul( t, c, t );	/*	t *= c;			*/
/*
	qtoe( t, (unsigned short *) &dtemp );
	printf( "%.4E\n", dtemp );
*/
if( t[1] > t1[1] )
	{
	max = qone[1] - t1[1];
	if( max < 3*NBITS/7 )
		printf( "qfres asymp prec %d bits\n", max );
	goto done5;
	}
qmov( t, t1 );		/*	t1 = t;			*/
if( sign < 0 )
	qsub( t, g, g ); /*	g -= t;			*/
else
	qadd( t, g, g ); /*	g += t;			*/
sign = -sign;
}


done5:
qmul( x3, q, d );	/*	g /= x3 * q;		*/
qdiv( d, g, g );
qdiv( q, f, f );	/*	f /= q;		*/
return 0;





nasymp:

qfresnl( x, lss, lcc );
qmov( qone, co );
co[1] -= 1;
qsub( lcc, co, lcc );
qsub( lss, co, lss );

qmul( qpi, x2, x2 );
x2[1] -= 1;
qcos( x2, co );
qsin( x2, si );

		/* f = (.5-ss)*cos - (.5-cc)*sin	*/

qmul( lss, co, f );
qmul( lcc, si, t );
qsub( t, f, f );

		/* g = (.5-cc)*cos + (.5-ss)*sin	*/
qmul( lcc, co, g );
qmul( lss, si, t );
qadd( t, g, g );
return 0;
}

