home *** CD-ROM | disk | FTP | other *** search
/ Power-Programmierung / CD2.mdf / c / library / dos / math / cephes / cprob / incbi.c < prev    next >
Encoding:
C/C++ Source or Header  |  1992-11-17  |  2.8 KB  |  188 lines
  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 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.  *    DEC       0,1     0,100      2500    2.9e-13   5.9e-15
  31.  *    IEEE      0,1     0,100      5000    8.3e-13   1.7e-14
  32.  *
  33.  * Larger errors were observed for a near zero and b large.
  34.  */
  35.  
  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 double MACHEP, MAXNUM, MAXLOG, MINLOG;
  46.  
  47. double incbi( aa, bb, yy0 )
  48. double aa, bb, yy0;
  49. {
  50. double a, b, y0;
  51. double d, y, x, x0, x1, lgm, yp, di;
  52. int i, n, rflg;
  53. double incbet();
  54. double ndtri(), exp(), fabs(), log(), sqrt(), lgam();
  55.  
  56.  
  57. if( yy0 <= 0 )
  58.     return(0.0);
  59. if( yy0 >= 1.0 )
  60.     return(1.0);
  61.  
  62. /* approximation to inverse function */
  63.  
  64. yp = -ndtri(yy0);
  65.  
  66. if( yy0 > 0.5 )
  67.     {
  68.     rflg = 1;
  69.     a = bb;
  70.     b = aa;
  71.     y0 = 1.0 - yy0;
  72.     yp = -yp;
  73.     }
  74. else
  75.     {
  76.     rflg = 0;
  77.     a = aa;
  78.     b = bb;
  79.     y0 = yy0;
  80.     }
  81.  
  82.  
  83. if( (aa <= 1.0) || (bb <= 1.0) )
  84.     {
  85.     y = yp * yp / 2.0;
  86.     }
  87. else
  88.     {
  89.     lgm = (yp * yp - 3.0)/6.0;
  90.     x0 = 2.0/( 1.0/(2.0*a-1.0)  +  1.0/(2.0*b-1.0) );
  91.     y = yp * sqrt( x0 + lgm ) / x0
  92.         - ( 1.0/(2.0*b-1.0) - 1.0/(2.0*a-1.0) )
  93.         * (lgm + 5.0/6.0 - 2.0/(3.0*x0));
  94.     y = 2.0 * y;
  95.     if( y < MINLOG )
  96.         {
  97.         x0 = 1.0;
  98.         goto under;
  99.         }
  100.     }
  101.  
  102. x = a/( a + b * exp(y) );
  103. y = incbet( a, b, x );
  104. yp = (y - y0)/y0;
  105. if( fabs(yp) < 1.0e-1 )
  106.     goto newt;
  107.  
  108. /* Resort to interval halving if not close enough */
  109. x0 = 0.0;
  110. x1 = 1.0;
  111. di = 0.5;
  112.  
  113. for( i=0; i<20; i++ )
  114.     {
  115.     if( i != 0 )
  116.         {
  117.         x = di * x1  + (1.0-di) * x0;
  118.         y = incbet( a, b, x );
  119.         yp = (y - y0)/y0;
  120.         if( fabs(yp) < 1.0e-3 )
  121.             goto newt;
  122.         }
  123.  
  124.     if( y < y0 )
  125.         {
  126.         x0 = x;
  127.         di = 0.5;
  128.         }
  129.     else
  130.         {
  131.         x1 = x;
  132.         di *= di;
  133.         if( di == 0.0 )
  134.             di = 0.5;
  135.         }
  136.     }
  137.  
  138. if( x0 == 0.0 )
  139.     {
  140. under:
  141.     mtherr( "incbi", UNDERFLOW );
  142.     goto done;
  143.     }
  144.  
  145. newt:
  146.  
  147. x0 = x;
  148. lgm = lgam(a+b) - lgam(a) - lgam(b);
  149.  
  150. for( i=0; i<10; i++ )
  151.     {
  152. /* compute the function at this point */
  153.     if( i != 0 )
  154.         y = incbet(a,b,x0);
  155. /* compute the derivative of the function at this point */
  156.     d = (a - 1.0) * log(x0) + (b - 1.0) * log(1.0-x0) + lgm;
  157.     if( d < MINLOG )
  158.         {
  159.         x0 = 0.0;
  160.         goto under;
  161.         }
  162.     d = exp(d);
  163. /* compute the step to the next approximation of x */
  164.     d = (y - y0)/d;
  165.     x = x0;
  166.     x0 = x0 - d;
  167.     if( x0 <= 0.0 )
  168.         {
  169.         x0 = 0.0;
  170.         goto under;
  171.         }
  172.     if( x0 >= 1.0 )
  173.         {
  174.         x0 = 1.0;
  175.         goto under;
  176.         }
  177.     if( i < 2 )
  178.         continue;
  179.     if( fabs(d/x0) < 256.0 * MACHEP )
  180.         goto done;
  181.     }
  182.  
  183. done:
  184. if( rflg )
  185.     x0 = 1.0 - x0;
  186. return( x0 );
  187. }
  188.