uClibc now has a math library. muahahahaha!

 -Erik

1 files changed, 225 insertions, 0 deletions

diff --git a/libm/ldouble/stdtrl.c b/libm/ldouble/stdtrl.c
new file mode 100644
index 000000000..4218d4133
--- /dev/null
+++ b/libm/ldouble/stdtrl.c
@@ -0,0 +1,225 @@
+/*							stdtrl.c
+ *
+ *	Student's t distribution
+ *
+ *
+ *
+ * SYNOPSIS:
+ *
+ * long double p, t, stdtrl();
+ * int k;
+ *
+ * p = stdtrl( k, t );
+ *
+ *
+ * DESCRIPTION:
+ *
+ * Computes the integral from minus infinity to t of the Student
+ * t distribution with integer k > 0 degrees of freedom:
+ *
+ *                                      t
+ *                                      -
+ *                                     | |
+ *              -                      |         2   -(k+1)/2
+ *             | ( (k+1)/2 )           |  (     x   )
+ *       ----------------------        |  ( 1 + --- )        dx
+ *                     -               |  (      k  )
+ *       sqrt( k pi ) | ( k/2 )        |
+ *                                   | |
+ *                                    -
+ *                                   -inf.
+ * 
+ * Relation to incomplete beta integral:
+ *
+ *        1 - stdtr(k,t) = 0.5 * incbet( k/2, 1/2, z )
+ * where
+ *        z = k/(k + t**2).
+ *
+ * For t < -1.6, this is the method of computation.  For higher t,
+ * a direct method is derived from integration by parts.
+ * Since the function is symmetric about t=0, the area under the
+ * right tail of the density is found by calling the function
+ * with -t instead of t.
+ * 
+ * ACCURACY:
+ *
+ * Tested at random 1 <= k <= 100.  The "domain" refers to t.
+ *                      Relative error:
+ * arithmetic   domain     # trials      peak         rms
+ *    IEEE     -100,-1.6    10000       5.7e-18     9.8e-19
+ *    IEEE     -1.6,100     10000       3.8e-18     1.0e-19
+ */
+
+/*							stdtril.c
+ *
+ *	Functional inverse of Student's t distribution
+ *
+ *
+ *
+ * SYNOPSIS:
+ *
+ * long double p, t, stdtril();
+ * int k;
+ *
+ * t = stdtril( k, p );
+ *
+ *
+ * DESCRIPTION:
+ *
+ * Given probability p, finds the argument t such that stdtrl(k,t)
+ * is equal to p.
+ * 
+ * ACCURACY:
+ *
+ * Tested at random 1 <= k <= 100.  The "domain" refers to p:
+ *                      Relative error:
+ * arithmetic   domain     # trials      peak         rms
+ *    IEEE       0,1        3500       4.2e-17     4.1e-18
+ */
+
+
+/*
+Cephes Math Library Release 2.3:  January, 1995
+Copyright 1984, 1995 by Stephen L. Moshier
+*/
+
+#include <math.h>
+
+extern long double PIL, MACHEPL, MAXNUML;
+#ifdef ANSIPROT
+extern long double sqrtl ( long double );
+extern long double atanl ( long double );
+extern long double incbetl ( long double, long double, long double );
+extern long double incbil ( long double, long double, long double );
+extern long double fabsl ( long double );
+#else
+long double sqrtl(), atanl(), incbetl(), incbil(), fabsl();
+#endif
+
+long double stdtrl( k, t )
+int k;
+long double t;
+{
+long double x, rk, z, f, tz, p, xsqk;
+int j;
+
+if( k <= 0 )
+	{
+	mtherr( "stdtrl", DOMAIN );
+	return(0.0L);
+	}
+
+if( t == 0.0L )
+	return( 0.5L );
+
+if( t < -1.6L )
+	{
+	rk = k;
+	z = rk / (rk + t * t);
+	p = 0.5L * incbetl( 0.5L*rk, 0.5L, z );
+	return( p );
+	}
+
+/*	compute integral from -t to + t */
+
+if( t < 0.0L )
+	x = -t;
+else
+	x = t;
+
+rk = k;	/* degrees of freedom */
+z = 1.0L + ( x * x )/rk;
+
+/* test if k is odd or even */
+if( (k & 1) != 0)
+	{
+
+	/*	computation for odd k	*/
+
+	xsqk = x/sqrtl(rk);
+	p = atanl( xsqk );
+	if( k > 1 )
+		{
+		f = 1.0L;
+		tz = 1.0L;
+		j = 3;
+		while(  (j<=(k-2)) && ( (tz/f) > MACHEPL )  )
+			{
+			tz *= (j-1)/( z * j );
+			f += tz;
+			j += 2;
+			}
+		p += f * xsqk/z;
+		}
+	p *= 2.0L/PIL;
+	}
+
+
+else
+	{
+
+	/*	computation for even k	*/
+
+	f = 1.0L;
+	tz = 1.0L;
+	j = 2;
+
+	while(  ( j <= (k-2) ) && ( (tz/f) > MACHEPL )  )
+		{
+		tz *= (j - 1)/( z * j );
+		f += tz;
+		j += 2;
+		}
+	p = f * x/sqrtl(z*rk);
+	}
+
+/*	common exit	*/
+
+
+if( t < 0.0L )
+	p = -p;	/* note destruction of relative accuracy */
+
+	p = 0.5L + 0.5L * p;
+return(p);
+}
+
+
+long double stdtril( k, p )
+int k;
+long double p;
+{
+long double t, rk, z;
+int rflg;
+
+if( k <= 0 || p <= 0.0L || p >= 1.0L )
+	{
+	mtherr( "stdtril", DOMAIN );
+	return(0.0L);
+	}
+
+rk = k;
+
+if( p > 0.25L && p < 0.75L )
+	{
+	if( p == 0.5L )
+		return( 0.0L );
+	z = 1.0L - 2.0L * p;
+	z = incbil( 0.5L, 0.5L*rk, fabsl(z) );
+	t = sqrtl( rk*z/(1.0L-z) );
+	if( p < 0.5L )
+		t = -t;
+	return( t );
+	}
+rflg = -1;
+if( p >= 0.5L)
+	{
+	p = 1.0L - p;
+	rflg = 1;
+	}
+z = incbil( 0.5L*rk, 0.5L, 2.0L*p );
+
+if( MAXNUML * z < rk )
+	return(rflg* MAXNUML);
+t = sqrtl( rk/z - rk );
+return( rflg * t );
+}