Category : C Source Code
Archive   : CEPHES22.ZIP
Filename : INCBET.C

 
Output of file : INCBET.C contained in archive : CEPHES22.ZIP
/* 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.
* If a < 1, the function calls itself recursively after a
* transformation to increase a to a+1.
*
* ACCURACY:
*
* Tested at random points (a,b,x) with a and b between 0
* and 100 and x between 0 and 1.
* Relative error (x ranges from 0 to 1):
* arithmetic domain # trials peak rms
* DEC 0,100 3300 3.5e-14 5.0e-15
* IEEE 0,100 10000 3.9e-13 5.2e-14
*
* Larger errors may occur for extreme ratios of a and b.
*
* ERROR MESSAGES:
* message condition value returned
* incbet domain x<0, x>1 0.0
*/


/*
Cephes Math Library, Release 2.2: July, 1992
Copyright 1984, 1987 by Stephen L. Moshier
Direct inquiries to 30 Frost Street, Cambridge, MA 02140
*/

#include "mconf.h"

#define BIG 1.44115188075855872E+17
extern double MACHEP, MINLOG, MAXLOG;

double incbet( aa, bb, xx )
double aa, bb, xx;
{
double ans, a, b, t, x, onemx;
double lgam(), exp(), log(), fabs();
double incbd(), incbcf();
short flag;

if( (xx <= 0.0) || ( xx >= 1.0) )
{
if( xx == 0.0 )
return(0.0);
if( xx == 1.0 )
return( 1.0 );
mtherr( "incbet", DOMAIN );
return( 0.0 );
}

onemx = 1.0 - xx;

/* transformation for small aa */
if( aa <= 1.0 )
{
ans = incbet( aa+1.0, bb, xx );
t = aa*log(xx) + bb*log( 1.0-xx )
+ lgam(aa+bb) - lgam(aa+1.0) - lgam(bb);
if( t > MINLOG )
ans += exp(t);
return( ans );
}

/* see if x is greater than the mean */

if( xx > (aa/(aa+bb)) )
{
flag = 1;
a = bb;
b = aa;
t = xx;
x = onemx;
}
else
{
flag = 0;
a = aa;
b = bb;
t = onemx;
x = xx;
}

/* Choose expansion for optimal convergence */

ans = x * (a+b-2.0)/(a-1.0);
if( ans < 1.0 )
{
ans = incbcf( a, b, x );
t = b * log( t );
}
else
{
ans = incbd( a, b, x );
t = (b-1.0) * log(t);
}

adone:
t += a*log(x) + lgam(a+b) - lgam(a) - lgam(b);
t += log( ans/a );

if( t < MINLOG )
{
if( flag == 0 )
{
mtherr( "incbet", UNDERFLOW );
return( 0.0 );
}
else
return(1.0);
}

t = exp(t);

if( flag == 1 )
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;
static double big = BIG;
double fabs();
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;

n = 0;
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 < MACHEP )
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 /= big;
pkm1 /= big;
qkm2 /= big;
qkm1 /= big;
}
if( (fabs(qk) < MACHEP) || (fabs(pk) < MACHEP) )
{
pkm2 *= big;
pkm1 *= big;
qkm2 *= big;
qkm1 *= big;
}
}
while( ++n < 100 );

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;
static double big = BIG;
double fabs();
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;

n = 0;
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 < MACHEP )
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 /= big;
pkm1 /= big;
qkm2 /= big;
qkm1 /= big;
}
if( (fabs(qk) < MACHEP) || (fabs(pk) < MACHEP) )
{
pkm2 *= big;
pkm1 *= big;
qkm2 *= big;
qkm1 *= big;
}
}
while( ++n < 100 );

cdone:
return(ans);
}


  3 Responses to “Category : C Source Code
Archive   : CEPHES22.ZIP
Filename : INCBET.C

  1. Very nice! Thank you for this wonderful archive. I wonder why I found it only now. Long live the BBS file archives!

  2. This is so awesome! 😀 I’d be cool if you could download an entire archive of this at once, though.

  3. But one thing that puzzles me is the “mtswslnkmcjklsdlsbdmMICROSOFT” string. There is an article about it here. It is definitely worth a read: http://www.os2museum.com/wp/mtswslnk/