summaryrefslogtreecommitdiff
path: root/libm/ldouble/stdtrl.c
diff options
context:
space:
mode:
authorEric Andersen <andersen@codepoet.org>2001-11-22 14:04:29 +0000
committerEric Andersen <andersen@codepoet.org>2001-11-22 14:04:29 +0000
commit7ce331c01ce6eb7b3f5c715a38a24359da9c6ee2 (patch)
tree3a7e8476e868ae15f4da1b7ce26b2db6f434468c /libm/ldouble/stdtrl.c
parentc117dd5fb183afb1a4790a6f6110d88704be6bf8 (diff)
Totally rework the math library, this time based on the MacOs X
math library (which is itself based on the math lib from FreeBSD). -Erik
Diffstat (limited to 'libm/ldouble/stdtrl.c')
-rw-r--r--libm/ldouble/stdtrl.c225
1 files changed, 0 insertions, 225 deletions
diff --git a/libm/ldouble/stdtrl.c b/libm/ldouble/stdtrl.c
deleted file mode 100644
index 4218d4133..000000000
--- a/libm/ldouble/stdtrl.c
+++ /dev/null
@@ -1,225 +0,0 @@
-/* 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 );
-}