lib/libm/incbi.c

   1 /*                                                      incbi()
   2  *
   3  *      Inverse of imcomplete beta integral
   4  *
   5  *
   6  *
   7  * SYNOPSIS:
   8  *
   9  * double a, b, x, y, incbi();
  10  *
  11  * x = incbi( a, b, y );
  12  *
  13  *
  14  *
  15  * DESCRIPTION:
  16  *
  17  * Given y, the function finds x such that
  18  *
  19  *  incbet( a, b, x ) = y .
  20  *
  21  * The routine performs interval halving or Newton iterations to find the
  22  * root of incbet(a,b,x) - y = 0.
  23  *
  24  *
  25  * ACCURACY:
  26  *
  27  *                      Relative error:
  28  *                x     a,b
  29  * arithmetic   domain  domain  # trials    peak       rms
  30  *    IEEE      0,1    .5,10000   50000    5.8e-12   1.3e-13
  31  *    IEEE      0,1   .25,100    100000    1.8e-13   3.9e-15
  32  *    IEEE      0,1     0,5       50000    1.1e-12   5.5e-15
  33  *    VAX       0,1    .5,100     25000    3.5e-14   1.1e-15
  34  * With a and b constrained to half-integer or integer values:
  35  *    IEEE      0,1    .5,10000   50000    5.8e-12   1.1e-13
  36  *    IEEE      0,1    .5,100    100000    1.7e-14   7.9e-16
  37  * With a = .5, b constrained to half-integer or integer values:
  38  *    IEEE      0,1    .5,10000   10000    8.3e-11   1.0e-11
  39  */
  40 \f
  41
  42 /*
  43 Cephes Math Library Release 2.8:  June, 2000
  44 Copyright 1984, 1996, 2000 by Stephen L. Moshier
  45 */
  46
  47 #include "mconf.h"
  48
  49 extern double MACHEP, MAXNUM, MAXLOG, MINLOG;
  50 #ifdef ANSIPROT
  51 extern double ndtri ( double );
  52 extern double exp ( double );
  53 extern double fabs ( double );
  54 extern double log ( double );
  55 extern double sqrt ( double );
  56 extern double lgam ( double );
  57 extern double incbet ( double, double, double );
  58 #else
  59 double ndtri(), exp(), fabs(), log(), sqrt(), lgam(), incbet();
  60 #endif
  61
  62 double incbi( aa, bb, yy0 )
  63 double aa, bb, yy0;
  64 {
  65 double a, b, y0, d, y, x, x0, x1, lgm, yp, di, dithresh, yl, yh, xt;
  66 int i, rflg, dir, nflg;
  67
  68
  69 i = 0;
  70 if( yy0 <= 0 )
  71         return(0.0);
  72 if( yy0 >= 1.0 )
  73         return(1.0);
  74 x0 = 0.0;
  75 yl = 0.0;
  76 x1 = 1.0;
  77 yh = 1.0;
  78 nflg = 0;
  79
  80 if( aa <= 1.0 || bb <= 1.0 )
  81         {
  82         dithresh = 1.0e-6;
  83         rflg = 0;
  84         a = aa;
  85         b = bb;
  86         y0 = yy0;
  87         x = a/(a+b);
  88         y = incbet( a, b, x );
  89         goto ihalve;
  90         }
  91 else
  92         {
  93         dithresh = 1.0e-4;
  94         }
  95 /* approximation to inverse function */
  96
  97 yp = -ndtri(yy0);
  98
  99 if( yy0 > 0.5 )
 100         {
 101         rflg = 1;
 102         a = bb;
 103         b = aa;
 104         y0 = 1.0 - yy0;
 105         yp = -yp;
 106         }
 107 else
 108         {
 109         rflg = 0;
 110         a = aa;
 111         b = bb;
 112         y0 = yy0;
 113         }
 114
 115 lgm = (yp * yp - 3.0)/6.0;
 116 x = 2.0/( 1.0/(2.0*a-1.0)  +  1.0/(2.0*b-1.0) );
 117 d = yp * sqrt( x + lgm ) / x
 118         - ( 1.0/(2.0*b-1.0) - 1.0/(2.0*a-1.0) )
 119         * (lgm + 5.0/6.0 - 2.0/(3.0*x));
 120 d = 2.0 * d;
 121 if( d < MINLOG )
 122         {
 123         x = 1.0;
 124         goto under;
 125         }
 126 x = a/( a + b * exp(d) );
 127 y = incbet( a, b, x );
 128 yp = (y - y0)/y0;
 129 if( fabs(yp) < 0.2 )
 130         goto newt;
 131
 132 /* Resort to interval halving if not close enough. */
 133 ihalve:
 134
 135 dir = 0;
 136 di = 0.5;
 137 for( i=0; i<100; i++ )
 138         {
 139         if( i != 0 )
 140                 {
 141                 x = x0  +  di * (x1 - x0);
 142                 if( x == 1.0 )
 143                         x = 1.0 - MACHEP;
 144                 if( x == 0.0 )
 145                         {
 146                         di = 0.5;
 147                         x = x0  +  di * (x1 - x0);
 148                         if( x == 0.0 )
 149                                 goto under;
 150                         }
 151                 y = incbet( a, b, x );
 152                 yp = (x1 - x0)/(x1 + x0);
 153                 if( fabs(yp) < dithresh )
 154                         goto newt;
 155                 yp = (y-y0)/y0;
 156                 if( fabs(yp) < dithresh )
 157                         goto newt;
 158                 }
 159         if( y < y0 )
 160                 {
 161                 x0 = x;
 162                 yl = y;
 163                 if( dir < 0 )
 164                         {
 165                         dir = 0;
 166                         di = 0.5;
 167                         }
 168                 else if( dir > 3 )
 169                         di = 1.0 - (1.0 - di) * (1.0 - di);
 170                 else if( dir > 1 )
 171                         di = 0.5 * di + 0.5;
 172                 else
 173                         di = (y0 - y)/(yh - yl);
 174                 dir += 1;
 175                 if( x0 > 0.75 )
 176                         {
 177                         if( rflg == 1 )
 178                                 {
 179                                 rflg = 0;
 180                                 a = aa;
 181                                 b = bb;
 182                                 y0 = yy0;
 183                                 }
 184                         else
 185                                 {
 186                                 rflg = 1;
 187                                 a = bb;
 188                                 b = aa;
 189                                 y0 = 1.0 - yy0;
 190                                 }
 191                         x = 1.0 - x;
 192                         y = incbet( a, b, x );
 193                         x0 = 0.0;
 194                         yl = 0.0;
 195                         x1 = 1.0;
 196                         yh = 1.0;
 197                         goto ihalve;
 198                         }
 199                 }
 200         else
 201                 {
 202                 x1 = x;
 203                 if( rflg == 1 && x1 < MACHEP )
 204                         {
 205                         x = 0.0;
 206                         goto done;
 207                         }
 208                 yh = y;
 209                 if( dir > 0 )
 210                         {
 211                         dir = 0;
 212                         di = 0.5;
 213                         }
 214                 else if( dir < -3 )
 215                         di = di * di;
 216                 else if( dir < -1 )
 217                         di = 0.5 * di;
 218                 else
 219                         di = (y - y0)/(yh - yl);
 220                 dir -= 1;
 221                 }
 222         }
 223 mtherr( "incbi", PLOSS );
 224 if( x0 >= 1.0 )
 225         {
 226         x = 1.0 - MACHEP;
 227         goto done;
 228         }
 229 if( x <= 0.0 )
 230         {
 231 under:
 232         mtherr( "incbi", UNDERFLOW );
 233         x = 0.0;
 234         goto done;
 235         }
 236
 237 newt:
 238
 239 if( nflg )
 240         goto done;
 241 nflg = 1;
 242 lgm = lgam(a+b) - lgam(a) - lgam(b);
 243
 244 for( i=0; i<8; i++ )
 245         {
 246         /* Compute the function at this point. */
 247         if( i != 0 )
 248                 y = incbet(a,b,x);
 249         if( y < yl )
 250                 {
 251                 x = x0;
 252                 y = yl;
 253                 }
 254         else if( y > yh )
 255                 {
 256                 x = x1;
 257                 y = yh;
 258                 }
 259         else if( y < y0 )
 260                 {
 261                 x0 = x;
 262                 yl = y;
 263                 }
 264         else
 265                 {
 266                 x1 = x;
 267                 yh = y;
 268                 }
 269         if( x == 1.0 || x == 0.0 )
 270                 break;
 271         /* Compute the derivative of the function at this point. */
 272         d = (a - 1.0) * log(x) + (b - 1.0) * log(1.0-x) + lgm;
 273         if( d < MINLOG )
 274                 goto done;
 275         if( d > MAXLOG )
 276                 break;
 277         d = exp(d);
 278         /* Compute the step to the next approximation of x. */
 279         d = (y - y0)/d;
 280         xt = x - d;
 281         if( xt <= x0 )
 282                 {
 283                 y = (x - x0) / (x1 - x0);
 284                 xt = x0 + 0.5 * y * (x - x0);
 285                 if( xt <= 0.0 )
 286                         break;
 287                 }
 288         if( xt >= x1 )
 289                 {
 290                 y = (x1 - x) / (x1 - x0);
 291                 xt = x1 - 0.5 * y * (x1 - x);
 292                 if( xt >= 1.0 )
 293                         break;
 294                 }
 295         x = xt;
 296         if( fabs(d/x) < 128.0 * MACHEP )
 297                 goto done;
 298         }
 299 /* Did not converge.  */
 300 dithresh = 256.0 * MACHEP;
 301 goto ihalve;
 302
 303 done:
 304
 305 if( rflg )
 306         {
 307         if( x <= MACHEP )
 308                 x = 1.0 - MACHEP;
 309         else
 310                 x = 1.0 - x;
 311         }
 312 return( x );
 313 }