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.**********************************************************************//**********************************************************************      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.**********************************************************************//**********************************************************************      void cdffnc(int *which,double *p,double *q,double *f,double *dfn,	          double *dfd,double *phonc,int *status,double *bound)                Cumulative Distribution Function               Non-central F distribution                              Function     Calculates any one parameter of the Non-central F     distribution given values for the others.                              Arguments     WHICH --> Integer indicating which of the next five argument               values is to be calculated from the others.               Legal range: 1..5               iwhich = 1 : Calculate P and Q from F,DFN,DFD and PNONC               iwhich = 2 : Calculate F from P,Q,DFN,DFD and PNONC               iwhich = 3 : Calculate DFN from P,Q,F,DFD and PNONC               iwhich = 4 : Calculate DFD from P,Q,F,DFN and PNONC               iwhich = 5 : Calculate PNONC from P,Q,F,DFN and DFD       P <--> The integral from 0 to F of the non-central f-density.              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.       F <--> Upper limit of integration of the non-central 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.               Must be in range: (0, +infinity).               Input range: (0, +infinity).               Search range: [ 1E-100, 1E100]     PNONC <-> The non-centrality parameter               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                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.20   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.                              WARNING     The  value  of the  cumulative  noncentral 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.**********************************************************************//**********************************************************************      void cdfgam(int *which,double *p,double *q,double *x,                  double *shape,double *scale,int *status,double *bound)               Cumulative Distribution Function                         GAMma Distribution                              Function     Calculates any one parameter of the gamma     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 X,SHAPE and SCALE               iwhich = 2 : Calculate X from P,Q,SHAPE and SCALE               iwhich = 3 : Calculate SHAPE from P,Q,X and SCALE               iwhich = 4 : Calculate SCALE from P,Q,X and SHAPE     P <--> The integral from 0 to X of the gamma density.            Input range: [0,1].     Q <--> 1-P.            Input range: (0, 1].            P + Q = 1.0.     X <--> The upper limit of integration of the gamma density.            Input range: [0, +infinity).            Search range: [0,1E100]     SHAPE <--> The shape parameter of the gamma density.                Input range: (0, +infinity).                Search range: [1E-100,1E100]     SCALE <--> The scale parameter of the gamma density.                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                10 if the gamma or inverse gamma routine cannot                   compute the answer.  Usually happens only for                   X and SHAPE very large (gt 1E10 or more)     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     Cumulative distribution function (P) is calculated directly by     the code associated with:     DiDinato, A. R. and Morris, A. H. Computation of the  incomplete     gamma function  ratios  and their  inverse.   ACM  Trans.  Math.     Softw. 12 (1986), 377-393.     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.                              Note     The gamma density is proportional to       T**(SHAPE - 1) * EXP(- SCALE * T)**********************************************************************//**********************************************************************      void cdfnbn(int *which,double *p,double *q,double *s,double *xn,	          double *pr,double *ompr,int *status,double *bound)               Cumulative Distribution Function               Negative BiNomial distribution                              Function     Calculates any one parameter of the negative binomial     distribution given values for the others.     The  cumulative  negative   binomial  distribution  returns  the     probability that there  will be  F or fewer failures before  the     XNth success in binomial trials each of which has probability of     success PR.     The individual term of the negative binomial is the probability of     S failures before XN successes and is          Choose( S, XN+S-1 ) * PR^(XN) * (1-PR)^S                              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 S,XN,PR and OMPR               iwhich = 2 : Calculate S from P,Q,XN,PR and OMPR               iwhich = 3 : Calculate XN from P,Q,S,PR and OMPR               iwhich = 4 : Calculate PR and OMPR from P,Q,S and XN     P <--> The cumulation from 0 to S of the  negative            binomial distribution.            Input range: [0,1].     Q <--> 1-P.            Input range: (0, 1].            P + Q = 1.0.     S <--> The upper limit of cumulation of the binomial distribution.            There are F or fewer failures before the XNth success.            Input range: [0, +infinity).            Search range: [0, 1E100]