1 files changed, 197 insertions, 0 deletions
diff --git a/libm/float/incbif.c b/libm/float/incbif.c
new file mode 100644
index 000000000..4d8c0652e
--- /dev/null
+++ b/libm/float/incbif.c
@@ -0,0 +1,197 @@
+/*							incbif()
+ *
+ *      Inverse of imcomplete beta integral
+ *
+ *
+ *
+ * SYNOPSIS:
+ *
+ * float a, b, x, y, incbif();
+ *
+ * x = incbif( a, b, y );
+ *
+ *
+ *
+ * DESCRIPTION:
+ *
+ * Given y, the function finds x such that
+ *
+ *  incbet( a, b, x ) = y.
+ *
+ * the routine performs up to 10 Newton iterations to find the
+ * root of incbet(a,b,x) - y = 0.
+ *
+ *
+ * ACCURACY:
+ *
+ *                      Relative error:
+ *                x     a,b
+ * arithmetic   domain  domain  # trials    peak       rms
+ *    IEEE      0,1     0,100     5000     2.8e-4    8.3e-6
+ *
+ * Overflow and larger errors may occur for one of a or b near zero
+ *  and the other large.
+ */
+
+
+/*
+Cephes Math Library Release 2.2:  July, 1992
+Copyright 1984, 1987, 1992 by Stephen L. Moshier
+Direct inquiries to 30 Frost Street, Cambridge, MA 02140
+*/
+
+#include <math.h>
+
+extern float MACHEPF, MINLOGF;
+
+#define fabsf(x) ( (x) < 0 ? -(x) : (x) )
+
+#ifdef ANSIC
+float incbetf(float, float, float);
+float ndtrif(float), expf(float), logf(float), sqrtf(float), lgamf(float);
+#else
+float incbetf();
+float ndtrif(), expf(), logf(), sqrtf(), lgamf();
+#endif
+
+float incbif( float aaa, float bbb, float yyy0 )
+{
+float aa, bb, yy0, a, b, y0;
+float d, y, x, x0, x1, lgm, yp, di;
+int i, rflg;
+
+
+aa = aaa;
+bb = bbb;
+yy0 = yyy0;
+if( yy0 <= 0 )
+	return(0.0);
+if( yy0 >= 1.0 )
+	return(1.0);
+
+/* approximation to inverse function */
+
+yp = -ndtrif(yy0);
+
+if( yy0 > 0.5 )
+	{
+	rflg = 1;
+	a = bb;
+	b = aa;
+	y0 = 1.0 - yy0;
+	yp = -yp;
+	}
+else
+	{
+	rflg = 0;
+	a = aa;
+	b = bb;
+	y0 = yy0;
+	}
+
+
+if( (aa <= 1.0) || (bb <= 1.0) )
+	{
+	y = 0.5 * yp * yp;
+	}
+else
+	{
+	lgm = (yp * yp - 3.0)* 0.16666666666666667;
+	x0 = 2.0/( 1.0/(2.0*a-1.0)  +  1.0/(2.0*b-1.0) );
+	y = yp * sqrtf( x0 + lgm ) / x0
+		- ( 1.0/(2.0*b-1.0) - 1.0/(2.0*a-1.0) )
+		* (lgm + 0.833333333333333333 - 2.0/(3.0*x0));
+	y = 2.0 * y;
+	if( y < MINLOGF )
+		{
+		x0 = 1.0;
+		goto under;
+		}
+	}
+
+x = a/( a + b * expf(y) );
+y = incbetf( a, b, x );
+yp = (y - y0)/y0;
+if( fabsf(yp) < 0.1 )
+	goto newt;
+
+/* Resort to interval halving if not close enough */
+x0 = 0.0;
+x1 = 1.0;
+di = 0.5;
+
+for( i=0; i<20; i++ )
+	{
+	if( i != 0 )
+		{
+		x = di * x1  + (1.0-di) * x0;
+		y = incbetf( a, b, x );
+		yp = (y - y0)/y0;
+		if( fabsf(yp) < 1.0e-3 )
+			goto newt;
+		}
+
+	if( y < y0 )
+		{
+		x0 = x;
+		di = 0.5;
+		}
+	else
+		{
+		x1 = x;
+		di *= di;
+		if( di == 0.0 )
+			di = 0.5;
+		}
+	}
+
+if( x0 == 0.0 )
+	{
+under:
+	mtherr( "incbif", UNDERFLOW );
+	goto done;
+	}
+
+newt:
+
+x0 = x;
+lgm = lgamf(a+b) - lgamf(a) - lgamf(b);
+
+for( i=0; i<10; i++ )
+	{
+/* compute the function at this point */
+	if( i != 0 )
+		y = incbetf(a,b,x0);
+/* compute the derivative of the function at this point */
+	d = (a - 1.0) * logf(x0) + (b - 1.0) * logf(1.0-x0) + lgm;
+	if( d < MINLOGF )
+		{
+		x0 = 0.0;
+		goto under;
+		}
+	d = expf(d);
+/* compute the step to the next approximation of x */
+	d = (y - y0)/d;
+	x = x0;
+	x0 = x0 - d;
+	if( x0 <= 0.0 )
+		{
+		x0 = 0.0;
+		goto under;
+		}
+	if( x0 >= 1.0 )
+		{
+		x0 = 1.0;
+		goto under;
+		}
+	if( i < 2 )
+		continue;
+	if( fabsf(d/x0) < 256.0 * MACHEPF )
+		goto done;
+	}
+
+done:
+if( rflg )
+	x0 = 1.0 - x0;
+return( x0 );
+}