From 7ce331c01ce6eb7b3f5c715a38a24359da9c6ee2 Mon Sep 17 00:00:00 2001 From: Eric Andersen Date: Thu, 22 Nov 2001 14:04:29 +0000 Subject: 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 --- libm/ldouble/bdtrl.c | 260 --------------------------------------------------- 1 file changed, 260 deletions(-) delete mode 100644 libm/ldouble/bdtrl.c (limited to 'libm/ldouble/bdtrl.c') diff --git a/libm/ldouble/bdtrl.c b/libm/ldouble/bdtrl.c deleted file mode 100644 index aca9577d1..000000000 --- a/libm/ldouble/bdtrl.c +++ /dev/null @@ -1,260 +0,0 @@ -/* 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 1 - */ - -/* bdtr() */ - - -/* -Cephes Math Library Release 2.3: March, 1995 -Copyright 1984, 1995 by Stephen L. Moshier -*/ - -#include -#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 ); -} -- cgit v1.2.3