#include "cdflib.h"
void cumf(double *f,double *dfn,double *dfd,double *cum,double *ccum)
/*
**********************************************************************
void cumf(double *f,double *dfn,double *dfd,double *cum,double *ccum)
CUMulative F distribution
Function
Computes the integral from 0 to F of the f-density with DFN
and DFD degrees of freedom.
Arguments
F --> Upper limit of integration of the f-density.
F is DOUBLE PRECISION
DFN --> Degrees of freedom of the numerator sum of squares.
DFN is DOUBLE PRECISI
DFD --> Degrees of freedom of the denominator sum of squares.
DFD is DOUBLE PRECISI
CUM <-- Cumulative f distribution.
CUM is DOUBLE PRECISI
CCUM <-- Compliment of Cumulative f distribution.
CCUM is DOUBLE PRECIS
Method
Formula 26.5.28 of Abramowitz and Stegun is used to reduce
the cumulative F to a cumulative beta distribution.
Note
If F is less than or equal to 0, 0 is returned.
**********************************************************************
*/
{
#define half 0.5e0
#define done 1.0e0
static double dsum,prod,xx,yy;
static int ierr;
static double T1,T2;
/*
..
.. Executable Statements ..
*/
if(!(*f <= 0.0e0)) goto S10;
*cum = 0.0e0;
*ccum = 1.0e0;
return;
S10:
prod = *dfn**f;
/*
XX is such that the incomplete beta with parameters
DFD/2 and DFN/2 evaluated at XX is 1 - CUM or CCUM
YY is 1 - XX
Calculate the smaller of XX and YY accurately
*/
dsum = *dfd+prod;
xx = *dfd/dsum;
if(xx > half) {
yy = prod/dsum;
xx = done-yy;
}
else yy = done-xx;
T1 = *dfd*half;
T2 = *dfn*half;
bratio(&T1,&T2,&xx,&yy,ccum,cum,&ierr);
return;
#undef half
#undef done
} /* END */