lib/libm/igami.c

   1 /*                                                      igami()
   2  *
   3  *      Inverse of complemented imcomplete gamma integral
   4  *
   5  *
   6  *
   7  * SYNOPSIS:
   8  *
   9  * double a, x, p, igami();
  10  *
  11  * x = igami( a, p );
  12  *
  13  * DESCRIPTION:
  14  *
  15  * Given p, the function finds x such that
  16  *
  17  *  igamc( a, x ) = p.
  18  *
  19  * Starting with the approximate value
  20  *
  21  *         3
  22  *  x = a t
  23  *
  24  *  where
  25  *
  26  *  t = 1 - d - ndtri(p) sqrt(d)
  27  *
  28  * and
  29  *
  30  *  d = 1/9a,
  31  *
  32  * the routine performs up to 10 Newton iterations to find the
  33  * root of igamc(a,x) - p = 0.
  34  *
  35  * ACCURACY:
  36  *
  37  * Tested at random a, p in the intervals indicated.
  38  *
  39  *                a        p                      Relative error:
  40  * arithmetic   domain   domain     # trials      peak         rms
  41  *    IEEE     0.5,100   0,0.5       100000       1.0e-14     1.7e-15
  42  *    IEEE     0.01,0.5  0,0.5       100000       9.0e-14     3.4e-15
  43  *    IEEE    0.5,10000  0,0.5        20000       2.3e-13     3.8e-14
  44  */
  45 \f
  46 /*
  47 Cephes Math Library Release 2.8:  June, 2000
  48 Copyright 1984, 1987, 1995, 2000 by Stephen L. Moshier
  49 */
  50
  51 #include "mconf.h"
  52
  53 extern double MACHEP, MAXNUM, MAXLOG, MINLOG;
  54 #ifdef ANSIPROT
  55 extern double igamc ( double, double );
  56 extern double ndtri ( double );
  57 extern double exp ( double );
  58 extern double fabs ( double );
  59 extern double log ( double );
  60 extern double sqrt ( double );
  61 extern double lgam ( double );
  62 #else
  63 double igamc(), ndtri(), exp(), fabs(), log(), sqrt(), lgam();
  64 #endif
  65
  66 double igami( a, y0 )
  67 double a, y0;
  68 {
  69 double x0, x1, x, yl, yh, y, d, lgm, dithresh;
  70 int i, dir;
  71
  72 /* bound the solution */
  73 x0 = MAXNUM;
  74 yl = 0;
  75 x1 = 0;
  76 yh = 1.0;
  77 dithresh = 5.0 * MACHEP;
  78
  79 /* approximation to inverse function */
  80 d = 1.0/(9.0*a);
  81 y = ( 1.0 - d - ndtri(y0) * sqrt(d) );
  82 x = a * y * y * y;
  83
  84 lgm = lgam(a);
  85
  86 for( i=0; i<10; i++ )
  87         {
  88         if( x > x0 || x < x1 )
  89                 goto ihalve;
  90         y = igamc(a,x);
  91         if( y < yl || y > yh )
  92                 goto ihalve;
  93         if( y < y0 )
  94                 {
  95                 x0 = x;
  96                 yl = y;
  97                 }
  98         else
  99                 {
 100                 x1 = x;
 101                 yh = y;
 102                 }
 103 /* compute the derivative of the function at this point */
 104         d = (a - 1.0) * log(x) - x - lgm;
 105         if( d < -MAXLOG )
 106                 goto ihalve;
 107         d = -exp(d);
 108 /* compute the step to the next approximation of x */
 109         d = (y - y0)/d;
 110         if( fabs(d/x) < MACHEP )
 111                 goto done;
 112         x = x - d;
 113         }
 114
 115 /* Resort to interval halving if Newton iteration did not converge. */
 116 ihalve:
 117
 118 d = 0.0625;
 119 if( x0 == MAXNUM )
 120         {
 121         if( x <= 0.0 )
 122                 x = 1.0;
 123         while( x0 == MAXNUM )
 124                 {
 125                 x = (1.0 + d) * x;
 126                 y = igamc( a, x );
 127                 if( y < y0 )
 128                         {
 129                         x0 = x;
 130                         yl = y;
 131                         break;
 132                         }
 133                 d = d + d;
 134                 }
 135         }
 136 d = 0.5;
 137 dir = 0;
 138
 139 for( i=0; i<400; i++ )
 140         {
 141         x = x1  +  d * (x0 - x1);
 142         y = igamc( a, x );
 143         lgm = (x0 - x1)/(x1 + x0);
 144         if( fabs(lgm) < dithresh )
 145                 break;
 146         lgm = (y - y0)/y0;
 147         if( fabs(lgm) < dithresh )
 148                 break;
 149         if( x <= 0.0 )
 150                 break;
 151         if( y >= y0 )
 152                 {
 153                 x1 = x;
 154                 yh = y;
 155                 if( dir < 0 )
 156                         {
 157                         dir = 0;
 158                         d = 0.5;
 159                         }
 160                 else if( dir > 1 )
 161                         d = 0.5 * d + 0.5;
 162                 else
 163                         d = (y0 - yl)/(yh - yl);
 164                 dir += 1;
 165                 }
 166         else
 167                 {
 168                 x0 = x;
 169                 yl = y;
 170                 if( dir > 0 )
 171                         {
 172                         dir = 0;
 173                         d = 0.5;
 174                         }
 175                 else if( dir < -1 )
 176                         d = 0.5 * d;
 177                 else
 178                         d = (y0 - yl)/(yh - yl);
 179                 dir -= 1;
 180                 }
 181         }
 182 if( x == 0.0 )
 183         mtherr( "igami", UNDERFLOW );
 184
 185 done:
 186 return( x );
 187 }