#include "cdflib.h"
void cumchn(double *x,double *df,double *pnonc,double *cum,
	    double *ccum)
/*
**********************************************************************
 
     void cumchn(double *x,double *df,double *pnonc,double *cum,
            double *ccum)

             CUMulative of the Non-central CHi-square distribution
 
 
                              Function
 
 
     Calculates     the       cumulative      non-central    chi-square
     distribution, i.e.,  the probability   that  a   random   variable
     which    follows  the  non-central chi-square  distribution,  with
     non-centrality  parameter    PNONC  and   continuous  degrees   of
     freedom DF, is less than or equal to X.
 
 
                              Arguments
 
 
     X       --> Upper limit of integration of the non-central
                 chi-square distribution.
                                                 X is DOUBLE PRECISION
 
     DF      --> Degrees of freedom of the non-central
                 chi-square distribution.
                                                 DF is DOUBLE PRECISION
 
     PNONC   --> Non-centrality parameter of the non-central
                 chi-square distribution.
                                                 PNONC is DOUBLE PRECIS
 
     CUM <-- Cumulative non-central chi-square distribution.
                                                 CUM is DOUBLE PRECISIO
 
     CCUM <-- Compliment of Cumulative non-central chi-square distribut
                                                 CCUM is DOUBLE PRECISI
 
 
                              Method
 
 
     Uses  formula  26.4.25   of  Abramowitz  and  Stegun, Handbook  of
     Mathematical    Functions,  US   NBS   (1966)    to calculate  the
     non-central chi-square.
 
 
                              Variables
 
 
     EPS     --- Convergence criterion.  The sum stops when a
                 term is less than EPS*SUM.
                                                 EPS is DOUBLE PRECISIO
 
     NTIRED  --- Maximum number of terms to be evaluated
                 in each sum.
                                                 NTIRED is INTEGER
 
     QCONV   --- .TRUE. if convergence achieved -
                 i.e., program did not stop on NTIRED criterion.
                                                 QCONV is LOGICAL
 
     CCUM <-- Compliment of Cumulative non-central
              chi-square distribution.
                                                 CCUM is DOUBLE PRECISI
 
**********************************************************************
*/
{
#define dg(i) (*df+2.0e0*(double)(i))
#define qsmall(xx) (int)(sum < 1.0e-20 || (xx) < eps*sum)
#define qtired(i) (int)((i) > ntired)
static double eps = 1.0e-5;
static int ntired = 1000;
static double adj,centaj,centwt,chid2,dfd2,lcntaj,lcntwt,lfact,pcent,pterm,sum,
    sumadj,term,wt,xnonc;
static int i,icent,iterb,iterf;
static double T1,T2,T3;
/*
     ..
     .. Executable Statements ..
*/
    if(!(*x <= 0.0e0)) goto S10;
    *cum = 0.0e0;
    *ccum = 1.0e0;
    return;
S10:
    if(!(*pnonc <= 1.0e-10)) goto S20;
/*
     When non-centrality parameter is (essentially) zero,
     use cumulative chi-square distribution
*/
    cumchi(x,df,cum,ccum);
    return;
S20:
    xnonc = *pnonc/2.0e0;
/*
**********************************************************************
     The following code calcualtes the weight, chi-square, and
     adjustment term for the central term in the infinite series.
     The central term is the one in which the poisson weight is
     greatest.  The adjustment term is the amount that must
     be subtracted from the chi-square to move up two degrees
     of freedom.
**********************************************************************
*/
    icent = fifidint(xnonc);
    if(icent == 0) icent = 1;
    chid2 = *x/2.0e0;
/*
     Calculate central weight term
*/
    T1 = (double)(icent+1);
    lfact = alngam(&T1);
    lcntwt = -xnonc+(double)icent*log(xnonc)-lfact;
    centwt = exp(lcntwt);
/*
     Calculate central chi-square
*/
    T2 = dg(icent);
    cumchi(x,&T2,&pcent,ccum);
/*
     Calculate central adjustment term
*/
    dfd2 = dg(icent)/2.0e0;
    T3 = 1.0e0+dfd2;
    lfact = alngam(&T3);
    lcntaj = dfd2*log(chid2)-chid2-lfact;
    centaj = exp(lcntaj);
    sum = centwt*pcent;
/*
**********************************************************************
     Sum backwards from the central term towards zero.
     Quit whenever either
     (1) the zero term is reached, or
     (2) the term gets small relative to the sum, or
     (3) More than NTIRED terms are totaled.
**********************************************************************
*/
    iterb = 0;
    sumadj = 0.0e0;
    adj = centaj;
    wt = centwt;
    i = icent;
    goto S40;
S30:
    if(qtired(iterb) || qsmall(term) || i == 0) goto S50;
S40:
    dfd2 = dg(i)/2.0e0;
/*
     Adjust chi-square for two fewer degrees of freedom.
     The adjusted value ends up in PTERM.
*/
    adj = adj*dfd2/chid2;
    sumadj += adj;
    pterm = pcent+sumadj;
/*
     Adjust poisson weight for J decreased by one
*/
    wt *= ((double)i/xnonc);
    term = wt*pterm;
    sum += term;
    i -= 1;
    iterb += 1;
    goto S30;
S50:
    iterf = 0;
/*
**********************************************************************
     Now sum forward from the central term towards infinity.
     Quit when either
     (1) the term gets small relative to the sum, or
     (2) More than NTIRED terms are totaled.
**********************************************************************
*/
    sumadj = adj = centaj;
    wt = centwt;
    i = icent;
    goto S70;
S60:
    if(qtired(iterf) || qsmall(term)) goto S80;
S70:
/*
     Update weights for next higher J
*/
    wt *= (xnonc/(double)(i+1));
/*
     Calculate PTERM and add term to sum
*/
    pterm = pcent-sumadj;
    term = wt*pterm;
    sum += term;
/*
     Update adjustment term for DF for next iteration
*/
    i += 1;
    dfd2 = dg(i)/2.0e0;
    adj = adj*chid2/dfd2;
    sumadj += adj;
    iterf += 1;
    goto S60;
S80:
    *cum = sum;
    *ccum = 0.5e0+(0.5e0-*cum);
    return;
#undef dg
#undef qsmall
#undef qtired
} /* END */