pq = *p+*q;    if(!(fabs(pq-0.5e0-0.5e0) > 3.0e0*spmpar(&K1))) goto S180;    if(!(pq < 0.0e0)) goto S160;    *bound = 0.0e0;    goto S170;S160:    *bound = 1.0e0;S170:    *status = 3;    return;S190:S180:    if(*which == 1) goto S220;/*     Select the minimum of P or Q*/    qporq = *p <= *q;    if(!qporq) goto S200;    porq = *p;    goto S210;S200:    porq = *q;S220:S210:/*     Calculate ANSWERS*/    if(1 == *which) {/*     Calculating P and Q*/        *status = 0;        cumchi(x,df,p,q);        if(porq > 1.5e0) {            *status = 10;            return;        }    }    else if(2 == *which) {/*     Calculating X*/        *x = 5.0e0;        T3 = inf;        T6 = atol;        T7 = tol;        dstinv(&K2,&T3,&K4,&K4,&K5,&T6,&T7);        *status = 0;        dinvr(status,x,&fx,&qleft,&qhi);S230:        if(!(*status == 1)) goto S270;        cumchi(x,df,&cum,&ccum);        if(!qporq) goto S240;        fx = cum-*p;        goto S250;S240:        fx = ccum-*q;S250:        if(!(fx+porq > 1.5e0)) goto S260;        *status = 10;        return;S260:        dinvr(status,x,&fx,&qleft,&qhi);        goto S230;S270:        if(!(*status == -1)) goto S300;        if(!qleft) goto S280;        *status = 1;        *bound = 0.0e0;        goto S290;S280:        *status = 2;        *bound = inf;S300:S290:        ;    }    else if(3 == *which) {/*     Calculating DF*/        *df = 5.0e0;        T8 = zero;        T9 = inf;        T10 = atol;        T11 = tol;        dstinv(&T8,&T9,&K4,&K4,&K5,&T10,&T11);        *status = 0;        dinvr(status,df,&fx,&qleft,&qhi);S310:        if(!(*status == 1)) goto S350;        cumchi(x,df,&cum,&ccum);        if(!qporq) goto S320;        fx = cum-*p;        goto S330;S320:        fx = ccum-*q;S330:        if(!(fx+porq > 1.5e0)) goto S340;        *status = 10;        return;S340:        dinvr(status,df,&fx,&qleft,&qhi);        goto S310;S350:        if(!(*status == -1)) goto S380;        if(!qleft) goto S360;        *status = 1;        *bound = zero;        goto S370;S360:        *status = 2;        *bound = inf;S370:        ;    }S380:    return;#undef tol#undef atol#undef zero#undef inf}void cdfchn(int *which,double *p,double *q,double *x,double *df,	    double *pnonc,int *status,double *bound)/**********************************************************************      void cdfchn(int *which,double *p,double *q,double *x,double *df,            double *pnonc,int *status,double *bound)               Cumulative Distribution Function               Non-central Chi-Square                              Function     Calculates any one parameter of the non-central chi-square     distribution given values for the others.                              Arguments     WHICH --> Integer indicating which of the next three argument               values is to be calculated from the others.               Input range: 1..4               iwhich = 1 : Calculate P and Q from X and DF               iwhich = 2 : Calculate X from P,DF and PNONC               iwhich = 3 : Calculate DF from P,X and PNONC               iwhich = 3 : Calculate PNONC from P,X and DF     P <--> The integral from 0 to X of the non-central chi-square            distribution.            Input range: [0, 1-1E-16).     Q <--> 1-P.            Q is not used by this subroutine and is only included            for similarity with other cdf* routines.     X <--> Upper limit of integration of the non-central            chi-square distribution.            Input range: [0, +infinity).            Search range: [0,1E100]     DF <--> Degrees of freedom of the non-central             chi-square distribution.             Input range: (0, +infinity).             Search range: [ 1E-100, 1E100]     PNONC <--> Non-centrality parameter of the non-central                chi-square distribution.                Input range: [0, +infinity).                Search range: [0,1E4]     STATUS <-- 0 if calculation completed correctly               -I if input parameter number I is out of range                1 if answer appears to be lower than lowest                  search bound                2 if answer appears to be higher than greatest                  search bound     BOUND <-- Undefined if STATUS is 0               Bound exceeded by parameter number I if STATUS               is negative.               Lower search bound if STATUS is 1.               Upper search bound if STATUS is 2.                              Method     Formula  26.4.25   of   Abramowitz   and   Stegun,  Handbook  of     Mathematical  Functions (1966) is used to compute the cumulative     distribution function.     Computation of other parameters involve a seach for a value that     produces  the desired  value  of P.   The search relies  on  the     monotinicity of P with the other parameter.                            WARNING     The computation time  required for this  routine is proportional     to the noncentrality  parameter  (PNONC).  Very large  values of     this parameter can consume immense  computer resources.  This is     why the search range is bounded by 10,000.**********************************************************************/{#define tent4 1.0e4#define tol 1.0e-8#define atol 1.0e-50#define zero 1.0e-100#define one ( 1.0e0 - 1.0e-16 )#define inf 1.0e100static double K1 = 0.0e0;static double K3 = 0.5e0;static double K4 = 5.0e0;static double fx,cum,ccum;static unsigned long qhi,qleft;static double T2,T5,T6,T7,T8,T9,T10,T11,T12,T13;/*     ..     .. Executable Statements ..*//*     Check arguments*/    if(!(*which < 1 || *which > 4)) goto S30;    if(!(*which < 1)) goto S10;    *bound = 1.0e0;    goto S20;S10:    *bound = 4.0e0;S20:    *status = -1;    return;S30:    if(*which == 1) goto S70;/*     P*/    if(!(*p < 0.0e0 || *p > one)) goto S60;    if(!(*p < 0.0e0)) goto S40;    *bound = 0.0e0;    goto S50;S40:    *bound = one;S50:    *status = -2;    return;S70:S60:    if(*which == 2) goto S90;/*     X*/    if(!(*x < 0.0e0)) goto S80;    *bound = 0.0e0;    *status = -4;    return;S90:S80:    if(*which == 3) goto S110;/*     DF*/    if(!(*df <= 0.0e0)) goto S100;    *bound = 0.0e0;    *status = -5;    return;S110:S100:    if(*which == 4) goto S130;/*     PNONC*/    if(!(*pnonc < 0.0e0)) goto S120;    *bound = 0.0e0;    *status = -6;    return;S130:S120:/*     Calculate ANSWERS*/    if(1 == *which) {/*     Calculating P and Q*/        cumchn(x,df,pnonc,p,q);        *status = 0;    }    else if(2 == *which) {/*     Calculating X*/        *x = 5.0e0;        T2 = inf;        T5 = atol;        T6 = tol;        dstinv(&K1,&T2,&K3,&K3,&K4,&T5,&T6);        *status = 0;        dinvr(status,x,&fx,&qleft,&qhi);S140:        if(!(*status == 1)) goto S150;        cumchn(x,df,pnonc,&cum,&ccum);        fx = cum-*p;        dinvr(status,x,&fx,&qleft,&qhi);        goto S140;S150:        if(!(*status == -1)) goto S180;        if(!qleft) goto S160;        *status = 1;        *bound = 0.0e0;        goto S170;S160:        *status = 2;        *bound = inf;S180:S170:        ;    }    else if(3 == *which) {/*     Calculating DF*/        *df = 5.0e0;        T7 = zero;        T8 = inf;        T9 = atol;        T10 = tol;        dstinv(&T7,&T8,&K3,&K3,&K4,&T9,&T10);        *status = 0;        dinvr(status,df,&fx,&qleft,&qhi);S190:        if(!(*status == 1)) goto S200;        cumchn(x,df,pnonc,&cum,&ccum);        fx = cum-*p;        dinvr(status,df,&fx,&qleft,&qhi);        goto S190;S200:        if(!(*status == -1)) goto S230;        if(!qleft) goto S210;        *status = 1;        *bound = zero;        goto S220;S210:        *status = 2;        *bound = inf;S230:S220:        ;    }    else if(4 == *which) {/*     Calculating PNONC*/        *pnonc = 5.0e0;        T11 = tent4;        T12 = atol;        T13 = tol;        dstinv(&K1,&T11,&K3,&K3,&K4,&T12,&T13);        *status = 0;        dinvr(status,pnonc,&fx,&qleft,&qhi);S240:        if(!(*status == 1)) goto S250;        cumchn(x,df,pnonc,&cum,&ccum);        fx = cum-*p;        dinvr(status,pnonc,&fx,&qleft,&qhi);        goto S240;S250:        if(!(*status == -1)) goto S280;        if(!qleft) goto S260;        *status = 1;        *bound = zero;        goto S270;S260:        *status = 2;        *bound = tent4;S270:        ;    }S280:    return;#undef tent4#undef tol#undef atol#undef zero#undef one#undef inf}void cdff(int *which,double *p,double *q,double *f,double *dfn,	  double *dfd,int *status,double *bound)/**********************************************************************      void cdff(int *which,double *p,double *q,double *f,double *dfn,          double *dfd,int *status,double *bound)               Cumulative Distribution Function               F distribution                              Function     Calculates any one parameter of the F distribution     given values for the others.                              Arguments     WHICH --> Integer indicating which of the next four argument               values is to be calculated from the others.               Legal range: 1..4               iwhich = 1 : Calculate P and Q from F,DFN and DFD               iwhich = 2 : Calculate F from P,Q,DFN and DFD               iwhich = 3 : Calculate DFN from P,Q,F and DFD               iwhich = 4 : Calculate DFD from P,Q,F and DFN       P <--> The integral from 0 to F of the f-density.              Input range: [0,1].       Q <--> 1-P.              Input range: (0, 1].              P + Q = 1.0.       F <--> Upper limit of integration of the f-density.              Input range: [0, +infinity).              Search range: [0,1E100]     DFN < --> Degrees of freedom of the numerator sum of squares.               Input range: (0, +infinity).               Search range: [ 1E-100, 1E100]     DFD < --> Degrees of freedom of the denominator sum of squares.               Input range: (0, +infinity).               Search range: [ 1E-100, 1E100]     STATUS <-- 0 if calculation completed correctly               -I if input parameter number I is out of range                1 if answer appears to be lower than lowest                  search bound                2 if answer appears to be higher than greatest                  search bound                3 if P + Q .ne. 1     BOUND <-- Undefined if STATUS is 0               Bound exceeded by parameter number I if STATUS               is negative.               Lower search bound if STATUS is 1.               Upper search bound if STATUS is 2.                              Method     Formula   26.6.2   of   Abramowitz   and   Stegun,  Handbook  of     Mathematical  Functions (1966) is used to reduce the computation     of the  cumulative  distribution function for the  F  variate to     that of an incomplete beta.     Computation of other parameters involve a seach for a value that     produces  the desired  value  of P.   The search relies  on  the     monotinicity of P with the other parameter.                              WARNING     The value of the  cumulative  F distribution is  not necessarily     monotone in  either degrees of freedom.  There  thus may  be two     values  that  provide a given CDF  value.   This routine assumes     monotonicity and will find an arbitrary one of the two values.**********************************************************************/{#define tol 1.0e-8#define atol 1.0e-50#define zero 1.0e-100#define inf 1.0e100static int K1 = 1;static double K2 = 0.0e0;static double K4 = 0.5e0;static double K5 = 5.0e0;static double pq,fx,cum,ccum;static unsigned long qhi,qleft,qporq;static double T3,T6,T7,T8,T9,T10,T11,T12,T13,T14,T15;/*     ..     .. Executable Statements ..*//*     Check arguments*/    if(!(*which < 1 || *which > 4)) goto S30;    if(!(*which < 1)) goto S10;    *bound = 1.0e0;    goto S20;S10:    *bound = 4.0e0;S20:    *status = -1;    return;S30:    if(*which == 1) goto S70;/*     P*/    if(!(*p < 0.0e0 || *p > 1.0e0)) goto S60;    if(!(*p < 0.0e0)) goto S40;    *bound = 0.0e0;    goto S50;S40:    *bound = 1.0e0;S50:    *status = -2;    return;S70:S60:    if(*which == 1) goto S110;/*     Q*/    if(!(*q <= 0.0e0 || *q > 1.0e0)) goto S100;    if(!(*q <= 0.0e0)) goto S80;    *bound = 0.0e0;    goto S90;S80:    *bound = 1.0e0;S90:    *status = -3;    return;S110:S100:    if(*which == 2) goto S130;/*     F*/    if(!(*f < 0.0e0)) goto S120;    *bound = 0.0e0;    *status = -4;    return;S130:S120:    if(*which == 3) goto S150;/*     DFN*/    if(!(*dfn <= 0.0e0)) goto S140;    *bound = 0.0e0;    *status = -5;    return;S150:S140:    if(*which == 4) goto S170;/*     D