summaryrefslogtreecommitdiff
path: root/libm/ldouble/stdtrl.c
diff options
context:
space:
mode:
authorEric Andersen <andersen@codepoet.org>2001-05-10 00:40:28 +0000
committerEric Andersen <andersen@codepoet.org>2001-05-10 00:40:28 +0000
commit1077fa4d772832f77a677ce7fb7c2d513b959e3f (patch)
tree579bee13fb0b58d2800206366ec2caecbb15f3fc /libm/ldouble/stdtrl.c
parent22358dd7ce7bb49792204b698f01a6f69b9c8e08 (diff)
uClibc now has a math library. muahahahaha!
-Erik
Diffstat (limited to 'libm/ldouble/stdtrl.c')
-rw-r--r--libm/ldouble/stdtrl.c225
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 );
+}