diff options
Diffstat (limited to 'libm/ldouble/stdtrl.c')
| -rw-r--r-- | libm/ldouble/stdtrl.c | 225 |
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 ); +} |