SUBROUTINE cdfchn(which,p,q,x,df,pnonc,status,bound)
C**********************************************************************
C
C SUBROUTINE CDFCHN( WHICH, P, Q, X, DF, PNONC, STATUS, BOUND )
C Cumulative Distribution Function
C Non-central Chi-Square
C
C
C Function
C
C
C Calculates any one parameter of the non-central chi-square
C distribution given values for the others.
C
C
C Arguments
C
C
C WHICH --> Integer indicating which of the next three argument
C values is to be calculated from the others.
C Input range: 1..4
C iwhich = 1 : Calculate P and Q from X and DF
C iwhich = 2 : Calculate X from P,DF and PNONC
C iwhich = 3 : Calculate DF from P,X and PNONC
C iwhich = 3 : Calculate PNONC from P,X and DF
C INTEGER WHICH
C
C P <--> The integral from 0 to X of the non-central chi-square
C distribution.
C Input range: [0, 1-1E-16).
C DOUBLE PRECISION P
C
C Q <--> 1-P.
C Q is not used by this subroutine and is only included
C for similarity with other cdf* routines.
C DOUBLE PRECISION Q
C
C X <--> Upper limit of integration of the non-central
C chi-square distribution.
C Input range: [0, +infinity).
C Search range: [0,1E300]
C DOUBLE PRECISION X
C
C DF <--> Degrees of freedom of the non-central
C chi-square distribution.
C Input range: (0, +infinity).
C Search range: [ 1E-300, 1E300]
C DOUBLE PRECISION DF
C
C PNONC <--> Non-centrality parameter of the non-central
C chi-square distribution.
C Input range: [0, +infinity).
C Search range: [0,1E4]
C DOUBLE PRECISION PNONC
C
C STATUS <-- 0 if calculation completed correctly
C -I if input parameter number I is out of range
C 1 if answer appears to be lower than lowest
C search bound
C 2 if answer appears to be higher than greatest
C search bound
C INTEGER STATUS
C
C BOUND <-- Undefined if STATUS is 0
C
C Bound exceeded by parameter number I if STATUS
C is negative.
C
C Lower search bound if STATUS is 1.
C
C Upper search bound if STATUS is 2.
C
C
C Method
C
C
C Formula 26.4.25 of Abramowitz and Stegun, Handbook of
C Mathematical Functions (1966) is used to compute the cumulative
C distribution function.
C
C Computation of other parameters involve a search for a value that
C produces the desired value of P. The search relies on the
C monotinicity of P with the other parameter.
C
C
C WARNING
C
C The computation time required for this routine is proportional
C to the noncentrality parameter (PNONC). Very large values of
C this parameter can consume immense computer resources. This is
C why the search range is bounded by 1e9.
C
C**********************************************************************
C .. Parameters ..
DOUBLE PRECISION tent9
PARAMETER (tent9=1.0D9)
DOUBLE PRECISION tol
PARAMETER (tol=1.0D-8)
DOUBLE PRECISION atol
PARAMETER (atol=1.0D-50)
DOUBLE PRECISION zero,one,inf
PARAMETER (zero=1.0D-300,one=1.0D0-1.0D-16,inf=1.0D300)
C ..
C .. Scalar Arguments ..
DOUBLE PRECISION bound,df,p,pnonc,q,x
INTEGER status,which
C ..
C .. Local Scalars ..
DOUBLE PRECISION ccum,cum,fx
LOGICAL qhi,qleft
C ..
C .. External Subroutines ..
EXTERNAL cumchn,dinvr,dstinv
C ..
IF (x.GT.inf) THEN
x = inf
END IF
IF (df.GT.inf) THEN
df = inf
END IF
IF (pnonc.GT.tent9) THEN
pnonc = tent9
END IF
IF (.NOT. ((which.LT.1).OR. (which.GT.4))) GO TO 30
IF (.NOT. (which.LT.1)) GO TO 10
bound = 1.0D0
GO TO 20
10 bound = 4.0D0
20 status = -1
RETURN
30 IF (which.EQ.1) GO TO 70
IF (.NOT. ((p.LT.0.0D0).OR. (p.GT.one))) GO TO 60
IF (.NOT. (p.LT.0.0D0)) GO TO 40
bound = 0.0D0
GO TO 50
40 bound = one
50 status = -2
RETURN
60 CONTINUE
70 IF (which.EQ.2) GO TO 90
IF (x.GE.0.0D0) GO TO 80
bound = 0.0D0
status = -4
RETURN
80 CONTINUE
90 IF (which.EQ.3) GO TO 110
IF (df.GT.0.0D0) GO TO 100
bound = 0.0D0
status = -5
RETURN
100 CONTINUE
110 IF (which.EQ.4) GO TO 130
IF (pnonc.GE.0.0D0) GO TO 120
bound = 0.0D0
status = -6
RETURN
120 CONTINUE
130 IF ((1).EQ. (which)) THEN
CALL cumchn(x,df,pnonc,p,q)
status = 0
ELSE IF ((2).EQ. (which)) THEN
x = 5.0D0
CALL dstinv(0.0D0,inf,0.5D0,0.5D0,5.0D0,atol,tol)
status = 0
CALL dinvr(status,x,fx,qleft,qhi)
140 IF (.NOT. (status.EQ.1)) GO TO 150
CALL cumchn(x,df,pnonc,cum,ccum)
fx = cum - p
CALL dinvr(status,x,fx,qleft,qhi)
GO TO 140
150 IF (.NOT. (status.EQ.-1)) GO TO 180
IF (.NOT. (qleft)) GO TO 160
status = 1
bound = 0.0D0
GO TO 170
160 status = 2
bound = inf
170 CONTINUE
180 CONTINUE
ELSE IF ((3).EQ. (which)) THEN
df = 5.0D0
CALL dstinv(zero,inf,0.5D0,0.5D0,5.0D0,atol,tol)
status = 0
CALL dinvr(status,df,fx,qleft,qhi)
190 IF (.NOT. (status.EQ.1)) GO TO 200
CALL cumchn(x,df,pnonc,cum,ccum)
fx = cum - p
CALL dinvr(status,df,fx,qleft,qhi)
GO TO 190
200 IF (.NOT. (status.EQ.-1)) GO TO 230
IF (.NOT. (qleft)) GO TO 210
status = 1
bound = zero
GO TO 220
210 status = 2
bound = inf
220 CONTINUE
230 CONTINUE
ELSE IF ((4).EQ. (which)) THEN
pnonc = 5.0D0
CALL dstinv(0.0D0,tent9,0.5D0,0.5D0,5.0D0,atol,tol)
status = 0
CALL dinvr(status,pnonc,fx,qleft,qhi)
240 IF (.NOT. (status.EQ.1)) GO TO 250
CALL cumchn(x,df,pnonc,cum,ccum)
fx = cum - p
CALL dinvr(status,pnonc,fx,qleft,qhi)
GO TO 240
250 IF (.NOT. (status.EQ.-1)) GO TO 280
IF (.NOT. (qleft)) GO TO 260
status = 1
bound = zero
GO TO 270
260 status = 2
bound = tent9
270 CONTINUE
280 END IF
RETURN
END