/*							igami()
 *
 *      Inverse of complemented imcomplete gamma integral
 *
 *
 *
 * SYNOPSIS:
 *
 * int qigami( a, p, x );
 * QELT *a, *p, *x;
 *
 * qigami( a, p, x );
 *
 * DESCRIPTION:
 *
 * The program refines an initial estimate generated by the
 * double precision routine igami to find the root of
 *
 *  igamc(a,x) - p = 0.
 *
 *
 * ACCURACY:
 *
 * Set to do just one Newton-Raphson iteration.
 *
 */

/*
Cephes Math Library Release 2.3:  March, 1995
Copyright 1984, 1987, 1995 by Stephen L. Moshier
*/

/* qigami.c */

#include "qhead.h"
extern QELT qone[];
static QELT d[NQ];
static QELT y[NQ];
static QELT x0[NQ];
static QELT t[NQ];
static QELT lgm[NQ];
int qlgam(), qigamc();

int qigami( a, y0, ans )
QELT a[], y0[], ans[];
{
double da, dy0, dx0;
int i;
double igami();

/* approximation to inverse function */
qtoe( a, (unsigned short *) &da );
qtoe( y0, (unsigned short *) &dy0 );
dx0 = igami( da, dy0 );
etoq( (unsigned short *) &dx0, x0 );

qlgam( a, lgm ); /* log(gamma(a)) */

for( i=0; i<1; i++ )
	{
/*
	if( x0 <= 0.0 )
		{
		printf( "x0= %.5E\n", x0 );
		goto loop;
		}
*/
/* compute the function at this point */
	qigamc( a, x0, y );
/* compute the derivative of the function at this point */
/*	d = (a - 1.0) * log(x0) - x0 - lgam(a); */
	qmov( a, d );
	qsub( qone, d, d );
	qlog( x0, t );
	qmul( t, d, d );
	qsub( x0, d, d );
	qsub( lgm, d, d );
	qexp( d, t );		/* d = -exp(d) */
	qneg(t);
/* compute the step to the next approximation of x */
	qsub( y0, y, y );	/* d = (y - y0)/d */
	qdiv( t, y, t );
	qsub( t, x0, x0 );	/* x0 = x0 - d */
	}
qmov( x0, ans );
return 0;
}
