/* incbetl.c
*
* Incomplete beta integral
*
*
* SYNOPSIS:
*
* long double a, b, x, y, incbetl();
*
* y = incbetl( 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 random points (a,b,x) with x between 0 and 1.
* arithmetic domain # trials peak rms
* IEEE 0,5 20000 4.5e-18 2.4e-19
* IEEE 0,100 100000 3.9e-17 1.0e-17
* Half-integer a, b:
* IEEE .5,10000 100000 3.9e-14 4.4e-15
* Outputs smaller than the IEEE gradual underflow threshold
* were excluded from these statistics.
*
* ERROR MESSAGES:
*
* message condition value returned
* incbetl domain x<0, x>1 0.0
*/
/*
Cephes Math Library, Release 2.3: January, 1995
Copyright 1984, 1995 by Stephen L. Moshier
*/
#include <math.h>
#define MAXGAML 1755.455L
static long double big = 9.223372036854775808e18L;
static long double biginv = 1.084202172485504434007e-19L;
extern long double MACHEPL, MINLOGL, MAXLOGL;
#ifdef ANSIPROT
extern long double gammal ( long double );
extern long double lgaml ( long double );
extern long double expl ( long double );
extern long double logl ( long double );
extern long double fabsl ( long double );
extern long double powl ( long double, long double );
static long double incbcfl( long double, long double, long double );
static long double incbdl( long double, long double, long double );
static long double pseriesl( long double, long double, long double );
#else
long double gammal(), lgaml(), expl(), logl(), fabsl(), powl();
static long double incbcfl(), incbdl(), pseriesl();
#endif
long double incbetl( aa, bb, xx )
long double aa, bb, xx;
{
long double a, b, t, x, w, xc, y;
int flag;
if( aa <= 0.0L || bb <= 0.0L )
goto domerr;
if( (xx <= 0.0L) || ( xx >= 1.0L) )
{
if( xx == 0.0L )
return( 0.0L );
if( xx == 1.0L )
return( 1.0L );
domerr:
mtherr( "incbetl", DOMAIN );
return( 0.0L );
}
flag = 0;
if( (bb * xx) <= 1.0L && xx <= 0.95L)
{
t = pseriesl(aa, bb, xx);
goto done;
}
w = 1.0L - 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.0L && x <= 0.95L)
{
t = pseriesl(a, b, x);
goto done;
}
/* Choose expansion for optimal convergence */
y = x * (a+b-2.0L) - (a-1.0L);
if( y < 0.0L )
w = incbcfl( a, b, x );
else
w = incbdl( a, b, x ) / xc;
/* Multiply w by the factor
a b _ _ _
x (1-x) | (a+b) / ( a | (a) | (b) ) . */
y = a * logl(x);
t = b * logl(xc);
if( (a+b) < MAXGAML && fabsl(y) < MAXLOGL && fabsl(t) < MAXLOGL )
{
t = powl(xc,b);
t *= powl(x,a);
t /= a;
t *= w;
t *= gammal(a+b) / (gammal(a) * gammal(b));
goto done;
}
else
{
/* Resort to logarithms. */
y += t + lgaml(a+b) - lgaml(a) - lgaml(b);
y += logl(w/a);
if( y < MINLOGL )
t = 0.0L;
else
t = expl(y);
}
done:
if( flag == 1 )
{
if( t <= MACHEPL )
t = 1.0L - MACHEPL;
else
t = 1.0L - t;
}
return( t );
}
/* Continued fraction expansion #1
* for incomplete beta integral
*/
static long double incbcfl( a, b, x )
long double a, b, x;
{
long double xk, pk, pkm1, pkm2, qk, qkm1, qkm2;
long double k1, k2, k3, k4, k5, k6, k7, k8;
long double r, t, ans, thresh;
int n;
k1 = a;
k2 = a + b;
k3 = a;
k4 = a + 1.0L;
k5 = 1.0L;
k6 = b - 1.0L;
k7 = k4;
k8 = a + 2.0L;
pkm2 = 0.0L;
qkm2 = 1.0L;
pkm1 = 1.0L;
qkm1 = 1.0L;
ans = 1.0L;
r = 1.0L;
n = 0;
thresh = 3.0L * MACHEPL;
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.0L )
r = pk/qk;
if( r != 0.0L )
{
t = fabsl( (ans - r)/r );
ans = r;
}
else
t = 1.0L;
if( t < thresh )
goto cdone;
k1 += 1.0L;
k2 += 1.0L;
k3 += 2.0L;
k4 += 2.0L;
k5 += 1.0L;
k6 -= 1.0L;
k7 += 2.0L;
k8 += 2.0L;
if( (fabsl(qk) + fabsl(pk)) > big )
{
pkm2 *= biginv;
pkm1 *= biginv;
qkm2 *= biginv;
qkm1 *= biginv;
}
if( (fabsl(qk) < biginv) || (fabsl(pk) < biginv) )
{
pkm2 *= big;
pkm1 *= big;
qkm2 *= big;
qkm1 *= big;
}
}
while( ++n < 400 );
mtherr( "incbetl", PLOSS );
cdone:
return(ans);
}
/* Continued fraction expansion #2
* for incomplete beta integral
*/
static long double incbdl( a, b, x )
long double a, b, x;
{
long double xk, pk, pkm1, pkm2, qk, qkm1, qkm2;
long double k1, k2, k3, k4, k5, k6, k7, k8;
long double r, t, ans, z, thresh;
int n;
k1 = a;
k2 = b - 1.0L;
k3 = a;
k4 = a + 1.0L;
k5 = 1.0L;
k6 = a + b;
k7 = a + 1.0L;
k8 = a + 2.0L;
pkm2 = 0.0L;
qkm2 = 1.0L;
pkm1 = 1.0L;
qkm1 = 1.0L;
z = x / (1.0L-x);
ans = 1.0L;
r = 1.0L;
n = 0;
thresh = 3.0L * MACHEPL;
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.0L )
r = pk/qk;
if( r != 0.0L )
{
t = fabsl( (ans - r)/r );
ans = r;
}
else
t = 1.0L;
if( t < thresh )
goto cdone;
k1 += 1.0L;
k2 -= 1.0L;
k3 += 2.0L;
k4 += 2.0L;
k5 += 1.0L;
k6 += 1.0L;
k7 += 2.0L;
k8 += 2.0L;
if( (fabsl(qk) + fabsl(pk)) > big )
{
pkm2 *= biginv;
pkm1 *= biginv;
qkm2 *= biginv;
qkm1 *= biginv;
}
if( (fabsl(qk) < biginv) || (fabsl(pk) < biginv) )
{
pkm2 *= big;
pkm1 *= big;
qkm2 *= big;
qkm1 *= big;
}
}
while( ++n < 400 );
mtherr( "incbetl", PLOSS );
cdone:
return(ans);
}
/* Power series for incomplete gamma integral.
Use when b*x is small. */
static long double pseriesl( a, b, x )
long double a, b, x;
{
long double s, t, u, v, n, t1, z, ai;
ai = 1.0L / a;
u = (1.0L - b) * x;
v = u / (a + 1.0L);
t1 = v;
t = u;
n = 2.0L;
s = 0.0L;
z = MACHEPL * ai;
while( fabsl(v) > z )
{
u = (n - b) * x / n;
t *= u;
v = t / (a + n);
s += v;
n += 1.0L;
}
s += t1;
s += ai;
u = a * logl(x);
if( (a+b) < MAXGAML && fabsl(u) < MAXLOGL )
{
t = gammal(a+b)/(gammal(a)*gammal(b));
s = s * t * powl(x,a);
}
else
{
t = lgaml(a+b) - lgaml(a) - lgaml(b) + u + logl(s);
if( t < MINLOGL )
s = 0.0L;
else
s = expl(t);
}
return(s);
}