uClibc now has a math library. muahahahaha!

 -Erik

1 files changed, 424 insertions, 0 deletions

diff --git a/libm/float/incbetf.c b/libm/float/incbetf.c
new file mode 100644
index 000000000..fed9aae4b
--- /dev/null
+++ b/libm/float/incbetf.c
@@ -0,0 +1,424 @@
+/*							incbetf.c
+ *
+ *	Incomplete beta integral
+ *
+ *
+ * SYNOPSIS:
+ *
+ * float a, b, x, y, incbetf();
+ *
+ * y = incbetf( a, b, x );
+ *
+ *
+ * DESCRIPTION:
+ *
+ * Returns incomplete beta integral of the arguments, evaluated
+ * from zero to x.  The function is defined as
+ *
+ *                  x
+ *     -            -
+ *    | (a+b)      | |  a-1     b-1
+ *  -----------    |   t   (1-t)   dt.
+ *   -     -     | |
+ *  | (a) | (b)   -
+ *                 0
+ *
+ * The domain of definition is 0 <= x <= 1.  In this
+ * implementation a and b are restricted to positive values.
+ * The integral from x to 1 may be obtained by the symmetry
+ * relation
+ *
+ *    1 - incbet( a, b, x )  =  incbet( b, a, 1-x ).
+ *
+ * The integral is evaluated by a continued fraction expansion.
+ * If a < 1, the function calls itself recursively after a
+ * transformation to increase a to a+1.
+ *
+ * ACCURACY:
+ *
+ * Tested at random points (a,b,x) with a and b in the indicated
+ * interval and x between 0 and 1.
+ *
+ * arithmetic   domain     # trials      peak         rms
+ * Relative error:
+ *    IEEE       0,30       10000       3.7e-5      5.1e-6
+ *    IEEE       0,100      10000       1.7e-4      2.5e-5
+ * The useful domain for relative error is limited by underflow
+ * of the single precision exponential function.
+ * Absolute error:
+ *    IEEE       0,30      100000       2.2e-5      9.6e-7
+ *    IEEE       0,100      10000       6.5e-5      3.7e-6
+ *
+ * Larger errors may occur for extreme ratios of a and b.
+ *
+ * ERROR MESSAGES:
+ *   message         condition      value returned
+ * incbetf domain     x<0, x>1          0.0
+ */
+
+
+/*
+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>
+
+#ifdef ANSIC
+float lgamf(float), expf(float), logf(float);
+static float incbdf(float, float, float);
+static float incbcff(float, float, float);
+float incbpsf(float, float, float);
+#else
+float lgamf(), expf(), logf();
+float incbpsf();
+static float incbcff(), incbdf();
+#endif
+
+#define fabsf(x) ( (x) < 0 ? -(x) : (x) )
+
+/* BIG = 1/MACHEPF */
+#define BIG   16777216.
+extern float MACHEPF, MAXLOGF;
+#define MINLOGF (-MAXLOGF)
+
+float incbetf( float aaa, float bbb, float xxx )
+{
+float aa, bb, xx, ans, a, b, t, x, onemx;
+int flag;
+
+aa = aaa;
+bb = bbb;
+xx = xxx;
+if( (xx <= 0.0) || ( xx >= 1.0) )
+	{
+	if( xx == 0.0 )
+		return(0.0);
+	if( xx == 1.0 )
+		return( 1.0 );
+	mtherr( "incbetf", DOMAIN );
+	return( 0.0 );
+	}
+
+onemx = 1.0 - xx;
+
+
+/* transformation for small aa */
+
+if( aa <= 1.0 )
+	{
+	ans = incbetf( aa+1.0, bb, xx );
+	t = aa*logf(xx) + bb*logf( 1.0-xx )
+		+ lgamf(aa+bb) - lgamf(aa+1.0) - lgamf(bb);
+	if( t > MINLOGF )
+		ans += expf(t);
+	return( ans );
+	}
+
+
+/* see if x is greater than the mean */
+
+if( xx > (aa/(aa+bb)) )
+	{
+	flag = 1;
+	a = bb;
+	b = aa;
+	t = xx;
+	x = onemx;
+	}
+else
+	{
+	flag = 0;
+	a = aa;
+	b = bb;
+	t = onemx;
+	x = xx;
+	}
+
+/* transformation for small aa */
+/*
+if( a <= 1.0 )
+	{
+	ans = a*logf(x) + b*logf( onemx )
+		+ lgamf(a+b) - lgamf(a+1.0) - lgamf(b);
+ 	t = incbetf( a+1.0, b, x );
+	if( ans > MINLOGF )
+		t += expf(ans);
+	goto bdone;
+	}
+*/
+/* Choose expansion for optimal convergence */
+
+
+if( b > 10.0 )
+	{
+if( fabsf(b*x/a) < 0.3 )
+	{
+	t = incbpsf( a, b, x );
+	goto bdone;
+	}
+	}
+
+ans = x * (a+b-2.0)/(a-1.0);
+if( ans < 1.0 )
+	{
+	ans = incbcff( a, b, x );
+	t = b * logf( t );
+	}
+else
+	{
+	ans = incbdf( a, b, x );
+	t = (b-1.0) * logf(t);
+	}
+
+t += a*logf(x) + lgamf(a+b) - lgamf(a) - lgamf(b);
+t += logf( ans/a );
+
+if( t < MINLOGF )
+	{
+	t = 0.0;
+	if( flag == 0 )
+		{
+		mtherr( "incbetf", UNDERFLOW );
+		}
+	}
+else
+	{
+	t = expf(t);
+	}
+bdone:
+
+if( flag )
+	t = 1.0 - t;
+
+return( t );
+}
+
+/* Continued fraction expansion #1
+ * for incomplete beta integral
+ */
+
+static float incbcff( float aa, float bb, float xx )
+{
+float a, b, x, xk, pk, pkm1, pkm2, qk, qkm1, qkm2;
+float k1, k2, k3, k4, k5, k6, k7, k8;
+float r, t, ans;
+static float big = BIG;
+int n;
+
+a = aa;
+b = bb;
+x = xx;
+k1 = a;
+k2 = a + b;
+k3 = a;
+k4 = a + 1.0;
+k5 = 1.0;
+k6 = b - 1.0;
+k7 = k4;
+k8 = a + 2.0;
+
+pkm2 = 0.0;
+qkm2 = 1.0;
+pkm1 = 1.0;
+qkm1 = 1.0;
+ans = 1.0;
+r = 0.0;
+n = 0;
+do
+	{
+	
+	xk = -( x * k1 * k2 )/( k3 * k4 );
+	pk = pkm1 +  pkm2 * xk;
+	qk = qkm1 +  qkm2 * xk;
+	pkm2 = pkm1;
+	pkm1 = pk;
+	qkm2 = qkm1;
+	qkm1 = qk;
+
+	xk = ( x * k5 * k6 )/( k7 * k8 );
+	pk = pkm1 +  pkm2 * xk;
+	qk = qkm1 +  qkm2 * xk;
+	pkm2 = pkm1;
+	pkm1 = pk;
+	qkm2 = qkm1;
+	qkm1 = qk;
+
+	if( qk != 0 )
+		r = pk/qk;
+	if( r != 0 )
+		{
+		t = fabsf( (ans - r)/r );
+		ans = r;
+		}
+	else
+		t = 1.0;
+
+	if( t < MACHEPF )
+		goto cdone;
+
+	k1 += 1.0;
+	k2 += 1.0;
+	k3 += 2.0;
+	k4 += 2.0;
+	k5 += 1.0;
+	k6 -= 1.0;
+	k7 += 2.0;
+	k8 += 2.0;
+
+	if( (fabsf(qk) + fabsf(pk)) > big )
+		{
+		pkm2 *= MACHEPF;
+		pkm1 *= MACHEPF;
+		qkm2 *= MACHEPF;
+		qkm1 *= MACHEPF;
+		}
+	if( (fabsf(qk) < MACHEPF) || (fabsf(pk) < MACHEPF) )
+		{
+		pkm2 *= big;
+		pkm1 *= big;
+		qkm2 *= big;
+		qkm1 *= big;
+		}
+	}
+while( ++n < 100 );
+
+cdone:
+return(ans);
+}
+
+
+/* Continued fraction expansion #2
+ * for incomplete beta integral
+ */
+
+static float incbdf( float aa, float bb, float xx )
+{
+float a, b, x, xk, pk, pkm1, pkm2, qk, qkm1, qkm2;
+float k1, k2, k3, k4, k5, k6, k7, k8;
+float r, t, ans, z;
+static float big = BIG;
+int n;
+
+a = aa;
+b = bb;
+x = xx;
+k1 = a;
+k2 = b - 1.0;
+k3 = a;
+k4 = a + 1.0;
+k5 = 1.0;
+k6 = a + b;
+k7 = a + 1.0;;
+k8 = a + 2.0;
+
+pkm2 = 0.0;
+qkm2 = 1.0;
+pkm1 = 1.0;
+qkm1 = 1.0;
+z = x / (1.0-x);
+ans = 1.0;
+r = 0.0;
+n = 0;
+do
+	{
+	
+	xk = -( z * k1 * k2 )/( k3 * k4 );
+	pk = pkm1 +  pkm2 * xk;
+	qk = qkm1 +  qkm2 * xk;
+	pkm2 = pkm1;
+	pkm1 = pk;
+	qkm2 = qkm1;
+	qkm1 = qk;
+
+	xk = ( z * k5 * k6 )/( k7 * k8 );
+	pk = pkm1 +  pkm2 * xk;
+	qk = qkm1 +  qkm2 * xk;
+	pkm2 = pkm1;
+	pkm1 = pk;
+	qkm2 = qkm1;
+	qkm1 = qk;
+
+	if( qk != 0 )
+		r = pk/qk;
+	if( r != 0 )
+		{
+		t = fabsf( (ans - r)/r );
+		ans = r;
+		}
+	else
+		t = 1.0;
+
+	if( t < MACHEPF )
+		goto cdone;
+
+	k1 += 1.0;
+	k2 -= 1.0;
+	k3 += 2.0;
+	k4 += 2.0;
+	k5 += 1.0;
+	k6 += 1.0;
+	k7 += 2.0;
+	k8 += 2.0;
+
+	if( (fabsf(qk) + fabsf(pk)) > big )
+		{
+		pkm2 *= MACHEPF;
+		pkm1 *= MACHEPF;
+		qkm2 *= MACHEPF;
+		qkm1 *= MACHEPF;
+		}
+	if( (fabsf(qk) < MACHEPF) || (fabsf(pk) < MACHEPF) )
+		{
+		pkm2 *= big;
+		pkm1 *= big;
+		qkm2 *= big;
+		qkm1 *= big;
+		}
+	}
+while( ++n < 100 );
+
+cdone:
+return(ans);
+}
+
+
+/* power series */
+float incbpsf( float aa, float bb, float xx )
+{
+float a, b, x, t, u, y, s;
+
+a = aa;
+b = bb;
+x = xx;
+
+y = a * logf(x) + (b-1.0)*logf(1.0-x) - logf(a);
+y -= lgamf(a) + lgamf(b);
+y += lgamf(a+b);
+
+
+t = x / (1.0 - x);
+s = 0.0;
+u = 1.0;
+do
+	{
+	b -= 1.0;
+	if( b == 0.0 )
+		break;
+	a += 1.0;
+	u *= t*b/a;
+	s += u;
+	}
+while( fabsf(u) > MACHEPF );
+
+if( y < MINLOGF )
+	{
+	mtherr( "incbetf", UNDERFLOW );
+	s = 0.0;
+	}
+else
+	s = expf(y) * (1.0 + s);
+/*printf( "incbpsf: %.4e\n", s );*/
+return(s);
+}