diff options
Diffstat (limited to 'libm/ldouble/bdtrl.c')
| -rw-r--r-- | libm/ldouble/bdtrl.c | 260 |
1 files changed, 260 insertions, 0 deletions
diff --git a/libm/ldouble/bdtrl.c b/libm/ldouble/bdtrl.c new file mode 100644 index 000000000..aca9577d1 --- /dev/null +++ b/libm/ldouble/bdtrl.c @@ -0,0 +1,260 @@ +/* bdtrl.c + * + * Binomial distribution + * + * + * + * SYNOPSIS: + * + * int k, n; + * long double p, y, bdtrl(); + * + * y = bdtrl( k, n, p ); + * + * + * + * DESCRIPTION: + * + * Returns the sum of the terms 0 through k of the Binomial + * probability density: + * + * k + * -- ( n ) j n-j + * > ( ) p (1-p) + * -- ( j ) + * j=0 + * + * The terms are not summed directly; instead the incomplete + * beta integral is employed, according to the formula + * + * y = bdtr( k, n, p ) = incbet( n-k, k+1, 1-p ). + * + * The arguments must be positive, with p ranging from 0 to 1. + * + * + * + * ACCURACY: + * + * Tested at random points (k,n,p) with a and b between 0 + * and 10000 and p between 0 and 1. + * Relative error: + * arithmetic domain # trials peak rms + * IEEE 0,10000 3000 1.6e-14 2.2e-15 + * + * ERROR MESSAGES: + * + * message condition value returned + * bdtrl domain k < 0 0.0 + * n < k + * x < 0, x > 1 + * + */ + /* bdtrcl() + * + * Complemented binomial distribution + * + * + * + * SYNOPSIS: + * + * int k, n; + * long double p, y, bdtrcl(); + * + * y = bdtrcl( k, n, p ); + * + * + * + * DESCRIPTION: + * + * Returns the sum of the terms k+1 through n of the Binomial + * probability density: + * + * n + * -- ( n ) j n-j + * > ( ) p (1-p) + * -- ( j ) + * j=k+1 + * + * The terms are not summed directly; instead the incomplete + * beta integral is employed, according to the formula + * + * y = bdtrc( k, n, p ) = incbet( k+1, n-k, p ). + * + * The arguments must be positive, with p ranging from 0 to 1. + * + * + * + * ACCURACY: + * + * See incbet.c. + * + * ERROR MESSAGES: + * + * message condition value returned + * bdtrcl domain x<0, x>1, n<k 0.0 + */ + /* bdtril() + * + * Inverse binomial distribution + * + * + * + * SYNOPSIS: + * + * int k, n; + * long double p, y, bdtril(); + * + * p = bdtril( k, n, y ); + * + * + * + * DESCRIPTION: + * + * Finds the event probability p such that the sum of the + * terms 0 through k of the Binomial probability density + * is equal to the given cumulative probability y. + * + * This is accomplished using the inverse beta integral + * function and the relation + * + * 1 - p = incbi( n-k, k+1, y ). + * + * ACCURACY: + * + * See incbi.c. + * Tested at random k, n between 1 and 10000. The "domain" refers to p: + * Relative error: + * arithmetic domain # trials peak rms + * IEEE 0,1 3500 2.0e-15 8.2e-17 + * + * ERROR MESSAGES: + * + * message condition value returned + * bdtril domain k < 0, n <= k 0.0 + * x < 0, x > 1 + */ + +/* bdtr() */ + + +/* +Cephes Math Library Release 2.3: March, 1995 +Copyright 1984, 1995 by Stephen L. Moshier +*/ + +#include <math.h> +#ifdef ANSIPROT +extern long double incbetl ( long double, long double, long double ); +extern long double incbil ( long double, long double, long double ); +extern long double powl ( long double, long double ); +extern long double expm1l ( long double ); +extern long double log1pl ( long double ); +#else +long double incbetl(), incbil(), powl(), expm1l(), log1pl(); +#endif + +long double bdtrcl( k, n, p ) +int k, n; +long double p; +{ +long double dk, dn; + +if( (p < 0.0L) || (p > 1.0L) ) + goto domerr; +if( k < 0 ) + return( 1.0L ); + +if( n < k ) + { +domerr: + mtherr( "bdtrcl", DOMAIN ); + return( 0.0L ); + } + +if( k == n ) + return( 0.0L ); +dn = n - k; +if( k == 0 ) + { + if( p < .01L ) + dk = -expm1l( dn * log1pl(-p) ); + else + dk = 1.0L - powl( 1.0L-p, dn ); + } +else + { + dk = k + 1; + dk = incbetl( dk, dn, p ); + } +return( dk ); +} + + + +long double bdtrl( k, n, p ) +int k, n; +long double p; +{ +long double dk, dn, q; + +if( (p < 0.0L) || (p > 1.0L) ) + goto domerr; +if( (k < 0) || (n < k) ) + { +domerr: + mtherr( "bdtrl", DOMAIN ); + return( 0.0L ); + } + +if( k == n ) + return( 1.0L ); + +q = 1.0L - p; +dn = n - k; +if( k == 0 ) + { + dk = powl( q, dn ); + } +else + { + dk = k + 1; + dk = incbetl( dn, dk, q ); + } +return( dk ); +} + + +long double bdtril( k, n, y ) +int k, n; +long double y; +{ +long double dk, dn, p; + +if( (y < 0.0L) || (y > 1.0L) ) + goto domerr; +if( (k < 0) || (n <= k) ) + { +domerr: + mtherr( "bdtril", DOMAIN ); + return( 0.0L ); + } + +dn = n - k; +if( k == 0 ) + { + if( y > 0.8L ) + p = -expm1l( log1pl(y-1.0L) / dn ); + else + p = 1.0L - powl( y, 1.0L/dn ); + } +else + { + dk = k + 1; + p = incbetl( dn, dk, y ); + if( p > 0.5 ) + p = incbil( dk, dn, 1.0L-y ); + else + p = 1.0 - incbil( dn, dk, y ); + } +return( p ); +} |