/* igamil()
 *
 *      Inverse of complemented imcomplete gamma integral
 *
 *
 *
 * SYNOPSIS:
 *
 * long double a, x, y, igamil();
 *
 * x = igamil( a, y );
 *
 *
 *
 * DESCRIPTION:
 *
 * Given y, the function finds x such that
 *
 *  igamc( a, x ) = y.
 *
 * Starting with the approximate value
 *
 *         3
 *  x = a t
 *
 *  where
 *
 *  t = 1 - d - ndtri(y) sqrt(d)
 * 
 * and
 *
 *  d = 1/9a,
 *
 * the routine performs up to 10 Newton iterations to find the
 * root of igamc(a,x) - y = 0.
 *
 *
 * ACCURACY:
 *
 * Tested for a ranging from 0.5 to 30 and x from 0 to 0.5.
 *
 *                      Relative error:
 * arithmetic   domain     # trials      peak         rms
 *    DEC       0,0.5         3400       8.8e-16     1.3e-16
 *    IEEE      0,0.5        10000       1.1e-14     1.0e-15
 *
 */

/*
Cephes Math Library Release 2.3:  March, 1995
Copyright 1984, 1995 by Stephen L. Moshier
*/

#include "mconf.h"

extern long double MACHEPL, MAXNUML, MAXLOGL, MINLOGL;
#ifdef ANSIPROT
extern long double ndtril ( long double );
extern long double expl ( long double );
extern long double fabsl ( long double );
extern long double logl ( long double );
extern long double sqrtl ( long double );
extern long double lgaml ( long double );
extern long double igamcl ( long double, long double );
#else
long double ndtril(), expl(), fabsl(), logl(), sqrtl(), lgaml();
long double igamcl();
#endif

long double igamil( a, y0 )
long double a, y0;
{
long double x0, x1, x, yl, yh, y, d, lgm, dithresh;
int i, dir;

/* bound the solution */
x0 = MAXNUML;
yl = 0.0L;
x1 = 0.0L;
yh = 1.0L;
dithresh = 4.0 * MACHEPL;

/* approximation to inverse function */
d = 1.0L/(9.0L*a);
y = ( 1.0L - d - ndtril(y0) * sqrtl(d) );
x = a * y * y * y;

lgm = lgaml(a);

for( i=0; i<10; i++ )
 {
 if( x > x0 || x < x1 )
  goto ihalve;
 y = igamcl(a,x);
 if( y < yl || y > yh )
  goto ihalve;
 if( y < y0 )
  {
  x0 = x;
  yl = y;
  }
 else
  {
  x1 = x;
  yh = y;
  }
/* compute the derivative of the function at this point */
 d = (a - 1.0L) * logl(x0) - x0 - lgm;
 if( d < -MAXLOGL )
  goto ihalve;
 d = -expl(d);
/* compute the step to the next approximation of x */
 d = (y - y0)/d;
 x = x - d;
 if( i < 3 )
  continue;
 if( fabsl(d/x) < dithresh )
  goto done;
 }

/* Resort to interval halving if Newton iteration did not converge. */
ihalve:

d = 0.0625L;
if( x0 == MAXNUML )
 {
 if( x <= 0.0L )
  x = 1.0L;
 while( x0 == MAXNUML )
  {
  x = (1.0L + d) * x;
  y = igamcl( a, x );
  if( y < y0 )
   {
   x0 = x;
   yl = y;
   break;
   }
  d = d + d;
  }
 }
d = 0.5L;
dir = 0;

for( i=0; i<400; i++ )
 {
 x = x1  +  d * (x0 - x1);
 y = igamcl( a, x );
 lgm = (x0 - x1)/(x1 + x0);
 if( fabsl(lgm) < dithresh )
  break;
 lgm = (y - y0)/y0;
 if( fabsl(lgm) < dithresh )
  break;
 if( x <= 0.0L )
  break;
 if( y > y0 )
  {
  x1 = x;
  yh = y;
  if( dir < 0 )
   {
   dir = 0;
   d = 0.5L;
   }
  else if( dir > 1 )
   d = 0.5L * d + 0.5L; 
  else
   d = (y0 - yl)/(yh - yl);
  dir += 1;
  }
 else
  {
  x0 = x;
  yl = y;
  if( dir > 0 )
   {
   dir = 0;
   d = 0.5L;
   }
  else if( dir < -1 )
   d = 0.5L * d;
  else
   d = (y0 - yl)/(yh - yl);
  dir -= 1;
  }
 }
if( x == 0.0L )
 mtherr( "igamil", UNDERFLOW );

done:
return( x );
}