lib/libm/incbif.c

   1 /*                                                      incbif()
   2  *
   3  *      Inverse of imcomplete beta integral
   4  *
   5  *
   6  *
   7  * SYNOPSIS:
   8  *
   9  * float a, b, x, y, incbif();
  10  *
  11  * x = incbif( 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 up to 10 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     0,100     5000     2.8e-4    8.3e-6
  31  *
  32  * Overflow and larger errors may occur for one of a or b near zero
  33  *  and the other large.
  34  */
  35 \f
  36
  37 /*
  38 Cephes Math Library Release 2.2:  July, 1992
  39 Copyright 1984, 1987, 1992 by Stephen L. Moshier
  40 Direct inquiries to 30 Frost Street, Cambridge, MA 02140
  41 */
  42
  43 #include "mconf.h"
  44
  45 extern float MACHEPF, MINLOGF;
  46
  47 #define fabsf(x) ( (x) < 0 ? -(x) : (x) )
  48
  49 #ifdef ANSIC
  50 float incbetf(float, float, float);
  51 float ndtrif(float), expf(float), logf(float), sqrtf(float), lgamf(float);
  52 #else
  53 float incbetf();
  54 float ndtrif(), expf(), logf(), sqrtf(), lgamf();
  55 #endif
  56
  57 #ifdef ANSIC
  58 float incbif( float aaa, float bbb, float yyy0 )
  59 #else
  60 float incbif( aaa, bbb, yyy0 )
  61 double aaa, bbb, yyy0;
  62 #endif
  63 {
  64 float aa, bb, yy0, a, b, y0;
  65 float d, y, x, x0, x1, lgm, yp, di;
  66 int i, rflg;
  67
  68
  69 aa = aaa;
  70 bb = bbb;
  71 yy0 = yyy0;
  72 if( yy0 <= 0 )
  73         return(0.0);
  74 if( yy0 >= 1.0 )
  75         return(1.0);
  76
  77 /* approximation to inverse function */
  78
  79 yp = -ndtrif(yy0);
  80
  81 if( yy0 > 0.5 )
  82         {
  83         rflg = 1;
  84         a = bb;
  85         b = aa;
  86         y0 = 1.0 - yy0;
  87         yp = -yp;
  88         }
  89 else
  90         {
  91         rflg = 0;
  92         a = aa;
  93         b = bb;
  94         y0 = yy0;
  95         }
  96
  97
  98 if( (aa <= 1.0) || (bb <= 1.0) )
  99         {
 100         y = 0.5 * yp * yp;
 101         }
 102 else
 103         {
 104         lgm = (yp * yp - 3.0)* 0.16666666666666667;
 105         x0 = 2.0/( 1.0/(2.0*a-1.0)  +  1.0/(2.0*b-1.0) );
 106         y = yp * sqrtf( x0 + lgm ) / x0
 107                 - ( 1.0/(2.0*b-1.0) - 1.0/(2.0*a-1.0) )
 108                 * (lgm + 0.833333333333333333 - 2.0/(3.0*x0));
 109         y = 2.0 * y;
 110         if( y < MINLOGF )
 111                 {
 112                 x0 = 1.0;
 113                 goto under;
 114                 }
 115         }
 116
 117 x = a/( a + b * expf(y) );
 118 y = incbetf( a, b, x );
 119 yp = (y - y0)/y0;
 120 if( fabsf(yp) < 0.1 )
 121         goto newt;
 122
 123 /* Resort to interval halving if not close enough */
 124 x0 = 0.0;
 125 x1 = 1.0;
 126 di = 0.5;
 127
 128 for( i=0; i<20; i++ )
 129         {
 130         if( i != 0 )
 131                 {
 132                 x = di * x1  + (1.0-di) * x0;
 133                 y = incbetf( a, b, x );
 134                 yp = (y - y0)/y0;
 135                 if( fabsf(yp) < 1.0e-3 )
 136                         goto newt;
 137                 }
 138
 139         if( y < y0 )
 140                 {
 141                 x0 = x;
 142                 di = 0.5;
 143                 }
 144         else
 145                 {
 146                 x1 = x;
 147                 di *= di;
 148                 if( di == 0.0 )
 149                         di = 0.5;
 150                 }
 151         }
 152
 153 if( x0 == 0.0 )
 154         {
 155 under:
 156         mtherr( "incbif", UNDERFLOW );
 157         goto done;
 158         }
 159
 160 newt:
 161
 162 x0 = x;
 163 lgm = lgamf(a+b) - lgamf(a) - lgamf(b);
 164
 165 for( i=0; i<10; i++ )
 166         {
 167 /* compute the function at this point */
 168         if( i != 0 )
 169                 y = incbetf(a,b,x0);
 170 /* compute the derivative of the function at this point */
 171         d = (a - 1.0) * logf(x0) + (b - 1.0) * logf(1.0-x0) + lgm;
 172         if( d < MINLOGF )
 173                 {
 174                 x0 = 0.0;
 175                 goto under;
 176                 }
 177         d = expf(d);
 178 /* compute the step to the next approximation of x */
 179         d = (y - y0)/d;
 180         x = x0;
 181         x0 = x0 - d;
 182         if( x0 <= 0.0 )
 183                 {
 184                 x0 = 0.0;
 185                 goto under;
 186                 }
 187         if( x0 >= 1.0 )
 188                 {
 189                 x0 = 1.0;
 190                 goto under;
 191                 }
 192         if( i < 2 )
 193                 continue;
 194         if( fabsf(d/x0) < 256.0 * MACHEPF )
 195                 goto done;
 196         }
 197
 198 done:
 199 if( rflg )
 200         x0 = 1.0 - x0;
 201 return( x0 );
 202 }