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/double/incbet.c | 409 --------------------------------------------------- 1 file changed, 409 deletions(-) delete mode 100644 libm/double/incbet.c (limited to 'libm/double/incbet.c') 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 - -#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); -} -- cgit v1.2.3