summaryrefslogtreecommitdiff
path: root/libm/ldouble/bdtrl.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/bdtrl.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/bdtrl.c')
-rw-r--r--libm/ldouble/bdtrl.c260
1 files changed, 0 insertions, 260 deletions
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<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 );
-}