summaryrefslogtreecommitdiff
path: root/libm/ldouble/nbdtrl.c
diff options
context:
space:
mode:
Diffstat (limited to 'libm/ldouble/nbdtrl.c')
-rw-r--r--libm/ldouble/nbdtrl.c197
1 files changed, 197 insertions, 0 deletions
diff --git a/libm/ldouble/nbdtrl.c b/libm/ldouble/nbdtrl.c
new file mode 100644
index 000000000..91593f544
--- /dev/null
+++ b/libm/ldouble/nbdtrl.c
@@ -0,0 +1,197 @@
+/* nbdtrl.c
+ *
+ * Negative binomial distribution
+ *
+ *
+ *
+ * SYNOPSIS:
+ *
+ * int k, n;
+ * long double p, y, nbdtrl();
+ *
+ * y = nbdtrl( k, n, p );
+ *
+ *
+ *
+ * DESCRIPTION:
+ *
+ * Returns the sum of the terms 0 through k of the negative
+ * binomial distribution:
+ *
+ * k
+ * -- ( n+j-1 ) n j
+ * > ( ) p (1-p)
+ * -- ( j )
+ * j=0
+ *
+ * In a sequence of Bernoulli trials, this is the probability
+ * that k or fewer failures precede the nth success.
+ *
+ * The terms are not computed individually; instead the incomplete
+ * beta integral is employed, according to the formula
+ *
+ * y = nbdtr( k, n, p ) = incbet( n, k+1, p ).
+ *
+ * The arguments must be positive, with p ranging from 0 to 1.
+ *
+ *
+ *
+ * ACCURACY:
+ *
+ * Tested at random points (k,n,p) with k and n between 1 and 10,000
+ * and p between 0 and 1.
+ *
+ * arithmetic domain # trials peak rms
+ * Absolute error:
+ * IEEE 0,10000 10000 9.8e-15 2.1e-16
+ *
+ */
+ /* nbdtrcl.c
+ *
+ * Complemented negative binomial distribution
+ *
+ *
+ *
+ * SYNOPSIS:
+ *
+ * int k, n;
+ * long double p, y, nbdtrcl();
+ *
+ * y = nbdtrcl( k, n, p );
+ *
+ *
+ *
+ * DESCRIPTION:
+ *
+ * Returns the sum of the terms k+1 to infinity of the negative
+ * binomial distribution:
+ *
+ * inf
+ * -- ( n+j-1 ) n j
+ * > ( ) p (1-p)
+ * -- ( j )
+ * j=k+1
+ *
+ * The terms are not computed individually; instead the incomplete
+ * beta integral is employed, according to the formula
+ *
+ * y = nbdtrc( k, n, p ) = incbet( k+1, n, 1-p ).
+ *
+ * The arguments must be positive, with p ranging from 0 to 1.
+ *
+ *
+ *
+ * ACCURACY:
+ *
+ * See incbetl.c.
+ *
+ */
+ /* nbdtril
+ *
+ * Functional inverse of negative binomial distribution
+ *
+ *
+ *
+ * SYNOPSIS:
+ *
+ * int k, n;
+ * long double p, y, nbdtril();
+ *
+ * p = nbdtril( k, n, y );
+ *
+ *
+ *
+ * DESCRIPTION:
+ *
+ * Finds the argument p such that nbdtr(k,n,p) is equal to y.
+ *
+ * ACCURACY:
+ *
+ * Tested at random points (a,b,y), with y between 0 and 1.
+ *
+ * a,b Relative error:
+ * arithmetic domain # trials peak rms
+ * IEEE 0,100
+ * See also incbil.c.
+ */
+
+/*
+Cephes Math Library Release 2.3: January,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 powl ( long double, long double );
+extern long double incbil ( long double, long double, long double );
+#else
+long double incbetl(), powl(), incbil();
+#endif
+
+long double nbdtrcl( 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 )
+ {
+domerr:
+ mtherr( "nbdtrl", DOMAIN );
+ return( 0.0L );
+ }
+dn = n;
+if( k == 0 )
+ return( 1.0L - powl( p, dn ) );
+
+dk = k+1;
+return( incbetl( dk, dn, 1.0L - p ) );
+}
+
+
+
+long double nbdtrl( 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 )
+ {
+domerr:
+ mtherr( "nbdtrl", DOMAIN );
+ return( 0.0L );
+ }
+dn = n;
+if( k == 0 )
+ return( powl( p, dn ) );
+
+dk = k+1;
+return( incbetl( dn, dk, p ) );
+}
+
+
+long double nbdtril( k, n, p )
+int k, n;
+long double p;
+{
+long double dk, dn, w;
+
+if( (p < 0.0L) || (p > 1.0L) )
+ goto domerr;
+if( k < 0 )
+ {
+domerr:
+ mtherr( "nbdtrl", DOMAIN );
+ return( 0.0L );
+ }
+dk = k+1;
+dn = n;
+w = incbil( dn, dk, p );
+return( w );
+}