/* 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.
 *
 * It is valid in the right-hand tail of the distribution, p < 0.5.
 *
 * 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();

double igami();


int qigami( a, y0, ans )
QELT a[], y0[], ans[];
{
double da, dy0, dx0;
int i;
union
  {
    unsigned short s[4];
    double d;
  } dz;

/* approximation to inverse function */
qtoe( a, dz.s );
da = dz.d;

qtoe( y0, dz.s );
dy0 = dz.d;

dx0 = igami( da, dy0 );

dz.d = dx0;
etoq( dz.s, 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;
}