summaryrefslogtreecommitdiff
path: root/libm/double/incbet.c
diff options
context:
space:
mode:
Diffstat (limited to 'libm/double/incbet.c')
-rw-r--r--libm/double/incbet.c409
1 files changed, 0 insertions, 409 deletions
diff --git a/libm/double/incbet.c b/libm/double/incbet.c
deleted file mode 100644
index ec236747d..000000000
--- a/libm/double/incbet.c
+++ /dev/null
@@ -1,409 +0,0 @@
-/* incbet.c
- *
- * Incomplete beta integral
- *
- *
- * SYNOPSIS:
- *
- * double a, b, x, y, incbet();
- *
- * y = incbet( a, b, x );
- *
- *
- * DESCRIPTION:
- *
- * Returns incomplete beta integral of the arguments, evaluated
- * from zero to x. The function is defined as
- *
- * x
- * - -
- * | (a+b) | | a-1 b-1
- * ----------- | t (1-t) dt.
- * - - | |
- * | (a) | (b) -
- * 0
- *
- * The domain of definition is 0 <= x <= 1. In this
- * implementation a and b are restricted to positive values.
- * The integral from x to 1 may be obtained by the symmetry
- * relation
- *
- * 1 - incbet( a, b, x ) = incbet( b, a, 1-x ).
- *
- * The integral is evaluated by a continued fraction expansion
- * or, when b*x is small, by a power series.
- *
- * ACCURACY:
- *
- * Tested at uniformly distributed random points (a,b,x) with a and b
- * in "domain" and x between 0 and 1.
- * Relative error
- * arithmetic domain # trials peak rms
- * IEEE 0,5 10000 6.9e-15 4.5e-16
- * IEEE 0,85 250000 2.2e-13 1.7e-14
- * IEEE 0,1000 30000 5.3e-12 6.3e-13
- * IEEE 0,10000 250000 9.3e-11 7.1e-12
- * IEEE 0,100000 10000 8.7e-10 4.8e-11
- * Outputs smaller than the IEEE gradual underflow threshold
- * were excluded from these statistics.
- *
- * ERROR MESSAGES:
- * message condition value returned
- * incbet domain x<0, x>1 0.0
- * incbet underflow 0.0
- */
-
-
-/*
-Cephes Math Library, Release 2.8: June, 2000
-Copyright 1984, 1995, 2000 by Stephen L. Moshier
-*/
-
-#include <math.h>
-
-#ifdef DEC
-#define MAXGAM 34.84425627277176174
-#else
-#define MAXGAM 171.624376956302725
-#endif
-
-extern double MACHEP, MINLOG, MAXLOG;
-#ifdef ANSIPROT
-extern double gamma ( double );
-extern double lgam ( double );
-extern double exp ( double );
-extern double log ( double );
-extern double pow ( double, double );
-extern double fabs ( double );
-static double incbcf(double, double, double);
-static double incbd(double, double, double);
-static double pseries(double, double, double);
-#else
-double gamma(), lgam(), exp(), log(), pow(), fabs();
-static double incbcf(), incbd(), pseries();
-#endif
-
-static double big = 4.503599627370496e15;
-static double biginv = 2.22044604925031308085e-16;
-
-
-double incbet( aa, bb, xx )
-double aa, bb, xx;
-{
-double a, b, t, x, xc, w, y;
-int flag;
-
-if( aa <= 0.0 || bb <= 0.0 )
- goto domerr;
-
-if( (xx <= 0.0) || ( xx >= 1.0) )
- {
- if( xx == 0.0 )
- return(0.0);
- if( xx == 1.0 )
- return( 1.0 );
-domerr:
- mtherr( "incbet", DOMAIN );
- return( 0.0 );
- }
-
-flag = 0;
-if( (bb * xx) <= 1.0 && xx <= 0.95)
- {
- t = pseries(aa, bb, xx);
- goto done;
- }
-
-w = 1.0 - xx;
-
-/* Reverse a and b if x is greater than the mean. */
-if( xx > (aa/(aa+bb)) )
- {
- flag = 1;
- a = bb;
- b = aa;
- xc = xx;
- x = w;
- }
-else
- {
- a = aa;
- b = bb;
- xc = w;
- x = xx;
- }
-
-if( flag == 1 && (b * x) <= 1.0 && x <= 0.95)
- {
- t = pseries(a, b, x);
- goto done;
- }
-
-/* Choose expansion for better convergence. */
-y = x * (a+b-2.0) - (a-1.0);
-if( y < 0.0 )
- w = incbcf( a, b, x );
-else
- w = incbd( a, b, x ) / xc;
-
-/* Multiply w by the factor
- a b _ _ _
- x (1-x) | (a+b) / ( a | (a) | (b) ) . */
-
-y = a * log(x);
-t = b * log(xc);
-if( (a+b) < MAXGAM && fabs(y) < MAXLOG && fabs(t) < MAXLOG )
- {
- t = pow(xc,b);
- t *= pow(x,a);
- t /= a;
- t *= w;
- t *= gamma(a+b) / (gamma(a) * gamma(b));
- goto done;
- }
-/* Resort to logarithms. */
-y += t + lgam(a+b) - lgam(a) - lgam(b);
-y += log(w/a);
-if( y < MINLOG )
- t = 0.0;
-else
- t = exp(y);
-
-done:
-
-if( flag == 1 )
- {
- if( t <= MACHEP )
- t = 1.0 - MACHEP;
- else
- t = 1.0 - t;
- }
-return( t );
-}
-
-/* Continued fraction expansion #1
- * for incomplete beta integral
- */
-
-static double incbcf( a, b, x )
-double a, b, x;
-{
-double xk, pk, pkm1, pkm2, qk, qkm1, qkm2;
-double k1, k2, k3, k4, k5, k6, k7, k8;
-double r, t, ans, thresh;
-int n;
-
-k1 = a;
-k2 = a + b;
-k3 = a;
-k4 = a + 1.0;
-k5 = 1.0;
-k6 = b - 1.0;
-k7 = k4;
-k8 = a + 2.0;
-
-pkm2 = 0.0;
-qkm2 = 1.0;
-pkm1 = 1.0;
-qkm1 = 1.0;
-ans = 1.0;
-r = 1.0;
-n = 0;
-thresh = 3.0 * MACHEP;
-do
- {
-
- xk = -( x * k1 * k2 )/( k3 * k4 );
- pk = pkm1 + pkm2 * xk;
- qk = qkm1 + qkm2 * xk;
- pkm2 = pkm1;
- pkm1 = pk;
- qkm2 = qkm1;
- qkm1 = qk;
-
- xk = ( x * k5 * k6 )/( k7 * k8 );
- pk = pkm1 + pkm2 * xk;
- qk = qkm1 + qkm2 * xk;
- pkm2 = pkm1;
- pkm1 = pk;
- qkm2 = qkm1;
- qkm1 = qk;
-
- if( qk != 0 )
- r = pk/qk;
- if( r != 0 )
- {
- t = fabs( (ans - r)/r );
- ans = r;
- }
- else
- t = 1.0;
-
- if( t < thresh )
- goto cdone;
-
- k1 += 1.0;
- k2 += 1.0;
- k3 += 2.0;
- k4 += 2.0;
- k5 += 1.0;
- k6 -= 1.0;
- k7 += 2.0;
- k8 += 2.0;
-
- if( (fabs(qk) + fabs(pk)) > big )
- {
- pkm2 *= biginv;
- pkm1 *= biginv;
- qkm2 *= biginv;
- qkm1 *= biginv;
- }
- if( (fabs(qk) < biginv) || (fabs(pk) < biginv) )
- {
- pkm2 *= big;
- pkm1 *= big;
- qkm2 *= big;
- qkm1 *= big;
- }
- }
-while( ++n < 300 );
-
-cdone:
-return(ans);
-}
-
-
-/* Continued fraction expansion #2
- * for incomplete beta integral
- */
-
-static double incbd( a, b, x )
-double a, b, x;
-{
-double xk, pk, pkm1, pkm2, qk, qkm1, qkm2;
-double k1, k2, k3, k4, k5, k6, k7, k8;
-double r, t, ans, z, thresh;
-int n;
-
-k1 = a;
-k2 = b - 1.0;
-k3 = a;
-k4 = a + 1.0;
-k5 = 1.0;
-k6 = a + b;
-k7 = a + 1.0;;
-k8 = a + 2.0;
-
-pkm2 = 0.0;
-qkm2 = 1.0;
-pkm1 = 1.0;
-qkm1 = 1.0;
-z = x / (1.0-x);
-ans = 1.0;
-r = 1.0;
-n = 0;
-thresh = 3.0 * MACHEP;
-do
- {
-
- xk = -( z * k1 * k2 )/( k3 * k4 );
- pk = pkm1 + pkm2 * xk;
- qk = qkm1 + qkm2 * xk;
- pkm2 = pkm1;
- pkm1 = pk;
- qkm2 = qkm1;
- qkm1 = qk;
-
- xk = ( z * k5 * k6 )/( k7 * k8 );
- pk = pkm1 + pkm2 * xk;
- qk = qkm1 + qkm2 * xk;
- pkm2 = pkm1;
- pkm1 = pk;
- qkm2 = qkm1;
- qkm1 = qk;
-
- if( qk != 0 )
- r = pk/qk;
- if( r != 0 )
- {
- t = fabs( (ans - r)/r );
- ans = r;
- }
- else
- t = 1.0;
-
- if( t < thresh )
- goto cdone;
-
- k1 += 1.0;
- k2 -= 1.0;
- k3 += 2.0;
- k4 += 2.0;
- k5 += 1.0;
- k6 += 1.0;
- k7 += 2.0;
- k8 += 2.0;
-
- if( (fabs(qk) + fabs(pk)) > big )
- {
- pkm2 *= biginv;
- pkm1 *= biginv;
- qkm2 *= biginv;
- qkm1 *= biginv;
- }
- if( (fabs(qk) < biginv) || (fabs(pk) < biginv) )
- {
- pkm2 *= big;
- pkm1 *= big;
- qkm2 *= big;
- qkm1 *= big;
- }
- }
-while( ++n < 300 );
-cdone:
-return(ans);
-}
-
-/* Power series for incomplete beta integral.
- Use when b*x is small and x not too close to 1. */
-
-static double pseries( a, b, x )
-double a, b, x;
-{
-double s, t, u, v, n, t1, z, ai;
-
-ai = 1.0 / a;
-u = (1.0 - b) * x;
-v = u / (a + 1.0);
-t1 = v;
-t = u;
-n = 2.0;
-s = 0.0;
-z = MACHEP * ai;
-while( fabs(v) > z )
- {
- u = (n - b) * x / n;
- t *= u;
- v = t / (a + n);
- s += v;
- n += 1.0;
- }
-s += t1;
-s += ai;
-
-u = a * log(x);
-if( (a+b) < MAXGAM && fabs(u) < MAXLOG )
- {
- t = gamma(a+b)/(gamma(a)*gamma(b));
- s = s * t * pow(x,a);
- }
-else
- {
- t = lgam(a+b) - lgam(a) - lgam(b) + u + log(s);
- if( t < MINLOG )
- s = 0.0;
- else
- s = exp(t);
- }
-return(s);
-}
generated by cgit v1.2.3 (git 2.25.1) at 2024-11-05 02:37:46 +0000