00074 {
00075 #define dg(i) (*df+2.0e0*(double)(i))
00076 #define qsmall(xx) (int)(sum < 1.0e-20 || (xx) < eps*sum)
00077 #define qtired(i) (int)((i) > ntired)
00078 static double eps = 1.0e-5;
00079 static int ntired = 1000;
00080 static double adj,centaj,centwt,chid2,dfd2,lcntaj,lcntwt,lfact,pcent,pterm,sum,
00081     sumadj,term,wt,xnonc;
00082 static int i,icent,iterb,iterf;
00083 static double T1,T2,T3;
00084 /*
00085      ..
00086      .. Executable Statements ..
00087 */
00088     if(!(*x <= 0.0e0)) goto S10;
00089     *cum = 0.0e0;
00090     *ccum = 1.0e0;
00091     return;
00092 S10:
00093     if(!(*pnonc <= 1.0e-10)) goto S20;
00094 /*
00095      When non-centrality parameter is (essentially) zero,
00096      use cumulative chi-square distribution
00097 */
00098     cumchi(x,df,cum,ccum);
00099     return;
00100 S20:
00101     xnonc = *pnonc/2.0e0;
00102 /*
00103 **********************************************************************
00104      The following code calcualtes the weight, chi-square, and
00105      adjustment term for the central term in the infinite series.
00106      The central term is the one in which the poisson weight is
00107      greatest.  The adjustment term is the amount that must
00108      be subtracted from the chi-square to move up two degrees
00109      of freedom.
00110 **********************************************************************
00111 */
00112     icent = fifidint(xnonc);
00113     if(icent == 0) icent = 1;
00114     chid2 = *x/2.0e0;
00115 /*
00116      Calculate central weight term
00117 */
00118     T1 = (double)(icent+1);
00119     lfact = alngam(&T1);
00120     lcntwt = -xnonc+(double)icent*log(xnonc)-lfact;
00121     centwt = exp(lcntwt);
00122 /*
00123      Calculate central chi-square
00124 */
00125     T2 = dg(icent);
00126     cumchi(x,&T2,&pcent,ccum);
00127 /*
00128      Calculate central adjustment term
00129 */
00130     dfd2 = dg(icent)/2.0e0;
00131     T3 = 1.0e0+dfd2;
00132     lfact = alngam(&T3);
00133     lcntaj = dfd2*log(chid2)-chid2-lfact;
00134     centaj = exp(lcntaj);
00135     sum = centwt*pcent;
00136 /*
00137 **********************************************************************
00138      Sum backwards from the central term towards zero.
00139      Quit whenever either
00140      (1) the zero term is reached, or
00141      (2) the term gets small relative to the sum, or
00142      (3) More than NTIRED terms are totaled.
00143 **********************************************************************
00144 */
00145     iterb = 0;
00146     sumadj = 0.0e0;
00147     adj = centaj;
00148     wt = centwt;
00149     i = icent;
00150     goto S40;
00151 S30:
00152     if(qtired(iterb) || qsmall(term) || i == 0) goto S50;
00153 S40:
00154     dfd2 = dg(i)/2.0e0;
00155 /*
00156      Adjust chi-square for two fewer degrees of freedom.
00157      The adjusted value ends up in PTERM.
00158 */
00159     adj = adj*dfd2/chid2;
00160     sumadj += adj;
00161     pterm = pcent+sumadj;
00162 /*
00163      Adjust poisson weight for J decreased by one
00164 */
00165     wt *= ((double)i/xnonc);
00166     term = wt*pterm;
00167     sum += term;
00168     i -= 1;
00169     iterb += 1;
00170     goto S30;
00171 S50:
00172     iterf = 0;
00173 /*
00174 **********************************************************************
00175      Now sum forward from the central term towards infinity.
00176      Quit when either
00177      (1) the term gets small relative to the sum, or
00178      (2) More than NTIRED terms are totaled.
00179 **********************************************************************
00180 */
00181     sumadj = adj = centaj;
00182     wt = centwt;
00183     i = icent;
00184     goto S70;
00185 S60:
00186     if(qtired(iterf) || qsmall(term)) goto S80;
00187 S70:
00188 /*
00189      Update weights for next higher J
00190 */
00191     wt *= (xnonc/(double)(i+1));
00192 /*
00193      Calculate PTERM and add term to sum
00194 */
00195     pterm = pcent-sumadj;
00196     term = wt*pterm;
00197     sum += term;
00198 /*
00199      Update adjustment term for DF for next iteration
00200 */
00201     i += 1;
00202     dfd2 = dg(i)/2.0e0;
00203     adj = adj*chid2/dfd2;
00204     sumadj += adj;
00205     iterf += 1;
00206     goto S60;
00207 S80:
00208     *cum = sum;
00209     *ccum = 0.5e0+(0.5e0-*cum);
00210     return;
00211 #undef dg
00212 #undef qsmall
00213 #undef qtired
00214 } /* END */