DOUBLE PRECISION FUNCTION ALGDIV(A,B) C----------------------------------------------------------------------- C C COMPUTATION OF LN(GAMMA(B)/GAMMA(A+B)) WHEN B .GE. 8 C C -------- C C IN THIS ALGORITHM, DEL(X) IS THE FUNCTION DEFINED BY C LN(GAMMA(X)) = (X - 0.5)*LN(X) - X + 0.5*LN(2*PI) + DEL(X). C C----------------------------------------------------------------------- C .. SCALAR ARGUMENTS .. DOUBLE PRECISION A,B C .. C .. LOCAL SCALARS .. DOUBLE PRECISION C,C0,C1,C2,C3,C4,C5,D,H,S11,S3,S5,S7,S9,T,U,V,W, + X,X2 C .. C .. EXTERNAL FUNCTIONS .. DOUBLE PRECISION ALNREL EXTERNAL ALNREL C .. C .. INTRINSIC FUNCTIONS .. INTRINSIC DLOG C .. C .. DATA STATEMENTS .. DATA C0/.833333333333333D-01/,C1/-.277777777760991D-02/, + C2/.793650666825390D-03/,C3/-.595202931351870D-03/, + C4/.837308034031215D-03/,C5/-.165322962780713D-02/ C .. C .. EXECUTABLE STATEMENTS .. C------------------------ IF (A.LE.B) GO TO 10 H = B/A C = 1.0D0/ (1.0D0+H) X = H/ (1.0D0+H) D = A + (B-0.5D0) GO TO 20 10 H = A/B C = H/ (1.0D0+H) X = 1.0D0/ (1.0D0+H) D = B + (A-0.5D0) C C SET SN = (1 - X**N)/(1 - X) C 20 X2 = X*X S3 = 1.0D0 + (X+X2) S5 = 1.0D0 + (X+X2*S3) S7 = 1.0D0 + (X+X2*S5) S9 = 1.0D0 + (X+X2*S7) S11 = 1.0D0 + (X+X2*S9) C C SET W = DEL(B) - DEL(A + B) C T = (1.0D0/B)**2 W = ((((C5*S11*T+C4*S9)*T+C3*S7)*T+C2*S5)*T+C1*S3)*T + C0 W = W* (C/B) C C COMBINE THE RESULTS C U = D*ALNREL(A/B) V = A* (DLOG(B)-1.0D0) IF (U.LE.V) GO TO 30 ALGDIV = (W-V) - U RETURN 30 ALGDIV = (W-U) - V RETURN END C------------------------------------------------------------ DOUBLE PRECISION FUNCTION ALNGAM(X) C C********************************************************************** C C DOUBLE PRECISION FUNCTION ALNGAM(X) C REAL LN OF THE GAMMA FUNCTION C C C FUNCTION C C C RETURNS THE NATURAL LOGARITHM OF GAMMA(X). C C C ARGUMENTS C C C X --> VALUE AT WHICH SCALED LOG GAMMA IS TO BE RETURNED C X IS DOUBLE PRECISION C C C METHOD C C C IF X .LE. 6.0, THEN USE RECURSION TO GET X BELOW 3 C THEN APPLY RATIONAL APPROXIMATION NUMBER 5236 OF C HART ET AL, COMPUTER APPROXIMATIONS, JOHN WILEY AND C SONS, NY, 1968. C C IF X .GT. 6.0, THEN USE RECURSION TO GET X TO AT LEAST 12 AND C THEN USE FORMULA 5423 OF THE SAME SOURCE. C C********************************************************************** C C .. PARAMETERS .. DOUBLE PRECISION HLN2PI PARAMETER (HLN2PI=0.91893853320467274178D0) C .. C .. SCALAR ARGUMENTS .. DOUBLE PRECISION X C .. C .. LOCAL SCALARS .. DOUBLE PRECISION OFFSET,PROD,XX INTEGER I,N C .. C .. LOCAL ARRAYS .. DOUBLE PRECISION COEF(5),SCOEFD(4),SCOEFN(9) C .. C .. EXTERNAL FUNCTIONS .. DOUBLE PRECISION REVLPL EXTERNAL REVLPL C .. C .. INTRINSIC FUNCTIONS .. INTRINSIC LOG,DBLE C .. C .. DATA STATEMENTS .. DATA SCOEFN(1)/0.62003838007127258804D2/, + SCOEFN(2)/0.36036772530024836321D2/, + SCOEFN(3)/0.20782472531792126786D2/, + SCOEFN(4)/0.6338067999387272343D1/, + SCOEFN(5)/0.215994312846059073D1/, + SCOEFN(6)/0.3980671310203570498D0/, + SCOEFN(7)/0.1093115956710439502D0/, + SCOEFN(8)/0.92381945590275995D-2/, + SCOEFN(9)/0.29737866448101651D-2/ DATA SCOEFD(1)/0.62003838007126989331D2/, + SCOEFD(2)/0.9822521104713994894D1/, + SCOEFD(3)/-0.8906016659497461257D1/, + SCOEFD(4)/0.1000000000000000000D1/ DATA COEF(1)/0.83333333333333023564D-1/, + COEF(2)/-0.27777777768818808D-2/, + COEF(3)/0.79365006754279D-3/,COEF(4)/-0.594997310889D-3/, + COEF(5)/0.8065880899D-3/ C .. C .. EXECUTABLE STATEMENTS .. IF (.NOT. (X.LE.6.0D0)) GO TO 70 PROD = 1.0D0 XX = X IF (.NOT. (X.GT.3.0D0)) GO TO 30 10 IF (.NOT. (XX.GT.3.0D0)) GO TO 20 XX = XX - 1.0D0 PROD = PROD*XX GO TO 10 20 CONTINUE 30 IF (.NOT. (X.LT.2.0D0)) GO TO 60 40 IF (.NOT. (XX.LT.2.0D0)) GO TO 50 PROD = PROD/XX XX = XX + 1.0D0 GO TO 40 50 CONTINUE 60 ALNGAM = REVLPL(SCOEFN,9,XX-2.0D0)/REVLPL(SCOEFD,4,XX-2.0D0) C C C COMPUTE RATIONAL APPROXIMATION TO GAMMA(X) C C ALNGAM = ALNGAM*PROD ALNGAM = LOG(ALNGAM) GO TO 110 70 OFFSET = HLN2PI C C C IF NECESSARY MAKE X AT LEAST 12 AND CARRY CORRECTION IN OFFSET C C N = 12.0D0 - X IF (.NOT. (N.GT.0)) GO TO 90 PROD = 1.0D0 DO 80,I = 1,N PROD = PROD* (X+DBLE(I-1)) 80 CONTINUE OFFSET = OFFSET - LOG(PROD) XX = X + DBLE(N) GO TO 100 90 XX = X C C C COMPUTE POWER SERIES C C 100 ALNGAM = REVLPL(COEF,5,1.0D0/(XX**2D0))/XX ALNGAM = ALNGAM + OFFSET + (XX-0.5D0)*LOG(XX) - XX 110 RETURN END C------------------------ DOUBLE PRECISION FUNCTION ALNREL(A) C----------------------------------------------------------------------- C EVALUATION OF THE FUNCTION LN(1 + A) C----------------------------------------------------------------------- C .. SCALAR ARGUMENTS .. DOUBLE PRECISION A C .. C .. LOCAL SCALARS .. DOUBLE PRECISION P1,P2,P3,Q1,Q2,Q3,T,T2,W,X C .. C .. INTRINSIC FUNCTIONS .. INTRINSIC ABS,DBLE,DLOG C .. C .. DATA STATEMENTS .. DATA P1/-.129418923021993D+01/,P2/.405303492862024D+00/, + P3/-.178874546012214D-01/ DATA Q1/-.162752256355323D+01/,Q2/.747811014037616D+00/, + Q3/-.845104217945565D-01/ C .. C .. EXECUTABLE STATEMENTS .. C-------------------------- IF (ABS(A).GT.0.375D0) GO TO 10 T = A/ (A+2.0D0) T2 = T*T W = (((P3*T2+P2)*T2+P1)*T2+1.0D0)/ (((Q3*T2+Q2)*T2+Q1)*T2+1.0D0) ALNREL = 2.0D0*T*W RETURN C 10 X = 1.D0 + DBLE(A) ALNREL = DLOG(X) RETURN END C--------------------------------------------------- DOUBLE PRECISION FUNCTION APSER(A,B,X,EPS) C----------------------------------------------------------------------- C APSER YIELDS THE INCOMPLETE BETA RATIO I(SUB(1-X))(B,A) FOR C A .LE. MIN(EPS,EPS*B), B*X .LE. 1, AND X .LE. 0.5. USED WHEN C A IS VERY SMALL. USE ONLY IF ABOVE INEQUALITIES ARE SATISFIED. C----------------------------------------------------------------------- C .. SCALAR ARGUMENTS .. DOUBLE PRECISION A,B,EPS,X C .. C .. LOCAL SCALARS .. DOUBLE PRECISION AJ,BX,C,G,J,S,T,TOL C .. C .. EXTERNAL FUNCTIONS .. DOUBLE PRECISION PSI EXTERNAL PSI C .. C .. INTRINSIC FUNCTIONS .. INTRINSIC ABS,DLOG C .. C .. DATA STATEMENTS .. C-------------------- DATA G/.577215664901533D0/ C .. C .. EXECUTABLE STATEMENTS .. C-------------------- BX = B*X T = X - BX IF (B*EPS.GT.2.D-2) GO TO 10 C = DLOG(X) + PSI(B) + G + T GO TO 20 10 C = DLOG(BX) + G + T C 20 TOL = 5.0D0*EPS*ABS(C) J = 1.0D0 S = 0.0D0 30 J = J + 1.0D0 T = T* (X-BX/J) AJ = T/J S = S + AJ IF (ABS(AJ).GT.TOL) GO TO 30 C APSER = -A* (C+S) RETURN END C------------------------------------------------------- DOUBLE PRECISION FUNCTION BASYM(A,B,LAMBDA,EPS) C----------------------------------------------------------------------- C ASYMPTOTIC EXPANSION FOR IX(A,B) FOR LARGE A AND B. C LAMBDA = (A + B)*Y - B AND EPS IS THE TOLERANCE USED. C IT IS ASSUMED THAT LAMBDA IS NONNEGATIVE AND THAT C A AND B ARE GREATER THAN OR EQUAL TO 15. C----------------------------------------------------------------------- C .. SCALAR ARGUMENTS .. DOUBLE PRECISION A,B,EPS,LAMBDA C .. C .. LOCAL SCALARS .. DOUBLE PRECISION BSUM,DSUM,E0,E1,F,H,H2,HN,J0,J1,R,R0,R1,S,SUM,T, + T0,T1,U,W,W0,Z,Z0,Z2,ZN,ZNM1 INTEGER I,IM1,IMJ,J,M,MM1,MMJ,N,NP1,NUM C .. C .. LOCAL ARRAYS .. DOUBLE PRECISION A0(21),B0(21),C(21),D(21) C .. C .. EXTERNAL FUNCTIONS .. DOUBLE PRECISION BCORR,ERFC1,RLOG1 EXTERNAL BCORR,ERFC1,RLOG1 C .. C .. INTRINSIC FUNCTIONS .. INTRINSIC ABS,EXP,SQRT C .. C .. DATA STATEMENTS .. C------------------------ C ****** NUM IS THE MAXIMUM VALUE THAT N CAN TAKE IN THE DO LOOP C ENDING AT STATEMENT 50. IT IS REQUIRED THAT NUM BE EVEN. C THE ARRAYS A0, B0, C, D HAVE DIMENSION NUM + 1. C C------------------------ C E0 = 2/SQRT(PI) C E1 = 2**(-3/2) C------------------------ DATA NUM/20/ DATA E0/1.12837916709551D0/,E1/.353553390593274D0/ C .. C .. EXECUTABLE STATEMENTS .. C------------------------ BASYM = 0.0D0 IF (A.GE.B) GO TO 10 H = A/B R0 = 1.0D0/ (1.0D0+H) R1 = (B-A)/B W0 = 1.0D0/SQRT(A* (1.0D0+H)) GO TO 20 10 H = B/A R0 = 1.0D0/ (1.0D0+H) R1 = (B-A)/A W0 = 1.0D0/SQRT(B* (1.0D0+H)) C 20 F = A*RLOG1(-LAMBDA/A) + B*RLOG1(LAMBDA/B) T = EXP(-F) IF (T.EQ.0.0D0) RETURN Z0 = SQRT(F) Z = 0.5D0* (Z0/E1) Z2 = F + F C A0(1) = (2.0D0/3.0D0)*R1 C(1) = -0.5D0*A0(1) D(1) = -C(1) J0 = (0.5D0/E0)*ERFC1(1,Z0) J1 = E1 SUM = J0 + D(1)*W0*J1 C S = 1.0D0 H2 = H*H HN = 1.0D0 W = W0 ZNM1 = Z ZN = Z2 DO 70 N = 2,NUM,2 HN = H2*HN A0(N) = 2.0D0*R0* (1.0D0+H*HN)/ (N+2.0D0) NP1 = N + 1 S = S + HN A0(NP1) = 2.0D0*R1*S/ (N+3.0D0) C DO 60 I = N,NP1 R = -0.5D0* (I+1.0D0) B0(1) = R*A0(1) DO 40 M = 2,I BSUM = 0.0D0 MM1 = M - 1 DO 30 J = 1,MM1 MMJ = M - J BSUM = BSUM + (J*R-MMJ)*A0(J)*B0(MMJ) 30 CONTINUE B0(M) = R*A0(M) + BSUM/M 40 CONTINUE C(I) = B0(I)/ (I+1.0D0) C DSUM = 0.0D0 IM1 = I - 1 DO 50 J = 1,IM1 IMJ = I - J DSUM = DSUM + D(IMJ)*C(J) 50 CONTINUE D(I) = - (DSUM+C(I)) 60 CONTINUE C J0 = E1*ZNM1 + (N-1.0D0)*J0 J1 = E1*ZN + N*J1 ZNM1 = Z2*ZNM1 ZN = Z2*ZN W = W0*W T0 = D(N)*W*J0 W = W0*W T1 = D(NP1)*W*J1 SUM = SUM + (T0+T1) IF ((ABS(T0)+ABS(T1)).LE.EPS*SUM) GO TO 80 70 CONTINUE C 80 U = EXP(-BCORR(A,B)) BASYM = E0*T*U*SUM RETURN END C------------------------------------------------------------ DOUBLE PRECISION FUNCTION BCORR(A0,B0) C----------------------------------------------------------------------- C C EVALUATION OF DEL(A0) + DEL(B0) - DEL(A0 + B0) WHERE C LN(GAMMA(A)) = (A - 0.5)*LN(A) - A + 0.5*LN(2*PI) + DEL(A). C IT IS ASSUMED THAT A0 .GE. 8 AND B0 .GE. 8. C C----------------------------------------------------------------------- C .. SCALAR ARGUMENTS .. DOUBLE PRECISION A0,B0 C .. C .. LOCAL SCALARS .. DOUBLE PRECISION A,B,C,C0,C1,C2,C3,C4,C5,H,S11,S3,S5,S7,S9,T,W,X, + X2 C .. C .. INTRINSIC FUNCTIONS .. INTRINSIC DMAX1,DMIN1 C .. C .. DATA STATEMENTS .. DATA C0/.833333333333333D-01/,C1/-.277777777760991D-02/, + C2/.793650666825390D-03/,C3/-.595202931351870D-03/, + C4/.837308034031215D-03/,C5/-.165322962780713D-02/ C .. C .. EXECUTABLE STATEMENTS .. C------------------------ A = DMIN1(A0,B0) B = DMAX1(A0,B0) C H = A/B C = H/ (1.0D0+H) X = 1.0D0/ (1.0D0+H) X2 = X*X C C SET SN = (1 - X**N)/(1 - X) C S3 = 1.0D0 + (X+X2) S5 = 1.0D0 + (X+X2*S3) S7 = 1.0D0 + (X+X2*S5) S9 = 1.0D0 + (X+X2*S7) S11 = 1.0D0 + (X+X2*S9) C C SET W = DEL(B) - DEL(A + B) C T = (1.0D0/B)**2 W = ((((C5*S11*T+C4*S9)*T+C3*S7)*T+C2*S5)*T+C1*S3)*T + C0 W = W* (C/B) C C COMPUTE DEL(A) + W C T = (1.0D0/A)**2 BCORR = (((((C5*T+C4)*T+C3)*T+C2)*T+C1)*T+C0)/A + W RETURN END C------------------------------------------------------------- DOUBLE PRECISION FUNCTION BETACN(A,B,X,IWHICH) IMPLICIT DOUBLE PRECISION (A-H,O-P,R-Z),INTEGER (I-N),LOGICAL (Q) C********************************************************************** C C DOUBLE PRECISION FUNCTION BETACN(A,B,X,IWHICH) C BETA CONSTANTS C C COMPUTES THE CONSTANTS NEEDED FOR THE CUMULATIVE BETA. WERE NESTE C PROCEDURES ALLOWED IN FORTRAN, THIS WOULD BE NESTED IN CUMBET. C C C FUNCTION C C C IF IWHICH .EQ. 1 C RETURNS: C A*LOG(X) + (B-1)*LOG(1-X) - DLNBET( A, B ) - LOG(A) C C IF IWHICH .EQ. 2 C RETURNS: C A*LOG(X) + B*LOG(1-X) - DLNBET( A, B ) - LOG(A) C C********************************************************************** C .. SCALAR ARGUMENTS .. DOUBLE PRECISION A,B,X INTEGER IWHICH C .. C .. LOCAL SCALARS .. DOUBLE PRECISION ALCON,ALCONR,ALNCON,AOLD,BOLD,TEMP,Z LOGICAL QREVRS,QSAME C .. C .. EXTERNAL FUNCTIONS .. DOUBLE PRECISION DLN1MX,DLNBET EXTERNAL DLN1MX,DLNBET C .. C .. INTRINSIC FUNCTIONS .. INTRINSIC LOG C .. C .. SAVE STATEMENT .. SAVE AOLD,BOLD,ALCON,ALCONR C .. C .. DATA STATEMENTS .. DATA AOLD/-1.0D0/,BOLD/-1.0D0/ C .. C .. EXECUTABLE STATEMENTS .. QSAME = AOLD .EQ. A .AND. BOLD .EQ. B QREVRS = AOLD .EQ. B .AND. BOLD .EQ. A IF (QSAME) THEN ALNCON = ALCON ELSE IF (QREVRS) THEN ALNCON = ALCONR ELSE QSAME = .TRUE. TEMP = DLNBET(A,B) ALCON = -TEMP - LOG(A) ALCONR = -TEMP - LOG(B) ALNCON = ALCON AOLD = A BOLD = B END IF IF (.NOT. (X.GT.0.1D0)) GO TO 10 Z = LOG(1.0D0-X) GO TO 20 10 Z = DLN1MX(X) 20 IF (.NOT. (IWHICH.EQ.1)) GO TO 30 BETACN = A*LOG(X) + (B-1.0D0)*Z + ALNCON GO TO 40 30 BETACN = A*LOG(X) + B*Z + ALNCON 40 RETURN END C-------------------------------------------------------- DOUBLE PRECISION FUNCTION BETALN(A0,B0) C----------------------------------------------------------------------- C EVALUATION OF THE LOGARITHM OF THE BETA FUNCTION C----------------------------------------------------------------------- C E = 0.5*LN(2*PI) C-------------------------- C .. SCALAR ARGUMENTS .. DOUBLE PRECISION A0,B0 C .. C .. LOCAL SCALARS .. DOUBLE PRECISION A,B,C,E,H,U,V,W,Z INTEGER I,N C .. C .. EXTERNAL FUNCTIONS .. DOUBLE PRECISION ALGDIV,ALNREL,BCORR,GAMLN,GSUMLN EXTERNAL ALGDIV,ALNREL,BCORR,GAMLN,GSUMLN C .. C .. INTRINSIC FUNCTIONS .. INTRINSIC DLOG,DMAX1,DMIN1 C .. C .. DATA STATEMENTS .. DATA E/.918938533204673D0/ C .. C .. EXECUTABLE STATEMENTS .. C-------------------------- A = DMIN1(A0,B0) B = DMAX1(A0,B0) IF (A.GE.8.0D0) GO TO 100 IF (A.GE.1.0D0) GO TO 20 C----------------------------------------------------------------------- C PROCEDURE WHEN A .LT. 1 C----------------------------------------------------------------------- IF (B.GE.8.0D0) GO TO 10 BETALN = GAMLN(A) + (GAMLN(B)-GAMLN(A+B)) RETURN 10 BETALN = GAMLN(A) + ALGDIV(A,B) RETURN C----------------------------------------------------------------------- C PROCEDURE WHEN 1 .LE. A .LT. 8 C----------------------------------------------------------------------- 20 IF (A.GT.2.0D0) GO TO 40 IF (B.GT.2.0D0) GO TO 30 BETALN = GAMLN(A) + GAMLN(B) - GSUMLN(A,B) RETURN 30 W = 0.0D0 IF (B.LT.8.0D0) GO TO 60 BETALN = GAMLN(A) + ALGDIV(A,B) RETURN C C REDUCTION OF A WHEN B .LE. 1000 C 40 IF (B.GT.1000.0D0) GO TO 80 N = A - 1.0D0 W = 1.0D0 DO 50 I = 1,N A = A - 1.0D0 H = A/B W = W* (H/ (1.0D0+H)) 50 CONTINUE W = DLOG(W) IF (B.LT.8.0D0) GO TO 60 BETALN = W + GAMLN(A) + ALGDIV(A,B) RETURN C C REDUCTION OF B WHEN B .LT. 8 C 60 N = B - 1.0D0 Z = 1.0D0 DO 70 I = 1,N B = B - 1.0D0 Z = Z* (B/ (A+B)) 70 CONTINUE BETALN = W + DLOG(Z) + (GAMLN(A)+ (GAMLN(B)-GSUMLN(A,B))) RETURN C C REDUCTION OF A WHEN B .GT. 1000 C 80 N = A - 1.0D0 W = 1.0D0 DO 90 I = 1,N A = A - 1.0D0 W = W* (A/ (1.0D0+A/B)) 90 CONTINUE BETALN = (DLOG(W)-N*DLOG(B)) + (GAMLN(A)+ALGDIV(A,B)) RETURN C----------------------------------------------------------------------- C PROCEDURE WHEN A .GE. 8 C----------------------------------------------------------------------- 100 W = BCORR(A,B) H = A/B C = H/ (1.0D0+H) U = - (A-0.5D0)*DLOG(C) V = B*ALNREL(H) IF (U.LE.V) GO TO 110 BETALN = (((-0.5D0*DLOG(B)+E)+W)-V) - U RETURN 110 BETALN = (((-0.5D0*DLOG(B)+E)+W)-U) - V RETURN END C----------------------------------------------------------------------- DOUBLE PRECISION FUNCTION BFRAC(A,B,X,Y,LAMBDA,EPS) C----------------------------------------------------------------------- C CONTINUED FRACTION EXPANSION FOR IX(A,B) WHEN A,B .GT. 1. C IT IS ASSUMED THAT LAMBDA = (A + B)*Y - B. C----------------------------------------------------------------------- C .. SCALAR ARGUMENTS .. DOUBLE PRECISION A,B,EPS,LAMBDA,X,Y C .. C .. LOCAL SCALARS .. DOUBLE PRECISION ALPHA,AN,ANP1,BETA,BN,BNP1,C,C0,C1,E,N,P,R,R0,S, + T,W,YP1 C .. C .. EXTERNAL FUNCTIONS .. DOUBLE PRECISION BRCOMP EXTERNAL BRCOMP C .. C .. INTRINSIC FUNCTIONS .. INTRINSIC ABS C .. C .. EXECUTABLE STATEMENTS .. C-------------------- BFRAC = BRCOMP(A,B,X,Y) IF (BFRAC.EQ.0.0D0) RETURN C C = 1.0D0 + LAMBDA C0 = B/A C1 = 1.0D0 + 1.0D0/A YP1 = Y + 1.0D0 C N = 0.0D0 P = 1.0D0 S = A + 1.0D0 AN = 0.0D0 BN = 1.0D0 ANP1 = 1.0D0 BNP1 = C/C1 R = C1/C C C CONTINUED FRACTION CALCULATION C 10 N = N + 1.0D0 T = N/A W = N* (B-N)*X E = A/S ALPHA = (P* (P+C0)*E*E)* (W*X) E = (1.0D0+T)/ (C1+T+T) BETA = N + W/S + E* (C+N*YP1) P = 1.0D0 + T S = S + 2.0D0 C C UPDATE AN, BN, ANP1, AND BNP1 C T = ALPHA*AN + BETA*ANP1 AN = ANP1 ANP1 = T T = ALPHA*BN + BETA*BNP1 BN = BNP1 BNP1 = T C R0 = R R = ANP1/BNP1 IF (ABS(R-R0).LE.EPS*R) GO TO 20 C C RESCALE AN, BN, ANP1, AND BNP1 C AN = AN/BNP1 BN = BN/BNP1 ANP1 = R BNP1 = 1.0D0 GO TO 10 C C TERMINATION C 20 BFRAC = BFRAC*R RETURN END C-------------------------------------------------------- SUBROUTINE BGRAT(A,B,X,Y,W,EPS,IERR) C----------------------------------------------------------------------- C ASYMPTOTIC EXPANSION FOR IX(A,B) WHEN A IS LARGER THAN B. C THE RESULT OF THE EXPANSION IS ADDED TO W. IT IS ASSUMED C THAT A .GE. 15 AND B .LE. 1. EPS IS THE TOLERANCE USED. C IERR IS A VARIABLE THAT REPORTS THE STATUS OF THE RESULTS. C----------------------------------------------------------------------- C .. SCALAR ARGUMENTS .. DOUBLE PRECISION A,B,EPS,W,X,Y INTEGER IERR C .. C .. LOCAL SCALARS .. DOUBLE PRECISION BM1,BP2N,CN,COEF,DJ,J,L,LNX,N2,NU,P,Q,R,S,SUM,T, + T2,U,V,Z INTEGER I,N,NM1 C .. C .. LOCAL ARRAYS .. DOUBLE PRECISION C(30),D(30) C .. C .. EXTERNAL FUNCTIONS .. DOUBLE PRECISION ALGDIV,ALNREL,GAM1 EXTERNAL ALGDIV,ALNREL,GAM1 C .. C .. EXTERNAL SUBROUTINES .. EXTERNAL GRAT1 C .. C .. INTRINSIC FUNCTIONS .. INTRINSIC ABS,DLOG,EXP C .. C .. EXECUTABLE STATEMENTS .. C BM1 = (B-0.5D0) - 0.5D0 NU = A + 0.5D0*BM1 IF (Y.GT.0.375D0) GO TO 10 LNX = ALNREL(-Y) GO TO 20 10 LNX = DLOG(X) 20 Z = -NU*LNX IF (B*Z.EQ.0.0D0) GO TO 70 C C COMPUTATION OF THE EXPANSION C SET R = EXP(-Z)*Z**B/GAMMA(B) C R = B* (1.0D0+GAM1(B))*EXP(B*DLOG(Z)) R = R*EXP(A*LNX)*EXP(0.5D0*BM1*LNX) U = ALGDIV(B,A) + B*DLOG(NU) U = R*EXP(-U) IF (U.EQ.0.0D0) GO TO 70 CALL GRAT1(B,Z,R,P,Q,EPS) C V = 0.25D0* (1.0D0/NU)**2 T2 = 0.25D0*LNX*LNX L = W/U J = Q/R SUM = J T = 1.0D0 CN = 1.0D0 N2 = 0.0D0 DO 50 N = 1,30 BP2N = B + N2 J = (BP2N* (BP2N+1.0D0)*J+ (Z+BP2N+1.0D0)*T)*V N2 = N2 + 2.0D0 T = T*T2 CN = CN/ (N2* (N2+1.0D0)) C(N) = CN S = 0.0D0 IF (N.EQ.1) GO TO 40 NM1 = N - 1 COEF = B - N DO 30 I = 1,NM1 S = S + COEF*C(I)*D(N-I) COEF = COEF + B 30 CONTINUE 40 D(N) = BM1*CN + S/N DJ = D(N)*J SUM = SUM + DJ IF (SUM.LE.0.0D0) GO TO 70 IF (ABS(DJ).LE.EPS* (SUM+L)) GO TO 60 50 CONTINUE C C ADD THE RESULTS TO W C 60 IERR = 0 W = W + U*SUM RETURN C C THE EXPANSION CANNOT BE COMPUTED C 70 IERR = 1 RETURN END C----------------------------------------------------------- DOUBLE PRECISION FUNCTION BINDEN(NS,N,P) C********************************************************************** C C DOUBLE PRECISION FUNCTION BINDEN(NS,N,P) C BINOMIAL DENSITY FUNCTION C C C FUNCTION C C C RETURNS THE PROBABILITY OF NS SUCCESSES IN N BINOMIAL TRIALS, C EACH OF WHICH HAS PROBABILITY P OF SUCCESS. C C C ARGUMENTS C C C NS --> NUMBER OF SUCCESSES FOR WHICH BINOMIAL DENSITY C IS CALCULATED. C NS IS INTEGER C C N --> NUMBER OF TRIALS C N IS INTEGER C C P --> PROBABILITY OF SUCCESS ON EACH BINOMIAL TRIAL. C P IS DOUBLE PRECISION C C********************************************************************** C C C VARIABLES C C DOUBLE PRECISION BICOEF DOUBLE PRECISION NLFACT,SLFACT,NMSLFC DOUBLE PRECISION XLFACT C********************************************************************** C C BICOEF --> BINOMIAL COEFFICIENT C(N,NS) C C NLFACT --> LOG(N!), CALCULATED BY XLFACT C C SLFACT --> LOG(NS!), CALCULATED BY XLFACT C C NMSLFC --> LOG((N-NS)!), CALCULATED BY XLFACT C C XLFACT --> RETURNS THE LOG OF THE FACTORIAL OF ITS ARGUMENT C C********************************************************************** C C C *** MAIN BODY OF PROGRAM C NLFACT = XLFACT(N) SLFACT = XLFACT(NS) NMSLFC = XLFACT(N-NS) C BICOEF = DEXP(NLFACT- (SLFACT+NMSLFC)) C BINDEN = BICOEF* (P**DBLE(NS))* ((1.0-P)**DBLE(N-NS)) C RETURN END C------------------------------------------------------------------ DOUBLE PRECISION FUNCTION BPSER(A,B,X,EPS) C----------------------------------------------------------------------- C POWER SERIES EXPANSION FOR EVALUATING IX(A,B) WHEN B .LE. 1 C OR B*X .LE. 0.7. EPS IS THE TOLERANCE USED. C----------------------------------------------------------------------- C .. SCALAR ARGUMENTS .. DOUBLE PRECISION A,B,EPS,X C .. C .. LOCAL SCALARS .. DOUBLE PRECISION A0,APB,B0,C,N,SUM,T,TOL,U,W,Z INTEGER I,M C .. C .. EXTERNAL FUNCTIONS .. DOUBLE PRECISION ALGDIV,BETALN,GAM1,GAMLN1 EXTERNAL ALGDIV,BETALN,GAM1,GAMLN1 C .. C .. INTRINSIC FUNCTIONS .. INTRINSIC ABS,DBLE,DLOG,DMAX1,DMIN1,EXP C .. C .. EXECUTABLE STATEMENTS .. C BPSER = 0.0D0 IF (X.EQ.0.0D0) RETURN C----------------------------------------------------------------------- C COMPUTE THE FACTOR X**A/(A*BETA(A,B)) C----------------------------------------------------------------------- A0 = DMIN1(A,B) IF (A0.LT.1.0D0) GO TO 10 Z = A*DLOG(X) - BETALN(A,B) BPSER = EXP(Z)/A GO TO 100 10 B0 = DMAX1(A,B) IF (B0.GE.8.0D0) GO TO 90 IF (B0.GT.1.0D0) GO TO 40 C C PROCEDURE FOR A0 .LT. 1 AND B0 .LE. 1 C BPSER = X**A IF (BPSER.EQ.0.0D0) RETURN C APB = A + B IF (APB.GT.1.0D0) GO TO 20 Z = 1.0D0 + GAM1(APB) GO TO 30 20 U = DBLE(A) + DBLE(B) - 1.D0 Z = (1.0D0+GAM1(U))/APB C 30 C = (1.0D0+GAM1(A))* (1.0D0+GAM1(B))/Z BPSER = BPSER*C* (B/APB) GO TO 100 C C PROCEDURE FOR A0 .LT. 1 AND 1 .LT. B0 .LT. 8 C 40 U = GAMLN1(A0) M = B0 - 1.0D0 IF (M.LT.1) GO TO 60 C = 1.0D0 DO 50 I = 1,M B0 = B0 - 1.0D0 C = C* (B0/ (A0+B0)) 50 CONTINUE U = DLOG(C) + U C 60 Z = A*DLOG(X) - U B0 = B0 - 1.0D0 APB = A0 + B0 IF (APB.GT.1.0D0) GO TO 70 T = 1.0D0 + GAM1(APB) GO TO 80 70 U = DBLE(A0) + DBLE(B0) - 1.D0 T = (1.0D0+GAM1(U))/APB 80 BPSER = EXP(Z)* (A0/A)* (1.0D0+GAM1(B0))/T GO TO 100 C C PROCEDURE FOR A0 .LT. 1 AND B0 .GE. 8 C 90 U = GAMLN1(A0) + ALGDIV(A0,B0) Z = A*DLOG(X) - U BPSER = (A0/A)*EXP(Z) 100 IF (BPSER.EQ.0.0D0 .OR. A.LE.0.1D0*EPS) RETURN C----------------------------------------------------------------------- C COMPUTE THE SERIES C----------------------------------------------------------------------- SUM = 0.0D0 N = 0.0D0 C = 1.0D0 TOL = EPS/A 110 N = N + 1.0D0 C = C* (0.5D0+ (0.5D0-B/N))*X W = C/ (A+N) SUM = SUM + W IF (ABS(W).GT.TOL) GO TO 110 BPSER = BPSER* (1.0D0+A*SUM) RETURN END C-------------------------------------------- SUBROUTINE BRATIO(A,B,X,Y,W,W1,IERR) C----------------------------------------------------------------------- C C EVALUATION OF THE INCOMPLETE BETA FUNCTION IX(A,B) C C -------------------- C C IT IS ASSUMED THAT A AND B ARE NONNEGATIVE, AND THAT X .LE. 1 C AND Y = 1 - X. BRATIO ASSIGNS W AND W1 THE VALUES C C W = IX(A,B) C W1 = 1 - IX(A,B) C C IERR IS A VARIABLE THAT REPORTS THE STATUS OF THE RESULTS. C IF NO INPUT ERRORS ARE DETECTED THEN IERR IS SET TO 0 AND C W AND W1 ARE COMPUTED. OTHERWISE, IF AN ERROR IS DETECTED, C THEN W AND W1 ARE ASSIGNED THE VALUE 0 AND IERR IS SET TO C ONE OF THE FOLLOWING VALUES ... C C IERR = 1 IF A OR B IS NEGATIVE C IERR = 2 IF A = B = 0 C IERR = 3 IF X .LT. 0 OR X .GT. 1 C IERR = 4 IF Y .LT. 0 OR Y .GT. 1 C IERR = 5 IF X + Y .NE. 1 C IERR = 6 IF X = A = 0 C IERR = 7 IF Y = B = 0 C C-------------------- C WRITTEN BY ALFRED H. MORRIS, JR. C NAVAL SURFACE WARFARE CENTER C DAHLGREN, VIRGINIA C REVISED ... NOV 1991 C----------------------------------------------------------------------- C .. SCALAR ARGUMENTS .. DOUBLE PRECISION A,B,W,W1,X,Y INTEGER IERR C .. C .. LOCAL SCALARS .. DOUBLE PRECISION A0,B0,EPS,LAMBDA,T,X0,Y0,Z INTEGER IERR1,IND,N C .. C .. EXTERNAL FUNCTIONS .. DOUBLE PRECISION APSER,BASYM,BFRAC,BPSER,BUP,FPSER,SPMPAR EXTERNAL APSER,BASYM,BFRAC,BPSER,BUP,FPSER,SPMPAR C .. C .. EXTERNAL SUBROUTINES .. EXTERNAL BGRAT C .. C .. INTRINSIC FUNCTIONS .. INTRINSIC ABS,DMAX1,DMIN1 C .. C .. EXECUTABLE STATEMENTS .. C----------------------------------------------------------------------- C C ****** EPS IS A MACHINE DEPENDENT CONSTANT. EPS IS THE SMALLEST C FLOATING POINT NUMBER FOR WHICH 1.0 + EPS .GT. 1.0 C EPS = SPMPAR(1) C C----------------------------------------------------------------------- W = 0.0D0 W1 = 0.0D0 IF (A.LT.0.0D0 .OR. B.LT.0.0D0) GO TO 270 IF (A.EQ.0.0D0 .AND. B.EQ.0.0D0) GO TO 280 IF (X.LT.0.0D0 .OR. X.GT.1.0D0) GO TO 290 IF (Y.LT.0.0D0 .OR. Y.GT.1.0D0) GO TO 300 Z = ((X+Y)-0.5D0) - 0.5D0 IF (ABS(Z).GT.3.0D0*EPS) GO TO 310 C IERR = 0 IF (X.EQ.0.0D0) GO TO 210 IF (Y.EQ.0.0D0) GO TO 230 IF (A.EQ.0.0D0) GO TO 240 IF (B.EQ.0.0D0) GO TO 220 C EPS = DMAX1(EPS,1.D-15) IF (DMAX1(A,B).LT.1.D-3*EPS) GO TO 260 C IND = 0 A0 = A B0 = B X0 = X Y0 = Y IF (DMIN1(A0,B0).GT.1.0D0) GO TO 40 C C PROCEDURE FOR A0 .LE. 1 OR B0 .LE. 1 C IF (X.LE.0.5D0) GO TO 10 IND = 1 A0 = B B0 = A X0 = Y Y0 = X C 10 IF (B0.LT.DMIN1(EPS,EPS*A0)) GO TO 90 IF (A0.LT.DMIN1(EPS,EPS*B0) .AND. B0*X0.LE.1.0D0) GO TO 100 IF (DMAX1(A0,B0).GT.1.0D0) GO TO 20 IF (A0.GE.DMIN1(0.2D0,B0)) GO TO 110 IF (X0**A0.LE.0.9D0) GO TO 110 IF (X0.GE.0.3D0) GO TO 120 N = 20 GO TO 140 C 20 IF (B0.LE.1.0D0) GO TO 110 IF (X0.GE.0.3D0) GO TO 120 IF (X0.GE.0.1D0) GO TO 30 IF ((X0*B0)**A0.LE.0.7D0) GO TO 110 30 IF (B0.GT.15.0D0) GO TO 150 N = 20 GO TO 140 C C PROCEDURE FOR A0 .GT. 1 AND B0 .GT. 1 C 40 IF (A.GT.B) GO TO 50 LAMBDA = A - (A+B)*X GO TO 60 50 LAMBDA = (A+B)*Y - B 60 IF (LAMBDA.GE.0.0D0) GO TO 70 IND = 1 A0 = B B0 = A X0 = Y Y0 = X LAMBDA = ABS(LAMBDA) C 70 IF (B0.LT.40.0D0 .AND. B0*X0.LE.0.7D0) GO TO 110 IF (B0.LT.40.0D0) GO TO 160 IF (A0.GT.B0) GO TO 80 IF (A0.LE.100.0D0) GO TO 130 IF (LAMBDA.GT.0.03D0*A0) GO TO 130 GO TO 200 80 IF (B0.LE.100.0D0) GO TO 130 IF (LAMBDA.GT.0.03D0*B0) GO TO 130 GO TO 200 C C EVALUATION OF THE APPROPRIATE ALGORITHM C 90 W = FPSER(A0,B0,X0,EPS) W1 = 0.5D0 + (0.5D0-W) GO TO 250 C 100 W1 = APSER(A0,B0,X0,EPS) W = 0.5D0 + (0.5D0-W1) GO TO 250 C 110 W = BPSER(A0,B0,X0,EPS) W1 = 0.5D0 + (0.5D0-W) GO TO 250 C 120 W1 = BPSER(B0,A0,Y0,EPS) W = 0.5D0 + (0.5D0-W1) GO TO 250 C 130 W = BFRAC(A0,B0,X0,Y0,LAMBDA,15.0D0*EPS) W1 = 0.5D0 + (0.5D0-W) GO TO 250 C 140 W1 = BUP(B0,A0,Y0,X0,N,EPS) B0 = B0 + N 150 CALL BGRAT(B0,A0,Y0,X0,W1,15.0D0*EPS,IERR1) W = 0.5D0 + (0.5D0-W1) GO TO 250 C 160 N = B0 B0 = B0 - N IF (B0.NE.0.0D0) GO TO 170 N = N - 1 B0 = 1.0D0 170 W = BUP(B0,A0,Y0,X0,N,EPS) IF (X0.GT.0.7D0) GO TO 180 W = W + BPSER(A0,B0,X0,EPS) W1 = 0.5D0 + (0.5D0-W) GO TO 250 C 180 IF (A0.GT.15.0D0) GO TO 190 N = 20 W = W + BUP(A0,B0,X0,Y0,N,EPS) A0 = A0 + N 190 CALL BGRAT(A0,B0,X0,Y0,W,15.0D0*EPS,IERR1) W1 = 0.5D0 + (0.5D0-W) GO TO 250 C 200 W = BASYM(A0,B0,LAMBDA,100.0D0*EPS) W1 = 0.5D0 + (0.5D0-W) GO TO 250 C C TERMINATION OF THE PROCEDURE C 210 IF (A.EQ.0.0D0) GO TO 320 220 W = 0.0D0 W1 = 1.0D0 RETURN C 230 IF (B.EQ.0.0D0) GO TO 330 240 W = 1.0D0 W1 = 0.0D0 RETURN C 250 IF (IND.EQ.0) RETURN T = W W = W1 W1 = T RETURN C C PROCEDURE FOR A AND B .LT. 1.E-3*EPS C 260 W = B/ (A+B) W1 = A/ (A+B) RETURN C C ERROR RETURN C 270 IERR = 1 RETURN 280 IERR = 2 RETURN 290 IERR = 3 RETURN 300 IERR = 4 RETURN 310 IERR = 5 RETURN 320 IERR = 6 RETURN 330 IERR = 7 RETURN END C------------------------------------------------------ DOUBLE PRECISION FUNCTION BRCMP1(MU,A,B,X,Y) C----------------------------------------------------------------------- C EVALUATION OF EXP(MU) * (X**A*Y**B/BETA(A,B)) C----------------------------------------------------------------------- C .. SCALAR ARGUMENTS .. DOUBLE PRECISION A,B,X,Y INTEGER MU C .. C .. LOCAL SCALARS .. DOUBLE PRECISION A0,APB,B0,C,CONST,E,H,LAMBDA,LNX,LNY,T,U,V,X0,Y0, + Z INTEGER I,N C .. C .. EXTERNAL FUNCTIONS .. DOUBLE PRECISION ALGDIV,ALNREL,BCORR,BETALN,ESUM,GAM1,GAMLN1,RLOG1 EXTERNAL ALGDIV,ALNREL,BCORR,BETALN,ESUM,GAM1,GAMLN1,RLOG1 C .. C .. INTRINSIC FUNCTIONS .. INTRINSIC ABS,DBLE,DLOG,DMAX1,DMIN1,EXP,SQRT C .. C .. DATA STATEMENTS .. C----------------- C CONST = 1/SQRT(2*PI) C----------------- DATA CONST/.398942280401433D0/ C .. C .. EXECUTABLE STATEMENTS .. C A0 = DMIN1(A,B) IF (A0.GE.8.0D0) GO TO 130 C IF (X.GT.0.375D0) GO TO 10 LNX = DLOG(X) LNY = ALNREL(-X) GO TO 30 10 IF (Y.GT.0.375D0) GO TO 20 LNX = ALNREL(-Y) LNY = DLOG(Y) GO TO 30 20 LNX = DLOG(X) LNY = DLOG(Y) C 30 Z = A*LNX + B*LNY IF (A0.LT.1.0D0) GO TO 40 Z = Z - BETALN(A,B) BRCMP1 = ESUM(MU,Z) RETURN C----------------------------------------------------------------------- C PROCEDURE FOR A .LT. 1 OR B .LT. 1 C----------------------------------------------------------------------- 40 B0 = DMAX1(A,B) IF (B0.GE.8.0D0) GO TO 120 IF (B0.GT.1.0D0) GO TO 70 C C ALGORITHM FOR B0 .LE. 1 C BRCMP1 = ESUM(MU,Z) IF (BRCMP1.EQ.0.0D0) RETURN C APB = A + B IF (APB.GT.1.0D0) GO TO 50 Z = 1.0D0 + GAM1(APB) GO TO 60 50 U = DBLE(A) + DBLE(B) - 1.D0 Z = (1.0D0+GAM1(U))/APB C 60 C = (1.0D0+GAM1(A))* (1.0D0+GAM1(B))/Z BRCMP1 = BRCMP1* (A0*C)/ (1.0D0+A0/B0) RETURN C C ALGORITHM FOR 1 .LT. B0 .LT. 8 C 70 U = GAMLN1(A0) N = B0 - 1.0D0 IF (N.LT.1) GO TO 90 C = 1.0D0 DO 80 I = 1,N B0 = B0 - 1.0D0 C = C* (B0/ (A0+B0)) 80 CONTINUE U = DLOG(C) + U C 90 Z = Z - U B0 = B0 - 1.0D0 APB = A0 + B0 IF (APB.GT.1.0D0) GO TO 100 T = 1.0D0 + GAM1(APB) GO TO 110 100 U = DBLE(A0) + DBLE(B0) - 1.D0 T = (1.0D0+GAM1(U))/APB 110 BRCMP1 = A0*ESUM(MU,Z)* (1.0D0+GAM1(B0))/T RETURN C C ALGORITHM FOR B0 .GE. 8 C 120 U = GAMLN1(A0) + ALGDIV(A0,B0) BRCMP1 = A0*ESUM(MU,Z-U) RETURN C----------------------------------------------------------------------- C PROCEDURE FOR A .GE. 8 AND B .GE. 8 C----------------------------------------------------------------------- 130 IF (A.GT.B) GO TO 140 H = A/B X0 = H/ (1.0D0+H) Y0 = 1.0D0/ (1.0D0+H) LAMBDA = A - (A+B)*X GO TO 150 140 H = B/A X0 = 1.0D0/ (1.0D0+H) Y0 = H/ (1.0D0+H) LAMBDA = (A+B)*Y - B C 150 E = -LAMBDA/A IF (ABS(E).GT.0.6D0) GO TO 160 U = RLOG1(E) GO TO 170 160 U = E - DLOG(X/X0) C 170 E = LAMBDA/B IF (ABS(E).GT.0.6D0) GO TO 180 V = RLOG1(E) GO TO 190 180 V = E - DLOG(Y/Y0) C 190 Z = ESUM(MU,- (A*U+B*V)) BRCMP1 = CONST*SQRT(B*X0)*Z*EXP(-BCORR(A,B)) RETURN END C---------------------------------------------------- DOUBLE PRECISION FUNCTION BRCOMP(A,B,X,Y) C----------------------------------------------------------------------- C EVALUATION OF X**A*Y**B/BETA(A,B) C----------------------------------------------------------------------- C .. SCALAR ARGUMENTS .. DOUBLE PRECISION A,B,X,Y C .. C .. LOCAL SCALARS .. DOUBLE PRECISION A0,APB,B0,C,CONST,E,H,LAMBDA,LNX,LNY,T,U,V,X0,Y0, + Z INTEGER I,N C .. C .. EXTERNAL FUNCTIONS .. DOUBLE PRECISION ALGDIV,ALNREL,BCORR,BETALN,GAM1,GAMLN1,RLOG1 EXTERNAL ALGDIV,ALNREL,BCORR,BETALN,GAM1,GAMLN1,RLOG1 C .. C .. INTRINSIC FUNCTIONS .. INTRINSIC ABS,DBLE,DLOG,DMAX1,DMIN1,EXP,SQRT C .. C .. DATA STATEMENTS .. C----------------- C CONST = 1/SQRT(2*PI) C----------------- DATA CONST/.398942280401433D0/ C .. C .. EXECUTABLE STATEMENTS .. C BRCOMP = 0.0D0 IF (X.EQ.0.0D0 .OR. Y.EQ.0.0D0) RETURN A0 = DMIN1(A,B) IF (A0.GE.8.0D0) GO TO 130 C IF (X.GT.0.375D0) GO TO 10 LNX = DLOG(X) LNY = ALNREL(-X) GO TO 30 10 IF (Y.GT.0.375D0) GO TO 20 LNX = ALNREL(-Y) LNY = DLOG(Y) GO TO 30 20 LNX = DLOG(X) LNY = DLOG(Y) C 30 Z = A*LNX + B*LNY IF (A0.LT.1.0D0) GO TO 40 Z = Z - BETALN(A,B) BRCOMP = EXP(Z) RETURN C----------------------------------------------------------------------- C PROCEDURE FOR A .LT. 1 OR B .LT. 1 C----------------------------------------------------------------------- 40 B0 = DMAX1(A,B) IF (B0.GE.8.0D0) GO TO 120 IF (B0.GT.1.0D0) GO TO 70 C C ALGORITHM FOR B0 .LE. 1 C BRCOMP = EXP(Z) IF (BRCOMP.EQ.0.0D0) RETURN C APB = A + B IF (APB.GT.1.0D0) GO TO 50 Z = 1.0D0 + GAM1(APB) GO TO 60 50 U = DBLE(A) + DBLE(B) - 1.D0 Z = (1.0D0+GAM1(U))/APB C 60 C = (1.0D0+GAM1(A))* (1.0D0+GAM1(B))/Z BRCOMP = BRCOMP* (A0*C)/ (1.0D0+A0/B0) RETURN C C ALGORITHM FOR 1 .LT. B0 .LT. 8 C 70 U = GAMLN1(A0) N = B0 - 1.0D0 IF (N.LT.1) GO TO 90 C = 1.0D0 DO 80 I = 1,N B0 = B0 - 1.0D0 C = C* (B0/ (A0+B0)) 80 CONTINUE U = DLOG(C) + U C 90 Z = Z - U B0 = B0 - 1.0D0 APB = A0 + B0 IF (APB.GT.1.0D0) GO TO 100 T = 1.0D0 + GAM1(APB) GO TO 110 100 U = DBLE(A0) + DBLE(B0) - 1.D0 T = (1.0D0+GAM1(U))/APB 110 BRCOMP = A0*EXP(Z)* (1.0D0+GAM1(B0))/T RETURN C C ALGORITHM FOR B0 .GE. 8 C 120 U = GAMLN1(A0) + ALGDIV(A0,B0) BRCOMP = A0*EXP(Z-U) RETURN C----------------------------------------------------------------------- C PROCEDURE FOR A .GE. 8 AND B .GE. 8 C----------------------------------------------------------------------- 130 IF (A.GT.B) GO TO 140 H = A/B X0 = H/ (1.0D0+H) Y0 = 1.0D0/ (1.0D0+H) LAMBDA = A - (A+B)*X GO TO 150 140 H = B/A X0 = 1.0D0/ (1.0D0+H) Y0 = H/ (1.0D0+H) LAMBDA = (A+B)*Y - B C 150 E = -LAMBDA/A IF (ABS(E).GT.0.6D0) GO TO 160 U = RLOG1(E) GO TO 170 160 U = E - DLOG(X/X0) C 170 E = LAMBDA/B IF (ABS(E).GT.0.6D0) GO TO 180 V = RLOG1(E) GO TO 190 180 V = E - DLOG(Y/Y0) C 190 Z = EXP(- (A*U+B*V)) BRCOMP = CONST*SQRT(B*X0)*Z*EXP(-BCORR(A,B)) RETURN END C------------------------------------------------------- DOUBLE PRECISION FUNCTION BUP(A,B,X,Y,N,EPS) C----------------------------------------------------------------------- C EVALUATION OF IX(A,B) - IX(A+N,B) WHERE N IS A POSITIVE INTEGER. C EPS IS THE TOLERANCE USED. C----------------------------------------------------------------------- C .. SCALAR ARGUMENTS .. DOUBLE PRECISION A,B,EPS,X,Y INTEGER N C .. C .. LOCAL SCALARS .. DOUBLE PRECISION AP1,APB,D,L,R,T,W INTEGER I,K,KP1,MU,NM1 C .. C .. EXTERNAL FUNCTIONS .. DOUBLE PRECISION BRCMP1,EXPARG EXTERNAL BRCMP1,EXPARG C .. C .. INTRINSIC FUNCTIONS .. INTRINSIC ABS,EXP C .. C .. EXECUTABLE STATEMENTS .. C C OBTAIN THE SCALING FACTOR EXP(-MU) AND C EXP(MU)*(X**A*Y**B/BETA(A,B))/A C APB = A + B AP1 = A + 1.0D0 MU = 0 D = 1.0D0 IF (N.EQ.1 .OR. A.LT.1.0D0) GO TO 10 IF (APB.LT.1.1D0*AP1) GO TO 10 MU = ABS(EXPARG(1)) K = EXPARG(0) IF (K.LT.MU) MU = K T = MU D = EXP(-T) C 10 BUP = BRCMP1(MU,A,B,X,Y)/A IF (N.EQ.1 .OR. BUP.EQ.0.0D0) RETURN NM1 = N - 1 W = D C C LET K BE THE INDEX OF THE MAXIMUM TERM C K = 0 IF (B.LE.1.0D0) GO TO 50 IF (Y.GT.1.D-4) GO TO 20 K = NM1 GO TO 30 20 R = (B-1.0D0)*X/Y - A IF (R.LT.1.0D0) GO TO 50 K = NM1 T = NM1 IF (R.LT.T) K = R C C ADD THE INCREASING TERMS OF THE SERIES C 30 DO 40 I = 1,K L = I - 1 D = ((APB+L)/ (AP1+L))*X*D W = W + D 40 CONTINUE IF (K.EQ.NM1) GO TO 70 C C ADD THE REMAINING TERMS OF THE SERIES C 50 KP1 = K + 1 DO 60 I = KP1,NM1 L = I - 1 D = ((APB+L)/ (AP1+L))*X*D W = W + D IF (D.LE.EPS*W) GO TO 70 60 CONTINUE C C TERMINATE THE PROCEDURE C 70 BUP = BUP*W RETURN END C------------------------------------------------------ SUBROUTINE CDFBET(WHICH,P,X,A,B,STATUS,BOUND) C********************************************************************** C C SUBROUTINE CDFBET( WHICH, P, X, A, B, STATUS, BOUND ) C CUMULATIVE DISTRIBUTION FUNCTION C BETA DISTRIBUTION C C C FUNCTION C C C CALCULATES ANY ONE PARAMETER OF THE BETA DISTRIBUTION GIVEN C VALUES FOR THE OTHERS. C C C ARGUMENTS C C C WHICH --> INTEGER INDICATING WHICH OF THE NEXT FOUR ARGUMENT C VALUES IS TO BE CALCULATED FROM THE OTHERS. C LEGAL RANGE: 1..4 C INTEGER WHICH C C P <--> THE INTEGRAL FROM 0 TO X OF THE CHI-SQUARE C DISTRIBUTION. C INPUT RANGE: [0, 1]. C DOUBLE PRECISION P C C X <--> UPPER LIMIT OF INTEGRATION OF BETA DENSITY. C INPUT RANGE: [0,1]. C SEARCH RANGE: [0,1] C DOUBLE PRECISION X C C A <--> THE FIRST PARAMETER OF THE BETA DENSITY. C INPUT RANGE: (0, +INFINITY). C SEARCH RANGE: [1E-30,1E30] C DOUBLE PRECISION A C C B <--> THE SECOND PARAMETER OF THE BETA DENSITY. C INPUT RANGE: (0, +INFINITY). C SEARCH RANGE: [1E-30,1E30] C DOUBLE PRECISION B 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 CUMULATIVE DISTRIBUTION FUNCTION (P) IS CALCULATED DIRECTLY BY C CODE ASSOCIATED WITH THE FOLLOWING REFERENCE. C C DIDINATO, A. R. AND MORRIS, A. H. ALGORITHM 708: SIGNIFICANT C DIGIT COMPUTATION OF THE INCOMPLETE BETA FUNCTION RATIOS. ACM C TRANS. MATH. SOFTW. 18 (1993), 360-373. C C COMPUTATION OF OTHER PARAMETERS INVOLVE A SEACH 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 NOTE C C C THE BETA DENSITY IS PROPORTIONAL TO C T^(A-1) * (1-T)^(B-1) C C********************************************************************** C .. PARAMETERS .. DOUBLE PRECISION TOL PARAMETER (TOL=1.0D-4) DOUBLE PRECISION ZERO,INF PARAMETER (ZERO=1.0D-30,INF=1.0D30) C .. C .. SCALAR ARGUMENTS .. DOUBLE PRECISION A,B,BOUND,P,X INTEGER STATUS,WHICH C .. C .. LOCAL SCALARS .. DOUBLE PRECISION AA,BB,FX,PP,XHI,XLO,ZZ LOGICAL QHI,QLEFT C .. C .. EXTERNAL FUNCTIONS .. DOUBLE PRECISION CUMBET EXTERNAL CUMBET C .. C .. EXTERNAL SUBROUTINES .. EXTERNAL INVR,STINVR,STZROR,ZROR C .. C .. STATEMENT FUNCTIONS .. DOUBLE PRECISION CUM C .. C .. STATEMENT FUNCTION DEFINITIONS .. CUM(ZZ,AA,BB,PP) = CUMBET(ZZ,AA,BB) - PP C .. C .. EXECUTABLE STATEMENTS .. C C CHECK ARGUMENTS C IF (.NOT. ((WHICH.LT.1).OR. (WHICH.GT.4))) GO TO 10 IF (WHICH.LT.1) THEN BOUND = 1.0D0 ELSE BOUND = 4.0D0 END IF STATUS = -1 RETURN 10 IF (WHICH.EQ.1) GO TO 30 C C P C IF (.NOT. ((P.LT.0.0).OR. (P.GT.1))) GO TO 20 IF (P.LT.0.0) THEN BOUND = 0.0D0 ELSE BOUND = 1.0D0 END IF STATUS = -2 RETURN 20 CONTINUE 30 IF (WHICH.EQ.2) GO TO 50 C C X C IF (.NOT. ((X.LT.0.0D0).OR. (X.GT.1D0))) GO TO 40 IF (X.LT.0.0) THEN BOUND = 0.0D0 ELSE BOUND = 1.0D0 END IF STATUS = -3 RETURN 40 CONTINUE 50 IF (WHICH.EQ.3) GO TO 70 C C A C IF (.NOT. (A.LE.0.0D0)) GO TO 60 BOUND = 0.0D0 STATUS = -4 RETURN 60 CONTINUE 70 IF (WHICH.EQ.4) GO TO 90 C C B C IF (.NOT. (B.LE.0.0D0)) GO TO 80 BOUND = 0.0D0 STATUS = -5 RETURN 80 CONTINUE C C CALCULATE ANSWERS C 90 IF ((1).EQ. (WHICH)) THEN C C CALCULATING P C P = CUMBET(X,A,B) STATUS = 0 ELSE IF ((2).EQ. (WHICH)) THEN C C CALCULATING X C CALL STZROR(0.0D0,1.0D0,TOL,TOL) STATUS = 0 CALL ZROR(STATUS,X,FX,XLO,XHI,QLEFT,QHI) 100 IF (.NOT. (STATUS.EQ.1)) GO TO 110 FX = CUM(X,A,B,P) CALL ZROR(STATUS,X,FX,XLO,XHI,QLEFT,QHI) GO TO 100 110 IF (.NOT. (STATUS.EQ.-1)) GO TO 140 IF (.NOT. (QLEFT)) GO TO 120 STATUS = 1 BOUND = 0.0D0 GO TO 130 120 STATUS = 2 BOUND = 1.0D0 130 CONTINUE 140 CONTINUE ELSE IF ((3).EQ. (WHICH)) THEN C C COMPUTING A C A = 5.0 CALL STINVR(ZERO,INF,0.5D0,0.5D0,5.0D0,TOL,TOL) STATUS = 0 CALL INVR(STATUS,A,FX,QLEFT,QHI) 150 IF (.NOT. (STATUS.EQ.1)) GO TO 160 FX = CUM(X,A,B,P) CALL INVR(STATUS,A,FX,QLEFT,QHI) GO TO 150 160 IF (.NOT. (STATUS.EQ.-1)) GO TO 190 IF (.NOT. (QLEFT)) GO TO 170 STATUS = 1 BOUND = ZERO GO TO 180 170 STATUS = 2 BOUND = INF 180 CONTINUE 190 CONTINUE ELSE IF ((4).EQ. (WHICH)) THEN C C COMPUTING B C B = 5.0D0 CALL STINVR(ZERO,INF,0.5D0,0.5D0,5.0D0,TOL,TOL) STATUS = 0 CALL INVR(STATUS,B,FX,QLEFT,QHI) 200 IF (.NOT. (STATUS.EQ.1)) GO TO 210 FX = CUM(X,A,B,P) CALL INVR(STATUS,B,FX,QLEFT,QHI) GO TO 200 210 IF (.NOT. (STATUS.EQ.-1)) GO TO 240 IF (.NOT. (QLEFT)) GO TO 220 STATUS = 1 BOUND = ZERO GO TO 230 220 STATUS = 2 BOUND = INF 230 CONTINUE 240 END IF RETURN END C------------------------------------------------------------- SUBROUTINE CDFBIN(WHICH,P,S,XN,PR,STATUS,BOUND) C********************************************************************** C C SUBROUTINE CDFBIN ( WHICH, P, S, XN, PR, STATUS, BOUND ) C CUMULATIVE DISTRIBUTION FUNCTION C BINOMIAL DISTRIBUTION C C C FUNCTION C C C CALCULATES ANY ONE PARAMETER OF THE BINOMIAL C DISTRIBUTION GIVEN VALUES FOR THE OTHERS. C C C ARGUMENTS C C C WHICH --> INTEGER INDICATING WHICH OF THE NEXT FOUR ARGUMENT C VALUES IS TO BE CALCULATED FROM THE OTHERS. C LEGAL RANGE: 1..4 C INTEGER WHICH C C P <--> THE CUMULATION FROM 0 TO S OF THE BINOMIAL DISTRIBUTION. C (PROBABLILITY OF S OR FEWER SUCCESSES IN XN TRIALS EACH C WITH PROBABILITY OF SUCCESS PR.) C INPUT RANGE: [0,1]. C DOUBLE PRECISION P C C S <--> THE NUMBER OF SUCCESSES OBSERVED. C INPUT RANGE: [0, XN] C SEARCH RANGE: [0, XN] C DOUBLE PRECISION S C C XN <--> THE NUMBER OF BINOMIAL TRIALS. C INPUT RANGE: (0, +INFINITY). C SEARCH RANGE: [1E-30, 1E30] C DOUBLE PRECISION XN C C PR <--> THE PROBABILITY OF SUCCESS IN EACH BINOMIAL TRIAL. C INPUT RANGE: [0,1]. C SEARCH RANGE: [0,1] C DOUBLE PRECISION PR 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.5.24 OF ABRAMOWITZ AND STEGUN, HANDBOOK OF C MATHEMATICAL FUNCTIONS (1966) IS USED TO REDUCE THE BINOMIAL C DISTRIBUTION TO THE CUMULATIVE INCOMPLETE BETA DISTRIBUTION. C C COMPUTATION OF OTHER PARAMETERS INVOLVE A SEACH 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********************************************************************** C .. PARAMETERS .. DOUBLE PRECISION TOL PARAMETER (TOL=1.0D-4) DOUBLE PRECISION ZERO,INF PARAMETER (ZERO=1.0D-30,INF=1.0D30) C .. C .. SCALAR ARGUMENTS .. DOUBLE PRECISION BOUND,P,PR,S,XN INTEGER STATUS,WHICH C .. C .. LOCAL SCALARS .. DOUBLE PRECISION FX,PP,PPR,SS,XHI,XLO,XNXN LOGICAL QHI,QLEFT C .. C .. EXTERNAL FUNCTIONS .. DOUBLE PRECISION CUMBIN EXTERNAL CUMBIN C .. C .. EXTERNAL SUBROUTINES .. EXTERNAL INVR,STINVR,STZROR,ZROR C .. C .. STATEMENT FUNCTIONS .. DOUBLE PRECISION CUM C .. C .. STATEMENT FUNCTION DEFINITIONS .. CUM(SS,XNXN,PPR,PP) = CUMBIN(SS,XNXN,PPR) - PP C .. C .. EXECUTABLE STATEMENTS .. C C CHECK ARGUMENTS C IF (.NOT. ((WHICH.LT.1).AND. (WHICH.GT.4))) GO TO 10 IF (WHICH.LT.1) THEN BOUND = 1.0D0 ELSE BOUND = 4.0D0 END IF STATUS = -1 RETURN 10 IF (WHICH.EQ.1) GO TO 30 C C P C IF (.NOT. ((P.LT.0.0D0).OR. (P.GT.1.0D0))) GO TO 20 IF (P.LT.0.0) THEN BOUND = 0.0D0 ELSE BOUND = 1.0D0 END IF STATUS = -2 RETURN 20 CONTINUE 30 IF (WHICH.EQ.3) GO TO 50 C C XN C IF (.NOT. (XN.LE.0.0D0)) GO TO 40 BOUND = 0.0D0 STATUS = -4 RETURN 40 CONTINUE 50 IF (WHICH.EQ.2) GO TO 70 C C S C IF (.NOT. ((S.LT.0.0D0).OR. ((WHICH.NE.3).AND. + (S.GT.XN)))) GO TO 60 IF (S.LT.0.0D0) THEN BOUND = 0.0D0 ELSE BOUND = XN END IF STATUS = -3 RETURN 60 CONTINUE 70 IF (WHICH.EQ.4) GO TO 90 C C PR C IF (.NOT. ((PR.LT.0.0D0).OR. (PR.GT.1.0D0))) GO TO 80 IF (PR.LT.0.0D0) THEN BOUND = 0.0D0 ELSE BOUND = 1.0D0 END IF STATUS = -5 RETURN 80 CONTINUE C C CALCULATE ANSWERS C 90 IF ((1).EQ. (WHICH)) THEN C C CALCULATING P C P = CUMBIN(S,XN,PR) STATUS = 0 ELSE IF ((2).EQ. (WHICH)) THEN C C CALCULATING S C S = 5.0D0 CALL STINVR(0.0D0,XN,0.5D0,0.5D0,5.0D0,TOL,TOL) STATUS = 0 CALL INVR(STATUS,S,FX,QLEFT,QHI) 100 IF (.NOT. (STATUS.EQ.1)) GO TO 110 FX = CUM(S,XN,PR,P) CALL INVR(STATUS,S,FX,QLEFT,QHI) GO TO 100 110 IF (.NOT. (STATUS.EQ.-1)) GO TO 140 IF (.NOT. (QLEFT)) GO TO 120 STATUS = 1 BOUND = 0.0D0 GO TO 130 120 STATUS = 2 BOUND = XN 130 CONTINUE 140 CONTINUE ELSE IF ((3).EQ. (WHICH)) THEN C C CALCULATING XN C XN = 5.0D0 CALL STINVR(ZERO,INF,0.5D0,0.5D0,5.0D0,TOL,TOL) STATUS = 0 CALL INVR(STATUS,XN,FX,QLEFT,QHI) 150 IF (.NOT. (STATUS.EQ.1)) GO TO 160 FX = CUM(S,XN,PR,P) CALL INVR(STATUS,XN,FX,QLEFT,QHI) GO TO 150 160 IF (.NOT. (STATUS.EQ.-1)) GO TO 190 IF (.NOT. (QLEFT)) GO TO 170 STATUS = 1 BOUND = ZERO GO TO 180 170 STATUS = 2 BOUND = INF 180 CONTINUE 190 CONTINUE ELSE IF ((4).EQ. (WHICH)) THEN C C CALCULATING PR C CALL STZROR(0.0D0,1.0D0,TOL,TOL) STATUS = 0 CALL ZROR(STATUS,PR,FX,XLO,XHI,QLEFT,QHI) 200 IF (.NOT. (STATUS.EQ.1)) GO TO 210 FX = CUM(S,XN,PR,P) CALL ZROR(STATUS,PR,FX,XLO,XHI,QLEFT,QHI) GO TO 200 210 IF (.NOT. (STATUS.EQ.-1)) GO TO 240 IF (.NOT. (QLEFT)) GO TO 220 STATUS = 1 BOUND = 0.0D0 GO TO 230 220 STATUS = 2 BOUND = 1.0D0 230 CONTINUE 240 END IF RETURN END C---------------------------------------------------------- SUBROUTINE CDFCHI(WHICH,P,X,DF,STATUS,BOUND) C********************************************************************** C C SUBROUTINE CDFCHI( WHICH, P, X, DF, STATUS, BOUND ) C CUMULATIVE DISTRIBUTION FUNCTION C CHI-SQUARE DISTRIBUTION C C C FUNCTION C C C CALCULATES ANY ONE PARAMETER OF THE 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 LEGAL RANGE: 1..3 C INTEGER WHICH C C P <--> THE INTEGRAL FROM 0 TO X OF THE CHI-SQUARE C DISTRIBUTION. C INPUT RANGE: [0, 1-1E-7). C DOUBLE PRECISION P C C X <--> UPPER LIMIT OF INTEGRATION OF THE NON-CENTRAL C CHI-SQUARE DISTRIBUTION. C INPUT RANGE: [0, +INFINITY). C SEARCH RANGE: [0,1E30] C DOUBLE PRECISION X C C DF <--> DEGREES OF FREEDOM OF THE C CHI-SQUARE DISTRIBUTION. C INPUT RANGE: (0, +INFINITY). C SEARCH RANGE: [ 1E-30, 1E30] C DOUBLE PRECISION DF 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 10 INDICATES ERROR RETURNED FROM CUMGAM. SEE C REFERENCES IN CDFGAM 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.19 OF ABRAMOWITZ AND STEGUN, HANDBOOK OF C MATHEMATICAL FUNCTIONS (1966) IS USED TO REDUCE THE CHISQURE C DISTRIBUTION TO THE INCOMPLETE DISTRIBUTION. C C COMPUTATION OF OTHER PARAMETERS INVOLVE A SEACH 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 .. PARAMETERS .. DOUBLE PRECISION TOL PARAMETER (TOL=1.0D-4) DOUBLE PRECISION ZERO,ONE,INF PARAMETER (ZERO=1.0D-30,ONE=1.0-1.0D-7,INF=1.0D30) C .. C .. SCALAR ARGUMENTS .. DOUBLE PRECISION BOUND,DF,P,X INTEGER STATUS,WHICH C .. C .. LOCAL SCALARS .. DOUBLE PRECISION DFDF,FX,PP,XX LOGICAL QHI,QLEFT C .. C .. EXTERNAL FUNCTIONS .. DOUBLE PRECISION CUMCHI EXTERNAL CUMCHI C .. C .. EXTERNAL SUBROUTINES .. EXTERNAL INVR,STINVR C .. C .. STATEMENT FUNCTIONS .. DOUBLE PRECISION CUM C .. C .. STATEMENT FUNCTION DEFINITIONS .. CUM(XX,DFDF,PP) = CUMCHI(XX,DFDF) - PP C .. C .. EXECUTABLE STATEMENTS .. C C CHECK ARGUMENTS C IF (.NOT. ((WHICH.LT.1).OR. (WHICH.GT.3))) GO TO 10 IF (WHICH.LT.1) THEN BOUND = 1.0D0 ELSE BOUND = 3.0D0 END IF STATUS = -1 RETURN 10 IF (WHICH.EQ.1) GO TO 30 C C P C IF (.NOT. ((P.LT.0.0D0).OR. (P.GT.ONE))) GO TO 20 IF (P.LT.0.0D0) THEN BOUND = 0.0D0 ELSE BOUND = ONE END IF STATUS = -2 RETURN 20 CONTINUE 30 IF (WHICH.EQ.2) GO TO 50 C C X C IF (.NOT. (X.LT.0.0D0)) GO TO 40 BOUND = 0.0D0 STATUS = -3 RETURN 40 CONTINUE 50 IF (WHICH.EQ.3) GO TO 70 C C DF C IF (.NOT. (DF.LE.0.0D0)) GO TO 60 BOUND = 0.0D0 STATUS = -4 RETURN 60 CONTINUE C C CALCULATE ANSWERS C 70 IF ((1).EQ. (WHICH)) THEN C C CALCULATING P C STATUS = 0 P = CUMCHI(X,DF) IF (P.GT.1.5D0) THEN STATUS = 10 RETURN END IF ELSE IF ((2).EQ. (WHICH)) THEN C C CALCULATING X C X = 5.0 CALL STINVR(0.0D0,INF,0.5D0,0.5D0,5.0D0,TOL,TOL) STATUS = 0 CALL INVR(STATUS,X,FX,QLEFT,QHI) 80 IF (.NOT. (STATUS.EQ.1)) GO TO 90 FX = CUM(X,DF,P) IF ((FX+P).GT.1.5) THEN STATUS = 10 RETURN END IF CALL INVR(STATUS,X,FX,QLEFT,QHI) GO TO 80 90 IF (.NOT. (STATUS.EQ.-1)) GO TO 120 IF (.NOT. (QLEFT)) GO TO 100 STATUS = 1 BOUND = 0.0 GO TO 110 100 STATUS = 2 BOUND = INF 110 CONTINUE 120 CONTINUE ELSE IF ((3).EQ. (WHICH)) THEN C C CALCULATING DF C DF = 5.0 CALL STINVR(ZERO,INF,0.5D0,0.5D0,5.0D0,TOL,TOL) STATUS = 0 CALL INVR(STATUS,DF,FX,QLEFT,QHI) 130 IF (.NOT. (STATUS.EQ.1)) GO TO 140 FX = CUM(X,DF,P) IF ((FX+P).GT.1.5) THEN STATUS = 10 RETURN END IF CALL INVR(STATUS,DF,FX,QLEFT,QHI) GO TO 130 140 IF (.NOT. (STATUS.EQ.-1)) GO TO 170 IF (.NOT. (QLEFT)) GO TO 150 STATUS = 1 BOUND = ZERO GO TO 160 150 STATUS = 2 BOUND = INF 160 CONTINUE 170 END IF RETURN END C----------------------------------------------------------- SUBROUTINE CDFCHN(WHICH,P,X,DF,PNONC,STATUS,BOUND) C********************************************************************** C C SUBROUTINE CDFCHN( WHICH, P, 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 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-7). C DOUBLE PRECISION P C C X <--> UPPER LIMIT OF INTEGRATION OF THE NON-CENTRAL C CHI-SQUARE DISTRIBUTION. C INPUT RANGE: [0, +INFINITY). C SEARCH RANGE: [0,1E30] 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-30, 1E30] 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 SEACH 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 10,000. C C********************************************************************** C .. PARAMETERS .. DOUBLE PRECISION TENT4 PARAMETER (TENT4=10000.0D0) DOUBLE PRECISION TOL PARAMETER (TOL=1.0D-4) DOUBLE PRECISION ZERO,ONE,INF PARAMETER (ZERO=1.0D-30,ONE=1.0D0-1.0D-7,INF=1.0D30) C .. C .. SCALAR ARGUMENTS .. DOUBLE PRECISION BOUND,DF,P,PNONC,X INTEGER STATUS,WHICH C .. C .. LOCAL SCALARS .. DOUBLE PRECISION DFDF,FX,PP,PPNONC,XX LOGICAL QHI,QLEFT C .. C .. EXTERNAL FUNCTIONS .. DOUBLE PRECISION CUMCHN EXTERNAL CUMCHN C .. C .. EXTERNAL SUBROUTINES .. EXTERNAL INVR,STINVR C .. C .. STATEMENT FUNCTIONS .. DOUBLE PRECISION CUM C .. C .. STATEMENT FUNCTION DEFINITIONS .. CUM(XX,DFDF,PPNONC,PP) = CUMCHN(XX,DFDF,PPNONC) - PP C .. C .. EXECUTABLE STATEMENTS .. C C CHECK ARGUMENTS C IF (.NOT. ((WHICH.LT.1).OR. (WHICH.GT.4))) GO TO 10 IF (WHICH.LT.1) THEN BOUND = 1.0D0 ELSE BOUND = 4.0D0 END IF STATUS = -1 RETURN 10 IF (WHICH.EQ.1) GO TO 30 C C P C IF (.NOT. ((P.LT.0.0D0).OR. (P.GT.ONE))) GO TO 20 IF (P.LT.0.0D0) THEN BOUND = 0.0D0 ELSE BOUND = ONE END IF STATUS = -2 RETURN 20 CONTINUE 30 IF (WHICH.EQ.2) GO TO 50 C C X C IF (.NOT. (X.LT.0.0D0)) GO TO 40 BOUND = 0.0D0 STATUS = -3 RETURN 40 CONTINUE 50 IF (WHICH.EQ.3) GO TO 70 C C DF C IF (.NOT. (DF.LE.0.0D0)) GO TO 60 BOUND = 0.0D0 STATUS = -4 RETURN 60 CONTINUE 70 IF (WHICH.EQ.4) GO TO 90 C C PNONC C IF (.NOT. (PNONC.LT.0.0D0)) GO TO 80 BOUND = 0.0D0 STATUS = -5 RETURN 80 CONTINUE C C CALCULATE ANSWERS C 90 IF ((1).EQ. (WHICH)) THEN C C CALCULATING P C P = CUMCHN(X,DF,PNONC) STATUS = 0 ELSE IF ((2).EQ. (WHICH)) THEN C C CALCULATING X C X = 5.0D0 CALL STINVR(0.0D0,INF,0.5D0,0.5D0,5.0D0,TOL,TOL) STATUS = 0 CALL INVR(STATUS,X,FX,QLEFT,QHI) 100 IF (.NOT. (STATUS.EQ.1)) GO TO 110 FX = CUM(X,DF,PNONC,P) CALL INVR(STATUS,X,FX,QLEFT,QHI) GO TO 100 110 IF (.NOT. (STATUS.EQ.-1)) GO TO 140 IF (.NOT. (QLEFT)) GO TO 120 STATUS = 1 BOUND = 0.0D0 GO TO 130 120 STATUS = 2 BOUND = INF 130 CONTINUE 140 CONTINUE ELSE IF ((3).EQ. (WHICH)) THEN C C CALCULATING DF C DF = 5.0D0 CALL STINVR(ZERO,INF,0.5D0,0.5D0,5.0D0,TOL,TOL) STATUS = 0 CALL INVR(STATUS,DF,FX,QLEFT,QHI) 150 IF (.NOT. (STATUS.EQ.1)) GO TO 160 FX = CUM(X,DF,PNONC,P) CALL INVR(STATUS,DF,FX,QLEFT,QHI) GO TO 150 160 IF (.NOT. (STATUS.EQ.-1)) GO TO 190 IF (.NOT. (QLEFT)) GO TO 170 STATUS = 1 BOUND = ZERO GO TO 180 170 STATUS = 2 BOUND = INF 180 CONTINUE 190 CONTINUE ELSE IF ((4).EQ. (WHICH)) THEN C C CALCULATING PNONC C PNONC = 5.0D0 CALL STINVR(0.0D0,TENT4,0.5D0,0.5D0,5.0D0,TOL,TOL) STATUS = 0 CALL INVR(STATUS,PNONC,FX,QLEFT,QHI) 200 IF (.NOT. (STATUS.EQ.1)) GO TO 210 FX = CUM(X,DF,PNONC,P) CALL INVR(STATUS,PNONC,FX,QLEFT,QHI) GO TO 200 210 IF (.NOT. (STATUS.EQ.-1)) GO TO 240 IF (.NOT. (QLEFT)) GO TO 220 STATUS = 1 BOUND = ZERO GO TO 230 220 STATUS = 2 BOUND = TENT4 230 CONTINUE 240 END IF RETURN END C----------------------------------------------------------------- SUBROUTINE CDFF(WHICH,P,F,DFN,DFD,STATUS,BOUND) C********************************************************************** C C SUBROUTINE CDFF( WHICH, P, F, DFN, DFD, STATUS, BOUND ) C CUMULATIVE DISTRIBUTION FUNCTION C F DISTRIBUTION C C C FUNCTION C C C CALCULATES ANY ONE PARAMETER OF THE F DISTRIBUTION C GIVEN VALUES FOR THE OTHERS. C C C ARGUMENTS C C C WHICH --> INTEGER INDICATING WHICH OF THE NEXT FOUR ARGUMENT C VALUES IS TO BE CALCULATED FROM THE OTHERS. C LEGAL RANGE: 1..4 C INTEGER WHICH C C P <--> THE INTEGRAL FROM 0 TO F OF THE F-DENSITY. C INPUT RANGE: [0,1-1E-7]. C DOUBLE PRECISION P C C F <--> UPPER LIMIT OF INTEGRATION OF THE F-DENSITY. C INPUT RANGE: [0, +INFINITY). C SEARCH RANGE: [0,1E30] C DOUBLE PRECISION F C C DFN < --> DEGREES OF FREEDOM OF THE NUMERATOR SUM OF SQUARES. C INPUT RANGE: (0, +INFINITY). C SEARCH RANGE: [ 1E-30, 1E30] C DOUBLE PRECISION DFN C C DFD < --> DEGREES OF FREEDOM OF THE DENOMINATOR SUM OF SQUARES. C INPUT RANGE: (0, +INFINITY). C SEARCH RANGE: [ 1E-30, 1E30] C DOUBLE PRECISION DFD 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.6.2 OF ABRAMOWITZ AND STEGUN, HANDBOOK OF C MATHEMATICAL FUNCTIONS (1966) IS USED TO REDUCE THE COMPUTATION C OF THE CUMULATIVE DISTRIBUTION FUNCTION FOR THE F VARIATE TO C THAT OF AN INCOMPLETE BETA. C C COMPUTATION OF OTHER PARAMETERS INVOLVE A SEACH 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 WARNING C C THE VALUE OF THE CUMULATIVE F DISTRIBUTION IS NOT NECESSARILY C MONOTONE IN EITHER DEGREES OF FREEDOM. THERE THUS MAY BE TWO C VALUES THAT PROVIDE A GIVEN CDF VALUE. THIS ROUTINE ASSUMES C MONOTONICITY AND WILL FIND AN ARBITRARY ONE OF THE TWO VALUES. C C********************************************************************** C .. PARAMETERS .. DOUBLE PRECISION TOL PARAMETER (TOL=1.0D-4) DOUBLE PRECISION ZERO,ONE,INF PARAMETER (ZERO=1.0D-30,ONE=1.0D0-1.0D-7,INF=1.0D30) C .. C .. SCALAR ARGUMENTS .. DOUBLE PRECISION BOUND,DFD,DFN,F,P INTEGER STATUS,WHICH C .. C .. LOCAL SCALARS .. DOUBLE PRECISION DFDD,DFNN,FF,FX,PP LOGICAL QHI,QLEFT C .. C .. EXTERNAL FUNCTIONS .. DOUBLE PRECISION CUMF EXTERNAL CUMF C .. C .. EXTERNAL SUBROUTINES .. EXTERNAL INVR,STINVR C .. C .. STATEMENT FUNCTIONS .. DOUBLE PRECISION CUM C .. C .. STATEMENT FUNCTION DEFINITIONS .. CUM(FF,DFNN,DFDD,PP) = CUMF(FF,DFNN,DFDD) - PP C .. C .. EXECUTABLE STATEMENTS .. C C CHECK ARGUMENTS C IF (.NOT. ((WHICH.LT.1).OR. (WHICH.GT.4))) GO TO 10 IF (WHICH.LT.1) THEN BOUND = 1.0D0 ELSE BOUND = 4.0D0 END IF STATUS = -1 RETURN 10 IF (WHICH.EQ.1) GO TO 30 C C P C IF (.NOT. ((P.LT.0.0D0).OR. (P.GT.ONE))) GO TO 20 IF (P.LT.0.0D0) THEN BOUND = 0.0D0 ELSE BOUND = ONE END IF STATUS = -2 RETURN 20 CONTINUE 30 IF (WHICH.EQ.2) GO TO 50 C C F C IF (.NOT. (F.LT.0.0D0)) GO TO 40 BOUND = 0.0D0 STATUS = -3 RETURN 40 CONTINUE 50 IF (WHICH.EQ.3) GO TO 70 C C DFN C IF (.NOT. (DFN.LE.0.0D0)) GO TO 60 BOUND = 0.0D0 STATUS = -4 RETURN 60 CONTINUE 70 IF (WHICH.EQ.4) GO TO 90 C C DFD C IF (.NOT. (DFD.LE.0.0D0)) GO TO 80 BOUND = 0.0D0 STATUS = -5 RETURN 80 CONTINUE C C CALCULATE ANSWERS C 90 IF ((1).EQ. (WHICH)) THEN C C CALCULATING P C P = CUMF(F,DFN,DFD) STATUS = 0 ELSE IF ((2).EQ. (WHICH)) THEN C C CALCULATING F C F = 5.0D0 CALL STINVR(0.0D0,INF,0.5D0,0.5D0,5.0D0,TOL,TOL) STATUS = 0 CALL INVR(STATUS,F,FX,QLEFT,QHI) 100 IF (.NOT. (STATUS.EQ.1)) GO TO 110 FX = CUM(F,DFN,DFD,P) CALL INVR(STATUS,F,FX,QLEFT,QHI) GO TO 100 110 IF (.NOT. (STATUS.EQ.-1)) GO TO 140 IF (.NOT. (QLEFT)) GO TO 120 STATUS = 1 BOUND = 0.0D0 GO TO 130 120 STATUS = 2 BOUND = INF 130 CONTINUE 140 CONTINUE ELSE IF ((3).EQ. (WHICH)) THEN C C CALCULATING DFN C DFN = 5.0D0 CALL STINVR(ZERO,INF,0.5D0,0.5D0,5.0D0,TOL,TOL) STATUS = 0 CALL INVR(STATUS,DFN,FX,QLEFT,QHI) 150 IF (.NOT. (STATUS.EQ.1)) GO TO 160 FX = CUM(F,DFN,DFD,P) CALL INVR(STATUS,DFN,FX,QLEFT,QHI) GO TO 150 160 IF (.NOT. (STATUS.EQ.-1)) GO TO 190 IF (.NOT. (QLEFT)) GO TO 170 STATUS = 1 BOUND = ZERO GO TO 180 170 STATUS = 2 BOUND = INF 180 CONTINUE 190 CONTINUE ELSE IF ((4).EQ. (WHICH)) THEN C C CALCULATING DFD C DFD = 5.0D0 CALL STINVR(ZERO,INF,0.5D0,0.5D0,5.0D0,TOL,TOL) STATUS = 0 CALL INVR(STATUS,DFD,FX,QLEFT,QHI) 200 IF (.NOT. (STATUS.EQ.1)) GO TO 210 FX = CUM(F,DFN,DFD,P) CALL INVR(STATUS,DFD,FX,QLEFT,QHI) GO TO 200 210 IF (.NOT. (STATUS.EQ.-1)) GO TO 240 IF (.NOT. (QLEFT)) GO TO 220 STATUS = 1 BOUND = ZERO GO TO 230 220 STATUS = 2 BOUND = INF 230 CONTINUE 240 END IF RETURN END C---------------------------------------------------------- SUBROUTINE CDFFNC(WHICH,P,F,DFN,DFD,PHONC,STATUS,BOUND) C********************************************************************** C C SUBROUTINE CDFFNC( WHICH, P, F, DFN, DFD, PNONC, STATUS, BOUND ) C CUMULATIVE DISTRIBUTION FUNCTION C NON-CENTRAL F DISTRIBUTION C C C FUNCTION C C C CALCULATES ANY ONE PARAMETER OF THE NON-CENTRAL F C DISTRIBUTION GIVEN VALUES FOR THE OTHERS. C C C ARGUMENTS C C C WHICH --> INTEGER INDICATING WHICH OF THE NEXT FOUR ARGUMENT C VALUES IS TO BE CALCULATED FROM THE OTHERS. C LEGAL RANGE: 1..5 C INTEGER WHICH C C P <--> THE INTEGRAL FROM 0 TO F OF THE NON-CENTRAL F-DENSITY. C INPUT RANGE: [0,1-1E-7]. C DOUBLE PRECISION P C C F <--> UPPER LIMIT OF INTEGRATION OF THE NON-CENTRAL F-DENSITY. C INPUT RANGE: [0, +INFINITY). C SEARCH RANGE: [0,1E30] C DOUBLE PRECISION F C C DFN < --> DEGREES OF FREEDOM OF THE NUMERATOR SUM OF SQUARES. C INPUT RANGE: (0, +INFINITY). C SEARCH RANGE: [ 1E-30, 1E30] C DOUBLE PRECISION DFN C C DFD < --> DEGREES OF FREEDOM OF THE DENOMINATOR SUM OF SQUARES. C MUST BE IN RANGE: (0, +INFINITY). C INPUT RANGE: (0, +INFINITY). C SEARCH RANGE: [ 1E-30, 1E30] C DOUBLE PRECISION DFD C C PNONC <-> THE NON-CENTRALITY PARAMETER C INPUT RANGE: [0,INFINITY) C SEARCH RANGE: [0,1E4] C DOUBLE PRECISION PHONC 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.6.20 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 SEACH 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 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 10,000. C C WARNING C C THE VALUE OF THE CUMULATIVE NONCENTRAL F DISTRIBUTION IS NOT C NECESSARILY MONOTONE IN EITHER DEGREES OF FREEDOM. THERE THUS C MAY BE TWO VALUES THAT PROVIDE A GIVEN CDF VALUE. THIS ROUTINE C ASSUMES MONOTONICITY AND WILL FIND AN ARBITRARY ONE OF THE TWO C VALUES. C C********************************************************************** C .. PARAMETERS .. DOUBLE PRECISION TENT4 PARAMETER (TENT4=10000.0D0) DOUBLE PRECISION TOL PARAMETER (TOL=1.0D-4) DOUBLE PRECISION ZERO,ONE,INF PARAMETER (ZERO=1.0D-30,ONE=1.0D0-1.0D-7,INF=1.0D30) C .. C .. SCALAR ARGUMENTS .. DOUBLE PRECISION BOUND,DFD,DFN,F,P,PHONC INTEGER STATUS,WHICH C .. C .. LOCAL SCALARS .. DOUBLE PRECISION DFDD,DFNN,FF,FX,PP,PPHONC LOGICAL QHI,QLEFT C .. C .. EXTERNAL FUNCTIONS .. DOUBLE PRECISION CUMFNC EXTERNAL CUMFNC C .. C .. EXTERNAL SUBROUTINES .. EXTERNAL INVR,STINVR C .. C .. STATEMENT FUNCTIONS .. DOUBLE PRECISION CUM C .. C .. STATEMENT FUNCTION DEFINITIONS .. CUM(FF,DFNN,DFDD,PP,PPHONC) = CUMFNC(FF,DFNN,DFDD,PPHONC) - PP C .. C .. EXECUTABLE STATEMENTS .. C C CHECK ARGUMENTS C IF (.NOT. ((WHICH.LT.1).OR. (WHICH.GT.5))) GO TO 10 IF (WHICH.LT.1) THEN BOUND = 1.0D0 ELSE BOUND = 5.0D0 END IF STATUS = -1 RETURN 10 IF (WHICH.EQ.1) GO TO 30 C C P C IF (.NOT. ((P.LT.0.0D0).OR. (P.GT.ONE))) GO TO 20 IF (P.LT.0.0D0) THEN BOUND = 0.0D0 ELSE BOUND = ONE END IF STATUS = -2 RETURN 20 CONTINUE 30 IF (WHICH.EQ.2) GO TO 50 C C F C IF (.NOT. (F.LT.0.0D0)) GO TO 40 BOUND = 0.0D0 STATUS = -3 RETURN 40 CONTINUE 50 IF (WHICH.EQ.3) GO TO 70 C C DFN C IF (.NOT. (DFN.LE.0.0D0)) GO TO 60 BOUND = 0.0D0 STATUS = -4 RETURN 60 CONTINUE 70 IF (WHICH.EQ.4) GO TO 90 C C DFD C IF (.NOT. (DFD.LE.0.0D0)) GO TO 80 BOUND = 0.0D0 STATUS = -5 RETURN 80 CONTINUE 90 IF (WHICH.EQ.5) GO TO 110 C C PHONC C IF (.NOT. (PHONC.LT.0.0D0)) GO TO 100 BOUND = 0.0D0 STATUS = -6 RETURN 100 CONTINUE C C CALCULATE ANSWERS C 110 IF ((1).EQ. (WHICH)) THEN C C CALCULATING P C P = CUMFNC(F,DFN,DFD,PHONC) STATUS = 0 ELSE IF ((2).EQ. (WHICH)) THEN C C CALCULATING F C F = 5.0D0 CALL STINVR(0.0D0,INF,0.5D0,0.5D0,5.0D0,TOL,TOL) STATUS = 0 CALL INVR(STATUS,F,FX,QLEFT,QHI) 120 IF (.NOT. (STATUS.EQ.1)) GO TO 130 FX = CUM(F,DFN,DFD,P,PHONC) CALL INVR(STATUS,F,FX,QLEFT,QHI) GO TO 120 130 IF (.NOT. (STATUS.EQ.-1)) GO TO 160 IF (.NOT. (QLEFT)) GO TO 140 STATUS = 1 BOUND = 0.0D0 GO TO 150 140 STATUS = 2 BOUND = INF 150 CONTINUE 160 CONTINUE ELSE IF ((3).EQ. (WHICH)) THEN C C CALCULATING DFN C DFN = 5.0D0 CALL STINVR(ZERO,INF,0.5D0,0.5D0,5.0D0,TOL,TOL) STATUS = 0 CALL INVR(STATUS,DFN,FX,QLEFT,QHI) 170 IF (.NOT. (STATUS.EQ.1)) GO TO 180 FX = CUM(F,DFN,DFD,P,PHONC) CALL INVR(STATUS,DFN,FX,QLEFT,QHI) GO TO 170 180 IF (.NOT. (STATUS.EQ.-1)) GO TO 210 IF (.NOT. (QLEFT)) GO TO 190 STATUS = 1 BOUND = ZERO GO TO 200 190 STATUS = 2 BOUND = INF 200 CONTINUE 210 CONTINUE ELSE IF ((4).EQ. (WHICH)) THEN C C CALCULATING DFD C DFD = 5.0D0 CALL STINVR(ZERO,INF,0.5D0,0.5D0,5.0D0,TOL,TOL) STATUS = 0 CALL INVR(STATUS,DFD,FX,QLEFT,QHI) 220 IF (.NOT. (STATUS.EQ.1)) GO TO 230 FX = CUM(F,DFN,DFD,P,PHONC) CALL INVR(STATUS,DFD,FX,QLEFT,QHI) GO TO 220 230 IF (.NOT. (STATUS.EQ.-1)) GO TO 260 IF (.NOT. (QLEFT)) GO TO 240 STATUS = 1 BOUND = ZERO GO TO 250 240 STATUS = 2 BOUND = INF 250 CONTINUE 260 CONTINUE ELSE IF ((5).EQ. (WHICH)) THEN C C CALCULATING PHONC C PHONC = 5.0D0 CALL STINVR(0.0D0,TENT4,0.5D0,0.5D0,5.0D0,TOL,TOL) STATUS = 0 CALL INVR(STATUS,PHONC,FX,QLEFT,QHI) 270 IF (.NOT. (STATUS.EQ.1)) GO TO 280 FX = CUM(F,DFN,DFD,P,PHONC) CALL INVR(STATUS,PHONC,FX,QLEFT,QHI) GO TO 270 280 IF (.NOT. (STATUS.EQ.-1)) GO TO 310 IF (.NOT. (QLEFT)) GO TO 290 STATUS = 1 BOUND = 0.0D0 GO TO 300 290 STATUS = 2 BOUND = TENT4 300 CONTINUE 310 END IF RETURN END C---------------------------------------------------- SUBROUTINE CDFGAM(WHICH,P,X,SHAPE,SCALE,STATUS,BOUND) C********************************************************************** C C SUBROUTINE CDFGAM( WHICH, P, X, SHAPE, SCALE, STATUS, BOUND ) C CUMULATIVE DISTRIBUTION FUNCTION C GAMMA DISTRIBUTION C C C FUNCTION C C C CALCULATES ANY ONE PARAMETER OF THE GAMMA C DISTRIBUTION GIVEN VALUES FOR THE OTHERS. C C C ARGUMENTS C C C WHICH --> INTEGER INDICATING WHICH OF THE NEXT FOUR ARGUMENT C VALUES IS TO BE CALCULATED FROM THE OTHERS. C LEGAL RANGE: 1..4 C INTEGER WHICH C C P <--> THE INTEGRAL FROM 0 TO X OF THE GAMMA DENSITY. C INPUT RANGE: [0,1-1E-7]. C DOUBLE PRECISION P C C X <--> THE UPPER LIMIT OF INTEGRATION OF THE GAMMA DENSITY. C INPUT RANGE: [0, +INFINITY). C SEARCH RANGE: [0,1E30] C DOUBLE PRECISION X C C SHAPE <--> THE SHAPE PARAMETER OF THE GAMMA DENSITY. C INPUT RANGE: (0, +INFINITY). C SEARCH RANGE: [1E-30,1E30] C DOUBLE PRECISION SHAPE C C C SCALE <--> THE SCALE PARAMETER OF THE GAMMA DENSITY. C INPUT RANGE: (0, +INFINITY). C SEARCH RANGE: (1E-30,1E30] C DOUBLE PRECISION SCALE 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 10 IF THE GAMMA OR INVERSE GAMMA ROUTINE CANNOT C COMPUTE THE ANSWER. USUALLY HAPPENS ONLY FOR C X AND SHAPE VERY LARGE (GT 1E10 OR MORE) 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 CUMULATIVE DISTRIBUTION FUNCTION (P) IS CALCULATED DIRECTLY BY C THE CODE ASSOCIATED WITH: C C DIDINATO, A. R. AND MORRIS, A. H. COMPUTATION OF THE INCOMPLETE C GAMMA FUNCTION RATIOS AND THEIR INVERSE. ACM TRANS. MATH. C SOFTW. 12 (1986), 377-393. C C COMPUTATION OF OTHER PARAMETERS INVOLVE A SEACH 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 NOTE C C C C THE GAMMA DENSITY IS PROPORTIONAL TO C T**(SHAPE - 1) * EXP(- SCALE * T) C C C********************************************************************** C .. PARAMETERS .. DOUBLE PRECISION TOL PARAMETER (TOL=1.0E-4) DOUBLE PRECISION ZERO,ONE,INF PARAMETER (ZERO=1.0E-30,ONE=1.0D0-1.0E-7,INF=1.0E30) C .. C .. SCALAR ARGUMENTS .. DOUBLE PRECISION BOUND,P,SCALE,SHAPE,X INTEGER STATUS,WHICH C .. C .. LOCAL SCALARS .. DOUBLE PRECISION XX DOUBLE PRECISION AA,FX,PP,XSCALE,ZZ INTEGER IERR LOGICAL QHI,QLEFT C .. C .. EXTERNAL FUNCTIONS .. DOUBLE PRECISION CUMGAM EXTERNAL CUMGAM C .. C .. EXTERNAL SUBROUTINES .. EXTERNAL GAMINV,INVR,STINVR C .. C .. INTRINSIC FUNCTIONS .. INTRINSIC DBLE C .. C .. STATEMENT FUNCTIONS .. DOUBLE PRECISION CUM C .. C .. STATEMENT FUNCTION DEFINITIONS .. CUM(ZZ,AA,PP) = CUMGAM(ZZ,AA) - PP C .. C .. EXECUTABLE STATEMENTS .. C C CHECK ARGUMENTS C IF (.NOT. ((WHICH.LT.1).OR. (WHICH.GT.4))) GO TO 10 IF (WHICH.LT.1) THEN BOUND = 1.0D0 ELSE BOUND = 4.0D0 END IF STATUS = -1 RETURN 10 IF (WHICH.EQ.1) GO TO 30 C C P C IF (.NOT. ((P.LT.0.0D0).OR. (P.GT.ONE))) GO TO 20 IF (P.LT.0.0D0) THEN BOUND = 0.0D0 ELSE BOUND = ONE END IF STATUS = -2 RETURN 20 CONTINUE 30 IF (WHICH.EQ.2) GO TO 50 C C X C IF (.NOT. (X.LT.0.0D0)) GO TO 40 BOUND = 0.0D0 STATUS = -3 RETURN 40 CONTINUE 50 IF (WHICH.EQ.3) GO TO 70 C C SHAPE C IF (.NOT. (SHAPE.LE.0.0D0)) GO TO 60 BOUND = 0.0D0 STATUS = -4 RETURN 60 CONTINUE 70 IF (WHICH.EQ.4) GO TO 90 C C SCALE C IF (.NOT. (SCALE.LE.0.0D0)) GO TO 80 BOUND = 0.0D0 STATUS = -5 RETURN 80 CONTINUE C C CALCULATE ANSWERS C 90 IF ((1).EQ. (WHICH)) THEN C C CALCULATING P C STATUS = 0 XSCALE = X*SCALE P = CUMGAM(XSCALE,SHAPE) IF (P.GT.1.5) STATUS = 10 ELSE IF ((2).EQ. (WHICH)) THEN C C COMPUTING X C CALL GAMINV(SHAPE,XX,-1.0D0,P,1.0D0-P,IERR) IF (IERR.LT.0.0D0) THEN STATUS = 10 RETURN ELSE X = XX/SCALE STATUS = 0 END IF ELSE IF ((3).EQ. (WHICH)) THEN C C COMPUTING SHAPE C SHAPE = 5.0D0 XSCALE = X*SCALE CALL STINVR(ZERO,INF,0.5D0,0.5D0,5.0D0,TOL,TOL) STATUS = 0 CALL INVR(STATUS,SHAPE,FX,QLEFT,QHI) 100 IF (.NOT. (STATUS.EQ.1)) GO TO 110 FX = CUM(XSCALE,SHAPE,P) IF ((FX+P).GT.1.5D0) THEN STATUS = 10 RETURN END IF CALL INVR(STATUS,SHAPE,FX,QLEFT,QHI) GO TO 100 110 IF (.NOT. (STATUS.EQ.-1)) GO TO 140 IF (.NOT. (QLEFT)) GO TO 120 STATUS = 1 BOUND = ZERO GO TO 130 120 STATUS = 2 BOUND = INF 130 CONTINUE 140 CONTINUE ELSE IF ((4).EQ. (WHICH)) THEN C C COMPUTING SCALE C CALL GAMINV(SHAPE,XX,-1.0D0,P,1.0D0-P,IERR) IF (IERR.LT.0.0D0) THEN STATUS = 10 RETURN ELSE SCALE = XX/X STATUS = 0 END IF END IF RETURN END C------------------------------------------------------ SUBROUTINE CDFNBN(WHICH,P,S,XN,PR,STATUS,BOUND) C********************************************************************** C C SUBROUTINE CDFNBN ( WHICH, P, S, XN, PR, STATUS, BOUND ) C CUMULATIVE DISTRIBUTION FUNCTION C NEGATIVE BINOMIAL DISTRIBUTION C C C FUNCTION C C C CALCULATES ANY ONE PARAMETER OF THE NEGATIVE BINOMIAL C DISTRIBUTION GIVEN VALUES FOR THE OTHERS. C C THE CUMULATIVE NEGATIVE BINOMIAL DISTRIBUTION RETURNS THE C PROBABILITY THAT THERE WILL BE F OR FEWER FAILURES BEFORE THE C XNTH SUCCESS IN BINOMIAL TRIALS EACH OF WHICH HAS PROBABILITY OF C SUCCESS PR. C C THE INDIVIDUAL TERM OF THE NEGATIVE BINOMIAL IS THE PROBABILITY OF C S FAILURES BEFORE XN SUCCESSES AND IS C CHOOSE( S, XN+S-1 ) * PR^(XN) * (1-PR)^S C C C ARGUMENTS C C C WHICH --> INTEGER INDICATING WHICH OF THE NEXT FOUR ARGUMENT C VALUES IS TO BE CALCULATED FROM THE OTHERS. C LEGAL RANGE: 1..4 C INTEGER WHICH C C P <--> THE CUMULATION FROM 0 TO S OF THE NEGATIVE C BINOMIAL DISTRIBUTION. C INPUT RANGE: [0,1-1E-7]. C DOUBLE PRECISION P C C S <--> THE UPPER LIMIT OF CUMULATION OF THE BINOMIAL DISTRIBUTION. C THERE ARE F OR FEWER FAILURES BEFORE THE XNTH SUCCESS. C INPUT RANGE: [0, +INFINITY). C SEARCH RANGE: [0, 1E30] C DOUBLE PRECISION S C C XN <--> THE NUMBER OF SUCCESSES. C INPUT RANGE: [0, +INFINITY). C SEARCH RANGE: [0, 1E30] C DOUBLE PRECISION XN C C PR <--> THE PROBABILITY OF SUCCESS IN EACH BINOMIAL TRIAL. C INPUT RANGE: [0,1]. C SEARCH RANGE: [0,1]. C DOUBLE PRECISION PR 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.5.26 OF ABRAMOWITZ AND STEGUN, HANDBOOK OF C MATHEMATICAL FUNCTIONS (1966) IS USED TO REDUCE CALCULATION OF C THE CUMULATIVE DISTRIBUTION FUNCTION TO THAT OF AN INCOMPLETE C BETA. C C COMPUTATION OF OTHER PARAMETERS INVOLVE A SEACH 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********************************************************************** C .. PARAMETERS .. DOUBLE PRECISION TOL PARAMETER (TOL=1.0D-4) DOUBLE PRECISION ONE,INF PARAMETER (ONE=1.0D0-1.0D-7,INF=1.0D30) C .. C .. SCALAR ARGUMENTS .. DOUBLE PRECISION BOUND,P,PR,S,XN INTEGER STATUS,WHICH C .. C .. LOCAL SCALARS .. DOUBLE PRECISION FX,PP,PPR,SS,XHI,XLO,XNXN LOGICAL QHI,QLEFT C .. C .. EXTERNAL FUNCTIONS .. DOUBLE PRECISION CUMNBN EXTERNAL CUMNBN C .. C .. EXTERNAL SUBROUTINES .. EXTERNAL INVR,STINVR,STZROR,ZROR C .. C .. STATEMENT FUNCTIONS .. DOUBLE PRECISION CUM C .. C .. STATEMENT FUNCTION DEFINITIONS .. C ..STATEMENT FUNCTIONS.. CUM(SS,XNXN,PPR,PP) = CUMNBN(SS,XNXN,PPR) - PP C .. C .. EXECUTABLE STATEMENTS .. C C CHECK ARGUMENTS C IF (.NOT. ((WHICH.LT.1).OR. (WHICH.GT.4))) GO TO 10 IF (WHICH.LT.1) THEN BOUND = 1.0D0 ELSE BOUND = 4.0D0 END IF STATUS = -1 RETURN 10 IF (WHICH.EQ.1) GO TO 30 C C P C IF (.NOT. ((P.LT.0.0D0).OR. (P.GT.ONE))) GO TO 20 IF (P.LT.0.0D0) THEN BOUND = 0.0D0 ELSE BOUND = ONE END IF STATUS = -2 RETURN 20 CONTINUE 30 IF (WHICH.EQ.2) GO TO 50 C C S C IF (.NOT. (S.LT.0.0D0)) GO TO 40 BOUND = 0.0D0 STATUS = -3 RETURN 40 CONTINUE 50 IF (WHICH.EQ.3) GO TO 70 C C XN C IF (.NOT. (XN.LT.0.0D0)) GO TO 60 BOUND = 0.0D0 STATUS = -4 RETURN 60 CONTINUE 70 IF (WHICH.EQ.4) GO TO 90 C C PR C IF (.NOT. ((PR.LT.0.0D0).OR. (PR.GT.1.0D0))) GO TO 80 IF (PR.LT.0.0D0) THEN BOUND = 0.0D0 ELSE BOUND = 1.0D0 END IF STATUS = -5 RETURN 80 CONTINUE C C CALCULATE ANSWERS C 90 IF ((1).EQ. (WHICH)) THEN C C CALCULATING P C P = CUMNBN(S,XN,PR) STATUS = 0 ELSE IF ((2).EQ. (WHICH)) THEN C C CALCULATING S C S = 5.0D0 CALL STINVR(0.0D0,INF,0.5D0,0.5D0,5.0D0,TOL,TOL) STATUS = 0 CALL INVR(STATUS,S,FX,QLEFT,QHI) 100 IF (.NOT. (STATUS.EQ.1)) GO TO 110 FX = CUM(S,XN,PR,P) CALL INVR(STATUS,S,FX,QLEFT,QHI) GO TO 100 110 IF (.NOT. (STATUS.EQ.-1)) GO TO 140 IF (.NOT. (QLEFT)) GO TO 120 STATUS = 1 BOUND = 0.0D0 GO TO 130 120 STATUS = 2 BOUND = INF 130 CONTINUE 140 CONTINUE ELSE IF ((3).EQ. (WHICH)) THEN C C CALCULATING XN C XN = 5.0D0 CALL STINVR(0.0D0,INF,0.5D0,0.5D0,5.0D0,TOL,TOL) STATUS = 0 CALL INVR(STATUS,XN,FX,QLEFT,QHI) 150 IF (.NOT. (STATUS.EQ.1)) GO TO 160 FX = CUM(S,XN,PR,P) CALL INVR(STATUS,XN,FX,QLEFT,QHI) GO TO 150 160 IF (.NOT. (STATUS.EQ.-1)) GO TO 190 IF (.NOT. (QLEFT)) GO TO 170 STATUS = 1 BOUND = 0.0D0 GO TO 180 170 STATUS = 2 BOUND = INF 180 CONTINUE 190 CONTINUE ELSE IF ((4).EQ. (WHICH)) THEN C C CALCULATING PR C CALL STZROR(0.0D0,1.0D0,TOL,TOL) STATUS = 0 CALL ZROR(STATUS,PR,FX,XLO,XHI,QLEFT,QHI) 200 IF (.NOT. (STATUS.EQ.1)) GO TO 210 FX = CUM(S,XN,PR,P) CALL ZROR(STATUS,PR,FX,XLO,XHI,QLEFT,QHI) GO TO 200 210 IF (.NOT. (STATUS.EQ.-1)) GO TO 240 IF (.NOT. (QLEFT)) GO TO 220 STATUS = 1 BOUND = 0.0D0 GO TO 230 220 STATUS = 2 BOUND = 1.0D0 230 CONTINUE 240 END IF RETURN END C----------------------------------------------------------- SUBROUTINE CDFNOR(WHICH,P,X,MEAN,SD,STATUS,BOUND) C********************************************************************** C C SUBROUTINE CDFNOR( WHICH, P, X, MEAN, SD, STATUS, BOUND ) C CUMULATIVE DISTRIBUTION FUNCTION C NORMAL DISTRIBUTION C C C FUNCTION C C C CALCULATES ANY ONE PARAMETER OF THE NORMAL C DISTRIBUTION GIVEN VALUES FOR THE OTHERS. C C C ARGUMENTS C C C WHICH --> INTEGER INDICATING WHICH OF THE NEXT PARAMETER C VALUES IS TO BE CALCULATED USING VALUES OF THE OTHERS. C LEGAL RANGE: 1..4 C INTEGER WHICH C C P <--> THE INTEGRAL FROM -INFINITY TO X OF THE NORMAL DENSITY. C INPUT RANGE: [1E-7,1-1E-7]. C DOUBLE PRECISION P C C X < --> UPPER LIMIT OF INTEGRATION OF THE NORMAL-DENSITY. C INPUT RANGE: ( -INFINITY, +INFINITY) C DOUBLE PRECISION X C C MEAN <--> THE MEAN OF THE NORMAL DENSITY. C INPUT RANGE: (-INFINITY, +INFINITY) C DOUBLE PRECISION MEAN C C SD <--> STANDARD DEVIATION OF THE NORMAL DENSITY. C INPUT RANGE: (0, +INFINITY). C DOUBLE PRECISION SD 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 THE RATIONAL FUNCTIONS FROM PAGES 90-95 OF KENNEDY AND GENTLE, C STATISTICAL COMPUTING, MARCEL DEKKER, NY, 1980 ARE USED TO C COMPUTE BOTH THE CUMULATIVE AND INVERSE STANDARD NORMAL. C THEREFORE NO SEARCHES ARE NECESSARY FOR ANY PARAMETER. C C FOR X < -15, THE ASYMPTOTIC EXPANSION FOR THE NORMAL IS USED. C THIS IS FORMULA 26.2.12 OF ABRAMOWITZ AND STEGUN. C C C NOTE C C C THE NORMAL DENSITY IS PROPORTIONAL TO C EXP( - 0.5 * (( X - MEAN)/SD)**2) C C C********************************************************************** C .. PARAMETERS .. DOUBLE PRECISION ZERO,ONE PARAMETER (ZERO=1.0D-7,ONE=1.0D0-1.0D-7) C .. C .. SCALAR ARGUMENTS .. DOUBLE PRECISION BOUND,MEAN,P,SD,X INTEGER STATUS,WHICH C .. C .. LOCAL SCALARS .. DOUBLE PRECISION Z C .. C .. EXTERNAL FUNCTIONS .. DOUBLE PRECISION CUMNOR,INVNOR EXTERNAL CUMNOR,INVNOR C .. C .. EXECUTABLE STATEMENTS .. C C CHECK ARGUMENTS C STATUS = 0 IF (.NOT. ((WHICH.LT.1).OR. (WHICH.GT.4))) GO TO 10 IF (WHICH.LT.1) THEN BOUND = 1.0D0 ELSE BOUND = 4.0D0 END IF STATUS = -1 RETURN 10 IF (WHICH.EQ.1) GO TO 30 C C P C IF (.NOT. ((P.LT.ZERO).OR. (P.GT.ONE))) GO TO 20 IF (P.LT.ZERO) THEN BOUND = ZERO ELSE BOUND = ONE END IF STATUS = -2 RETURN 20 CONTINUE 30 IF (WHICH.EQ.4) GO TO 50 C C SD C IF (.NOT. (SD.LE.0.0D0)) GO TO 40 BOUND = 0.0D0 STATUS = -5 RETURN 40 CONTINUE C C CALCULATE ANSWERS C 50 IF ((1).EQ. (WHICH)) THEN C C COMPUTING P C Z = (X-MEAN)/SD P = CUMNOR(Z) ELSE IF ((2).EQ. (WHICH)) THEN C C COMPUTING X C Z = INVNOR(P) X = SD*Z + MEAN ELSE IF ((3).EQ. (WHICH)) THEN C C COMPUTING THE MEAN C Z = INVNOR(P) MEAN = X - SD*Z ELSE IF ((4).EQ. (WHICH)) THEN C C COMPUTING SD C Z = INVNOR(P) SD = (X-MEAN)/Z END IF RETURN END C--------------------------------------------------------- SUBROUTINE CDFPOI(WHICH,P,S,XLAM,STATUS,BOUND) C********************************************************************** C C SUBROUTINE CDFPOI( WHICH, P, S, XLAM, STATUS, BOUND ) C CUMULATIVE DISTRIBUTION FUNCTION C POISSON DISTRIBUTION C C C FUNCTION C C C CALCULATES ANY ONE PARAMETER OF THE POISSON C DISTRIBUTION GIVEN VALUES FOR THE OTHERS. C C C ARGUMENTS C C C WHICH --> INTEGER INDICATING WHICH ARGUMENT C VALUE IS TO BE CALCULATED FROM THE OTHERS. C LEGAL RANGE: 1..3 C INTEGER WHICH C C P <--> THE CUMULATION FROM 0 TO S OF THE POISSON DENSITY. C INPUT RANGE: [0,1-1E-7]. C DOUBLE PRECISION P C C S <--> UPPER LIMIT OF CUMULATION OF THE POISSON. C INPUT RANGE: [0, +INFINITY). C SEARCH RANGE: [0,1E30] C DOUBLE PRECISION S C C XLAM <--> MEAN OF THE POISSON DISTRIBUTION. C INPUT RANGE: [0, +INFINITY). C SEARCH RANGE: [0,1E30] C DOUBLE PRECISION XLAM C C C METHOD C C C FORMULA 26.4.21 OF ABRAMOWITZ AND STEGUN, HANDBOOK OF C MATHEMATICAL FUNCTIONS (1966) IS USED TO REDUCE THE COMPUTATION C OF THE CUMULATIVE DISTRIBUTION FUNCTION TO THAT OF COMPUTING A C CHI-SQUARE, HENCE AN INCOMPLETE GAMMA FUNCTION. C C CUMULATIVE DISTRIBUTION FUNCTION (P) IS CALCULATED DIRECTLY. C COMPUTATION OF OTHER PARAMETERS INVOLVE A SEACH 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********************************************************************** C .. PARAMETERS .. DOUBLE PRECISION TOL PARAMETER (TOL=1.0D-4) DOUBLE PRECISION ONE,INF PARAMETER (ONE=1.0D0-1.0D-7,INF=1.0D30) C .. C .. SCALAR ARGUMENTS .. DOUBLE PRECISION BOUND,P,S,XLAM INTEGER STATUS,WHICH C .. C .. LOCAL SCALARS .. DOUBLE PRECISION FX,PP,SS,XXLAM LOGICAL QHI,QLEFT C .. C .. EXTERNAL FUNCTIONS .. DOUBLE PRECISION CUMPOI EXTERNAL CUMPOI C .. C .. EXTERNAL SUBROUTINES .. EXTERNAL INVR,STINVR C .. C .. STATEMENT FUNCTIONS .. DOUBLE PRECISION CUM C .. C .. STATEMENT FUNCTION DEFINITIONS .. CUM(SS,XXLAM,PP) = CUMPOI(SS,XXLAM) - PP C .. C .. EXECUTABLE STATEMENTS .. C C CHECK ARGUMENTS C IF (.NOT. ((WHICH.LT.1).OR. (WHICH.GT.3))) GO TO 10 IF (WHICH.LT.1) THEN BOUND = 1.0D0 ELSE BOUND = 3.0D0 END IF STATUS = -1 RETURN 10 IF (WHICH.EQ.1) GO TO 30 C C P C IF (.NOT. ((P.LT.0.0D0).OR. (P.GT.ONE))) GO TO 20 IF (P.LT.0.0D0) THEN BOUND = 0.0D0 ELSE BOUND = ONE END IF STATUS = -2 RETURN 20 CONTINUE 30 IF (WHICH.EQ.2) GO TO 50 C C S C IF (.NOT. (S.LT.0.0D0)) GO TO 40 BOUND = 0.0D0 STATUS = -3 RETURN 40 CONTINUE 50 IF (WHICH.EQ.3) GO TO 70 C C XLAM C IF (.NOT. (XLAM.LT.0.0D0)) GO TO 60 BOUND = 0.0D0 STATUS = -4 RETURN 60 CONTINUE C C CALCULATE ANSWERS C 70 IF ((1).EQ. (WHICH)) THEN C C CALCULATING P C P = CUMPOI(S,XLAM) STATUS = 0 ELSE IF ((2).EQ. (WHICH)) THEN C C CALCULATING S C S = 5.0D0 CALL STINVR(0.0D0,INF,0.5D0,0.5D0,5.0D0,TOL,TOL) STATUS = 0 CALL INVR(STATUS,S,FX,QLEFT,QHI) 80 IF (.NOT. (STATUS.EQ.1)) GO TO 90 FX = CUM(S,XLAM,P) CALL INVR(STATUS,S,FX,QLEFT,QHI) GO TO 80 90 IF (.NOT. (STATUS.EQ.-1)) GO TO 120 IF (.NOT. (QLEFT)) GO TO 100 STATUS = 1 BOUND = 0.0D0 GO TO 110 100 STATUS = 2 BOUND = INF 110 CONTINUE 120 CONTINUE ELSE IF ((3).EQ. (WHICH)) THEN C C CALCULATING XLAM C XLAM = 5.0D0 CALL STINVR(0.0D0,INF,0.5D0,0.5D0,5.0D0,TOL,TOL) STATUS = 0 CALL INVR(STATUS,XLAM,FX,QLEFT,QHI) 130 IF (.NOT. (STATUS.EQ.1)) GO TO 140 FX = CUM(S,XLAM,P) CALL INVR(STATUS,XLAM,FX,QLEFT,QHI) GO TO 130 140 IF (.NOT. (STATUS.EQ.-1)) GO TO 170 IF (.NOT. (QLEFT)) GO TO 150 STATUS = 1 BOUND = 0.0D0 GO TO 160 150 STATUS = 2 BOUND = INF 160 CONTINUE 170 END IF RETURN END C------------------------------------------------------------ SUBROUTINE CDFT(WHICH,P,T,DF,STATUS,BOUND) C********************************************************************** C C SUBROUTINE CDFT( WHICH, P, T, DF, STATUS, BOUND ) C CUMULATIVE DISTRIBUTION FUNCTION C T DISTRIBUTION C C C FUNCTION C C C CALCULATES ANY ONE PARAMETER OF THE T DISTRIBUTION GIVEN C VALUES FOR THE OTHERS. C C C ARGUMENTS C C C WHICH --> INTEGER INDICATING WHICH ARGUMENT C VALUES IS TO BE CALCULATED FROM THE OTHERS. C LEGAL RANGE: 1..3 C INTEGER WHICH C C P <--> THE INTEGRAL FROM -INFINITY TO T OF THE T-DENSITY. C INPUT RANGE: [1E-7,1-1E-7]. C DOUBLE PRECISION P C C T <--> UPPER LIMIT OF INTEGRATION OF THE T-DENSITY. C INPUT RANGE: ( -INFINITY, +INFINITY). C SEARCH RANGE: [ -1E30, 1E30 ] C DOUBLE PRECISION T C C DF <--> DEGREES OF FREEDOM OF THE T-DISTRIBUTION. C INPUT RANGE: (0 , +INFINITY). C SEARCH RANGE: [1E-30, 1E30] C DOUBLE PRECISION DF 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.5.27 OF ABRAMOWITZ AND STEGUN, HANDBOOK OF C MATHEMATICAL FUNCTIONS (1966) IS USED TO REDUCE THE COMPUTATION C OF THE CUMULATIVE DISTRIBUTION FUNCTION TO THAT OF AN INCOMPLETE C BETA. C C COMPUTATION OF OTHER PARAMETERS INVOLVE A SEACH 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 .. PARAMETERS .. DOUBLE PRECISION TOL PARAMETER (TOL=1.0D-4) DOUBLE PRECISION ZERO,ONE,INF PARAMETER (ZERO=1.0D-30,ONE=1.0D0-1.0D-7,INF=1.0D30) C .. C .. SCALAR ARGUMENTS .. DOUBLE PRECISION BOUND,DF,P,T INTEGER STATUS,WHICH C .. C .. LOCAL SCALARS .. DOUBLE PRECISION DFDF,FX,PP,TT LOGICAL QHI,QLEFT C .. C .. EXTERNAL FUNCTIONS .. DOUBLE PRECISION CUMT,T1 EXTERNAL CUMT,T1 C .. C .. EXTERNAL SUBROUTINES .. EXTERNAL INVR,STINVR C .. C .. STATEMENT FUNCTIONS .. DOUBLE PRECISION CUM C .. C .. STATEMENT FUNCTION DEFINITIONS .. CUM(TT,DFDF,PP) = CUMT(TT,DFDF) - PP C .. C .. EXECUTABLE STATEMENTS .. C C CHECK ARGUMENTS C IF (.NOT. ((WHICH.LT.1).OR. (WHICH.GT.3))) GO TO 10 IF (WHICH.LT.1) THEN BOUND = 1.0D0 ELSE BOUND = 3.0D0 END IF STATUS = -1 RETURN 10 IF (WHICH.EQ.1) GO TO 30 C C P C IF (.NOT. ((P.LT.ZERO).OR. (P.GT.ONE))) GO TO 20 IF (P.LT.ZERO) THEN BOUND = ZERO ELSE BOUND = ONE END IF STATUS = -2 RETURN 20 CONTINUE 30 IF (WHICH.EQ.3) GO TO 50 C C DF C IF (.NOT. (DF.LE.0.0D0)) GO TO 40 BOUND = 0.0D0 STATUS = -4 RETURN 40 CONTINUE C C CALCULATE ANSWERS C 50 IF ((1).EQ. (WHICH)) THEN C C COMPUTING P C P = CUMT(T,DF) STATUS = 0 ELSE IF ((2).EQ. (WHICH)) THEN C C COMPUTING T C C .. GET INITIAL APPROXIMATION FOR T C T = T1(P,DF) CALL STINVR(-INF,INF,0.5D0,0.5D0,5.0D0,TOL,TOL) STATUS = 0 CALL INVR(STATUS,T,FX,QLEFT,QHI) 60 IF (.NOT. (STATUS.EQ.1)) GO TO 70 FX = CUM(T,DF,P) CALL INVR(STATUS,T,FX,QLEFT,QHI) GO TO 60 70 IF (.NOT. (STATUS.EQ.-1)) GO TO 100 IF (.NOT. (QLEFT)) GO TO 80 STATUS = 1 BOUND = -INF GO TO 90 80 STATUS = 2 BOUND = INF 90 CONTINUE 100 CONTINUE ELSE IF ((3).EQ. (WHICH)) THEN C C COMPUTING DF C DF = 5.0D0 CALL STINVR(ZERO,INF,0.5D0,0.5D0,5.0D0,TOL,TOL) STATUS = 0 CALL INVR(STATUS,DF,FX,QLEFT,QHI) 110 IF (.NOT. (STATUS.EQ.1)) GO TO 120 FX = CUM(T,DF,P) CALL INVR(STATUS,DF,FX,QLEFT,QHI) GO TO 110 120 IF (.NOT. (STATUS.EQ.-1)) GO TO 150 IF (.NOT. (QLEFT)) GO TO 130 STATUS = 1 BOUND = ZERO GO TO 140 130 STATUS = 2 BOUND = INF 140 CONTINUE 150 END IF RETURN END C--------------------------------------------------- DOUBLE PRECISION FUNCTION CUMBET(XX,XA,XB) C********************************************************************** C C DOUBLE PRECISION FUNCTION CUMBET(X,A,B) C CUMULATIVE INCOMPLETE BETA DISTRIBUTION C C C FUNCTION C C C RETURNS THE CDF TO X OF THE INCOMPLETE BETA DISTRIBUTION C WITH PARAMETERS A AND B. THIS IS THE INTEGRAL FROM 0 TO X C OF (1/B(A,B))*F(T)) WHERE F(T) = T**(A-1) * (1-T)**(B-1) C C C ARGUMENTS C C C X --> UPPER LIMIT OF INTEGRATION. C X IS DOUBLE PRECISION C C A --> FIRST PARAMETER OF THE BETA DISTRIBUTION. C A IS DOUBLE PRECISION C C B --> SECOND PARAMETER OF THE BETA DISTRIBUTION. C B IS DOUBLE PRECISION C C C METHOD C C C C CALLS THE DOUBLE PRECISION VERSION OF THIS ROUTINE DUMBET. C C********************************************************************** C .. C .. SCALAR ARGUMENTS .. DOUBLE PRECISION XA,XB,XX C .. C .. EXTERNAL FUNCTIONS .. DOUBLE PRECISION DUMBET EXTERNAL DUMBET C .. C .. EXECUTABLE STATEMENTS .. CUMBET = DUMBET(XX,XA,XB) RETURN END C-------------------------------------------------- DOUBLE PRECISION FUNCTION CUMBIN(S,XN,PR) C C********************************************************************** C C DOUBLE PRECISION FUNCTION CUMBIN(S,XN,PBIN) C CUMULATIVE BINOMIAL DISTRIBUTION C C C FUNCTION C C C RETURNS THE PROBABILITY OF 0 TO S SUCCESSES IN XN BINOMIAL C TRIALS, EACH OF WHICH HAS A PROBABILITY OF SUCCESS, PBIN. C C C ARGUMENTS C C C S --> THE UPPER LIMIT OF CUMULATION OF THE BINOMIAL DISTRIBUTION. C S IS DOUBLE PRECISION C C XN --> THE NUMBER OF BINOMIAL TRIALS. C XN IS DOUBLE PRECISION C C PBIN --> THE PROBABILITY OF SUCCESS IN EACH BINOMIAL TRIAL. C PBIN IS DOUBLE PRECISION C C C METHOD C C C FORMULA 26.5.24 OF ABRAMOWITZ AND STEGUN, HANDBOOK OF C MATHEMATICAL FUNCTIONS (1966) IS USED TO REDUCE THE BINOMIAL C DISTRIBUTION TO THE CUMULATIVE BETA DISTRIBUTION. C C********************************************************************** C .. SCALAR ARGUMENTS .. DOUBLE PRECISION PR,S,XN C .. C .. EXTERNAL FUNCTIONS .. DOUBLE PRECISION DUMBET EXTERNAL DUMBET C .. C .. EXECUTABLE STATEMENTS .. IF (.NOT. (S.LT.XN)) GO TO 10 CUMBIN =1.D0-DUMBET(PR,S+1.D0,XN-S) GO TO 20 10 CUMBIN = 1.0D0 20 RETURN END C----------------------------------------------- DOUBLE PRECISION FUNCTION CUMCHI(X,DF) C********************************************************************** C C DOUBLE PRECISION FUNCTION CUMCHI(X,DF) C CUMULATIVE OF THE CHI-SQUARE DISTRIBUTION C C C FUNCTION C C C CALCULATES THE CUMULATIVE CHI-SQUARE DISTRIBUTION. C C C ARGUMENTS C C C X --> UPPER LIMIT OF INTEGRATION OF THE C CHI-SQUARE DISTRIBUTION. C X IS DOUBLE PRECISION C C DF --> DEGREES OF FREEDOM OF THE C CHI-SQUARE DISTRIBUTION. C DF IS DOUBLE PRECISION C C C METHOD C C C CALLS INCOMPLETE GAMMA FUNCTION (CUMGAM) C C********************************************************************** C .. SCALAR ARGUMENTS .. DOUBLE PRECISION DF,X C .. C .. LOCAL SCALARS .. DOUBLE PRECISION A,XX C .. C .. EXTERNAL FUNCTIONS .. DOUBLE PRECISION CUMGAM EXTERNAL CUMGAM C .. C .. EXECUTABLE STATEMENTS .. A = DF*0.5 XX = X*0.5 CUMCHI = CUMGAM(XX,A) RETURN END C-------------------------------------------------- DOUBLE PRECISION FUNCTION CUMCHN(X,DF,PNONC) C C********************************************************************** C C DOUBLE PRECISION FUNCTION CUMCHN(X,DF,PNONC) C CUMULATIVE OF THE NON-CENTRAL CHI-SQUARE DISTRIBUTION C C C FUNCTION C C C CALCULATES THE CUMULATIVE NON-CENTRAL CHI-SQUARE C DISTRIBUTION, I.E., THE PROBABILITY THAT A RANDOM VARIABLE C WHICH FOLLOWS THE NON-CENTRAL CHI-SQUARE DISTRIBUTION, WITH C NON-CENTRALITY PARAMETER PNONC AND CONTINUOUS DEGREES OF C FREEDOM DF, IS LESS THAN OR EQUAL TO X. C C C ARGUMENTS C C C X --> UPPER LIMIT OF INTEGRATION OF THE NON-CENTRAL C CHI-SQUARE DISTRIBUTION. C X IS DOUBLE PRECISION C C DF --> DEGREES OF FREEDOM OF THE NON-CENTRAL C CHI-SQUARE DISTRIBUTION. C DF IS DOUBLE PRECISION C C PNONC --> NON-CENTRALITY PARAMETER OF THE NON-CENTRAL C CHI-SQUARE DISTRIBUTION. C PNONC IS DOUBLE PRECISION C C C METHOD C C C USES FORMULA 26.4.25 OF ABRAMOWITZ AND STEGUN, HANDBOOK OF C MATHEMATICAL FUNCTIONS, US NBS (1966) TO CALCULATE THE C NON-CENTRAL CHI-SQUARE. C C C VARIABLES C C C EPS --- CONVERGENCE CRITERION. THE SUM STOPS WHEN A C TERM IS LESS THAN EPS*SUM. C EPS IS DOUBLE PRECISION C C NTIRED --- MAXIMUM NUMBER OF TERMS TO BE EVALUATED C IN EACH SUM. C NTIRED IS INTEGER C C QCONV --- .TRUE. IF CONVERGENCE ACHIEVED - C I.E., PROGRAM DID NOT STOP ON NTIRED CRITERION. C QCONV IS LOGICAL C C********************************************************************** C C C .. SCALAR ARGUMENTS .. DOUBLE PRECISION DF,PNONC,X C .. C .. LOCAL SCALARS .. DOUBLE PRECISION ADJ,CENTAJ,CENTWT,CHID2,DFD2,EPS,LCNTAJ, +LCNTWT,LFACT,PCENT,PTERM,SUM,SUMADJ,TERM,WT,XNONC,XX INTEGER I,ICENT,ITERB,ITERF,NTIRED C .. C .. EXTERNAL FUNCTIONS .. DOUBLE PRECISION ALNGAM,CUMCHI EXTERNAL ALNGAM,CUMCHI C .. C .. INTRINSIC FUNCTIONS .. INTRINSIC ALOG,EXP,FLOAT C .. C .. STATEMENT FUNCTIONS .. DOUBLE PRECISION DG LOGICAL QSMALL,QTIRED C .. C .. DATA STATEMENTS .. DATA NTIRED,EPS/1000D0,1.0D-5/ C .. C .. STATEMENT FUNCTION DEFINITIONS .. QTIRED(I) = I .GT. NTIRED QSMALL(XX) = SUM .LT. 1.0D-20 .OR. XX .LT. EPS*SUM DG(I) = DF + 2.0D0*FLOAT(I) C .. C .. EXECUTABLE STATEMENTS .. C C C C WHEN NON-CENTRALITY PARAMETER IS (ESSENTIALLY) ZERO, C USE CUMULATIVE CHI-SQUARE DISTRIBUTION C C IF (.NOT. (PNONC.LE.1.0D-10)) GO TO 10 CUMCHN = CUMCHI(X,DF) RETURN 10 XNONC = PNONC/2.0D0 C********************************************************************** C C THE FOLLOWING CODE CALCUALTES THE WEIGHT, CHI-SQUARE, AND C ADJUSTMENT TERM FOR THE CENTRAL TERM IN THE INFINITE SERIES. C THE CENTRAL TERM IS THE ONE IN WHICH THE POISSON WEIGHT IS C GREATEST. THE ADJUSTMENT TERM IS THE AMOUNT THAT MUST C BE SUBTRACTED FROM THE CHI-SQUARE TO MOVE UP TWO DEGREES C OF FREEDOM. C C********************************************************************** ICENT = XNONC IF (ICENT.EQ.0) ICENT = 1 CHID2 = X/2.0D0 C C C CALCULATE CENTRAL WEIGHT TERM C C LFACT = ALNGAM(FLOAT(ICENT+1)) LCNTWT = -XNONC + ICENT*DLOG(XNONC) - LFACT CENTWT = DEXP(LCNTWT) C C C CALCULATE CENTRAL CHI-SQUARE C C PCENT = CUMCHI(X,DG(ICENT)) C C C CALCULATE CENTRAL ADJUSTMENT TERM C C DFD2 = DG(ICENT)/2.0D0 LFACT = ALNGAM(1.0D0+DFD2) LCNTAJ = DFD2*DLOG(CHID2) - CHID2 - LFACT CENTAJ = DEXP(LCNTAJ) SUM = CENTWT*PCENT C********************************************************************** C C SUM BACKWARDS FROM THE CENTRAL TERM TOWARDS ZERO. C QUIT WHENEVER EITHER C (1) THE ZERO TERM IS REACHED, OR C (2) THE TERM GETS SMALL RELATIVE TO THE SUM, OR C (3) MORE THAN NTIRED TERMS ARE TOTALED. C C********************************************************************** ITERB = 0 SUMADJ = 0.0D0 ADJ = CENTAJ WT = CENTWT I = ICENT C GO TO 30 20 IF (QTIRED(ITERB) .OR. QSMALL(TERM) .OR. I.EQ.0) GO TO 40 30 DFD2 = DG(I)/2.0D0 C C C ADJUST CHI-SQUARE FOR TWO FEWER DEGREES OF FREEDOM. C THE ADJUSTED VALUE ENDS UP IN PTERM. C C ADJ = ADJ*DFD2/CHID2 SUMADJ = SUMADJ + ADJ PTERM = PCENT + SUMADJ C C C ADJUST POISSON WEIGHT FOR J DECREASED BY ONE C C WT = WT* (I/XNONC) TERM = WT*PTERM SUM = SUM + TERM I = I - 1 ITERB = ITERB + 1 GO TO 20 40 ITERF = 0 C********************************************************************** C C NOW SUM FORWARD FROM THE CENTRAL TERM TOWARDS INFINITY. C QUIT WHEN EITHER C (1) THE TERM GETS SMALL RELATIVE TO THE SUM, OR C (2) MORE THAN NTIRED TERMS ARE TOTALED. C C********************************************************************** SUMADJ = CENTAJ ADJ = CENTAJ WT = CENTWT I = ICENT C GO TO 60 50 IF (QTIRED(ITERF) .OR. QSMALL(TERM)) GO TO 70 C C C UPDATE WEIGHTS FOR NEXT HIGHER J C C 60 WT = WT* (XNONC/ (I+1)) C C C CALCULATE PTERM AND ADD TERM TO SUM C C PTERM = PCENT - SUMADJ TERM = WT*PTERM SUM = SUM + TERM C C C UPDATE ADJUSTMENT TERM FOR DF FOR NEXT ITERATION C C I = I + 1 DFD2 = DG(I)/2.0D0 ADJ = ADJ*CHID2/DFD2 SUMADJ = SUMADJ + ADJ ITERF = ITERF + 1 GO TO 50 70 CUMCHN = SUM C RETURN END C-------------------------------------------------------------------- DOUBLE PRECISION FUNCTION CUMF(F,DFN,DFD) C C********************************************************************** C C DOUBLE PRECISION FUNCTION CUMF(F,DFN,DFD) C CUMULATIVE F DISTRIBUTION C C C FUNCTION C C C COMPUTES THE INTEGRAL FROM 0 TO F OF THE F-DENSITY WITH DFN C AND DFD DEGREES OF FREEDOM. C C C ARGUMENTS C C C F --> UPPER LIMIT OF INTEGRATION OF THE F-DENSITY. C F IS DOUBLE PRECISION C C DFN --> DEGREES OF FREEDOM OF THE NUMERATOR SUM OF SQUARES. C DFN IS DOUBLE PRECISION C C DFD --> DEGREES OF FREEDOM OF THE DENOMINATOR SUM OF SQUARES. C DFD IS DOUBLE PRECISION C C C METHOD C C C FORMULA 26.5.28 OF ABRAMOWITZ AND STEGUN IS USED TO REDUCE C THE CUMULATIVE F TO A CUMULATIVE BETA DISTRIBUTION. C C C NOTE C C C IF F IS LESS THAN OR EQUAL TO 0, 0 IS RETURNED. C C********************************************************************** C .. SCALAR ARGUMENTS .. DOUBLE PRECISION DFD,DFN,F C .. C .. LOCAL SCALARS .. DOUBLE PRECISION DDFD,DDFN,DF,DSUM,DUMF,DUMMY,PROD,XX,YY INTEGER IERR C .. C .. INTRINSIC FUNCTIONS .. INTRINSIC DBLE,REAL C .. C .. PARAMETERS .. DOUBLE PRECISION HALF PARAMETER (HALF=0.5D0) DOUBLE PRECISION DONE PARAMETER (DONE=1.0D0) C .. C .. EXTERNAL SUBROUTINES .. EXTERNAL BRATIO C .. C .. EXECUTABLE STATEMENTS .. IF (.NOT. (F.LE.0.0D0)) GO TO 10 CUMF = 0.0D0 RETURN 10 DDFN = DFN C C XX IS SUCH THAT THE INCOMPLETE BETA WITH PARAMETERS C DFD/2 AND DFN/2 EVALUATED AT XX IS 1 - CUMF C C YY IS 1 - XX C C CALCULATE THE SMALLER OF XX AND YY ACCURATELY C DDFD = DFD DF = F PROD = DDFN*DF DSUM = DDFD + PROD XX = DDFD/DSUM IF (XX.GT.HALF) THEN YY = PROD/DSUM XX = DONE - YY ELSE YY = DONE - XX END IF CALL BRATIO(DDFD*HALF,DDFN*HALF,XX,YY,DUMMY,DUMF,IERR) CUMF =DUMF RETURN END C------------------------------------------------------------------- DOUBLE PRECISION FUNCTION CUMFNC(X,DFN,DFD,PNONC) C********************************************************************** C C F -NON- -C-ENTRAL F DISTRIBUTION C C C C FUNCTION C C C COMPUTES NONCENTRAL F DISTRIBUTION WITH DFN AND DFD C DEGREES OF FREEDOM AND NONCENTRALITY PARAMETER PNONC C C C ARGUMENTS C C C X --> UPPER LIMIT OF INTEGRATION OF NONCENTRAL F IN EQUATION C C DFN --> DEGREES OF FREEDOM OF NUMERATOR C C DFD --> DEGREES OF FREEDOM OF DENOMINATOR C C PNONC --> NONCENTRALITY PARAMETER. C C C METHOD C C C USES FORMULA 26.6.20 OF REFERENCE FOR INFINITE SERIES. C SERIES IS CALCULATED BACKWARD AND FORWARD FROM J = LAMBDA/2 C (THIS IS THE TERM WITH THE LARGEST POISSON WEIGHT) UNTIL C THE CONVERGENCE CRITERION IS MET. C C FOR SPEED, THE INCOMPLETE BETA FUNCTIONS ARE EVALUATED C BY FORMULA 26.5.16. C C C REFERENCE C C C HANDBOOD OF MATHEMATICAL FUNCTIONS C EDITED BY MILTON ABRAMOWITZ AND IRENE A. STEGUN C NATIONAL BUREAU OF STANDARDS APPLIED MATEMATICS SERIES - 55 C MARCH 1965 C P 947, EQUATIONS 26.6.17, 26.6.18 C C C NOTE C C C THE SUM CONTINUES UNTIL A SUCCEEDING TERM IS LESS THAN EPS C TIMES THE SUM (OR THE SUM IS LESS THAN 1.0E-20). EPS IS C SET TO 1.0E-4 IN A DATA STATEMENT WHICH CAN BE CHANGED. C C********************************************************************** C .. SCALAR ARGUMENTS .. DOUBLE PRECISION DFD,DFN,PNONC,X C .. C .. LOCAL SCALARS .. DOUBLE PRECISION DBETDN,DDFD,DDFN,DF,DSUM,DUMMY,PROD,XX,YY DOUBLE PRECISION ADN,ARG,AUP,B,BETDN,BETUP,CENTWT,DNTERM,EPS, +OMARG,SUM,UPTERM, XMULT,XNONC INTEGER I,ICENT,IERR C .. C .. EXTERNAL FUNCTIONS .. DOUBLE PRECISION ALNGAM,CUMF EXTERNAL ALNGAM,CUMF C .. C .. INTRINSIC FUNCTIONS .. INTRINSIC DLOG,DEXP,DBLE C .. C .. STATEMENT FUNCTIONS .. LOGICAL QSMALL C .. C .. EXTERNAL SUBROUTINES .. EXTERNAL BRATIO C .. C .. PARAMETERS .. DOUBLE PRECISION HALF PARAMETER (HALF=0.5D0) DOUBLE PRECISION DONE PARAMETER (DONE=1.0D0) C .. C .. DATA STATEMENTS .. DATA EPS/1.0D-4/ C .. C .. STATEMENT FUNCTION DEFINITIONS .. QSMALL(X) = SUM .LT. 1.0D-20 .OR. X .LT. EPS*SUM C .. C .. EXECUTABLE STATEMENTS .. C IF (.NOT. (X.LE.0.0D0)) GO TO 10 CUMFNC = 0.0D0 RETURN 10 IF (.NOT. (PNONC.LT.1.0D0-10)) GO TO 20 C C HANDLE CASE IN WHICH THE NON-CENTRALITY PARAMETER IS C (ESSENTIALLY) ZERO. CUMFNC = CUMF(X,DFN,DFD) RETURN 20 XNONC = PNONC/2.0D0 C CALCULATE THE CENTRAL TERM OF THE POISSON WEIGHTING FACTOR. ICENT = XNONC IF (ICENT.EQ.0) ICENT = 1 C COMPUTE CENTRAL WEIGHT TERM CENTWT = DEXP(-XNONC+ICENT*DLOG(XNONC)-ALNGAM(DBLE(ICENT+1))) C COMPUTE CENTRAL INCOMPLETE BETA TERM C ASSURE THAT MINIMUM OF ARG TO BETA AND 1 - ARG IS COMPUTED C ACCURATELY. DDFN = DFN DDFD = DFD DF = X PROD = DDFN*DF DSUM = DDFD + PROD YY = DDFD/DSUM IF (YY.GT.HALF) THEN XX = PROD/DSUM YY = DONE - XX ELSE XX = DONE - YY END IF CALL BRATIO(DDFN*HALF+DBLE(ICENT),DDFD*HALF,XX,YY,DBETDN,DUMMY, + IERR) BETDN =DBETDN ARG = XX OMARG =YY ADN = DFN/2.0D0 + ICENT AUP = ADN B = DFD/2.0D0 BETUP = BETDN SUM = CENTWT*BETDN C NOW SUM TERMS BACKWARD FROM ICENT UNTIL CONVERGENCE OR ALL DONE XMULT = CENTWT I = ICENT DNTERM = DEXP(ALNGAM(ADN+B)-ALNGAM(ADN+1.0D0)-ALNGAM(B)+ + ADN*DLOG(ARG)+B*DLOG(OMARG)) 30 IF (QSMALL(XMULT*BETDN) .OR. I.LE.0) GO TO 40 XMULT = XMULT* (I/XNONC) I = I - 1 ADN = ADN - 1 DNTERM = (ADN+1)/ ((ADN+B)*ARG)*DNTERM BETDN = BETDN + DNTERM SUM = SUM + XMULT*BETDN GO TO 30 40 I = ICENT + 1 C NOW SUM FORWARDS UNTIL CONVERGENCE XMULT = CENTWT IF ((AUP-1D0+B).EQ.0D0) THEN UPTERM = DEXP(-ALNGAM(AUP)-ALNGAM(B)+ (AUP-1D0) +*DLOG(ARG)+B*DLOG(OMARG)) ELSE UPTERM =DEXP(ALNGAM(AUP-1D0+B)-ALNGAM(AUP)-ALNGAM(B)+ + (AUP-1D0)*DLOG(ARG)+B*DLOG(OMARG)) END IF GO TO 60 50 IF (QSMALL(XMULT*BETUP)) GO TO 70 60 XMULT = XMULT* (XNONC/I) I = I + 1 AUP = AUP + 1 UPTERM = (AUP+B-2.0D0)*ARG/ (AUP-1D0)*UPTERM BETUP = BETUP - UPTERM SUM = SUM + XMULT*BETUP GO TO 50 70 CUMFNC = SUM RETURN END C----------------------------------------------------------- DOUBLE PRECISION FUNCTION CUMGAM(XX,AA) C C********************************************************************** C C DOUBLE PRECISION FUNCTION CUMGAM(X,A) C CUMULATIVE INCOMPLETE GAMMA DISTRIBUTION C C C FUNCTION C C C COMPUTES THE CUMULATIVE OF THE INCOMPLETE GAMMA C DISTRIBUTION, I.E., THE INTEGRAL FROM 0 TO X OF C (1/GAM(A))*EXP(-T)*T**(A-1) DT C WHERE GAM(A) IS THE COMPLETE GAMMA FUNCTION OF A, I.E., C GAM(A) = INTEGRAL FROM 0 TO INFINITY OF C EXP(-T)*T**(A-1) DT C C C ARGUMENTS C C C X --> THE UPPER LIMIT OF INTEGRATION OF THE INCOMPLETE GAMMA. C X IS DOUBLE PRECISION C C A --> THE SHAPE PARAMETER OF THE INCOMPLETE GAMMA. C A IS DOUBLE PRECISION C C C METHOD C C C CALLS THE ROUTINE GRATIO. C C********************************************************************** C C .. C .. SCALAR ARGUMENTS .. DOUBLE PRECISION AA,XX C .. C .. LOCAL SCALARS .. DOUBLE PRECISION ANS,QANS C .. C .. EXTERNAL ROUTINES .. EXTERNAL GRATIO C .. C .. EXECUTABLE STATEMENTS .. IF (.NOT. (XX.LE.0.0D0)) GO TO 10 CUMGAM = 0.0D0 RETURN C CALL GRATIO ROUTINE 10 CALL GRATIO(AA,XX,ANS,QANS,1) CUMGAM =ANS RETURN END C----------------------------------------------------------- DOUBLE PRECISION FUNCTION CUMHYP(K,N,M,L) C********************************************************************** C C DOUBLE PRECISION FUNCTION CUMHYP(K,N,M,L) C CUMULATIVE HYPERGEOMETRIC DISTRIBUTION C C C FUNCTION C C C CALCULATES THE CUMULATIVE HYPERGEOMETRIC DISTRIBUTION, I.E., C THE PROBABILITY THAT A RANDOM VARIABLE WHICH FOLLOWS THE C HYPERGEOMETRIC DISTRIBUTION WITH SAMPLE SIZE N, LOT SIZE L C AND M DEFECTIVES IN THE LOT, IS LESS THAN OR EQUAL TO K. C C C ARGUMENTS C C C K --> ARGUMENT FOR WHICH HYPERGEOMETRIC DISTRIBUTION C IS TO BE EVALUATED. C K IS INTEGER C C N --> SAMPLE SIZE. C N IS INTEGER C C M --> NUMBER OF DEFECTIVES IN THE LOT. C M IS INTEGER C C L --> LOT SIZE (MUST BE GREATER THAN OR EQUAL TO C BOTH N AND M). C L IS INTEGER C C C METHOD C C C PROBABILITY THAT A HYPERGEOMETRIC RANDOM VARIABLE X HAS THE C VALUE J+1 IS EQUAL TO THE PROBABILITY THAT X HAS THE VALUE J C TIMES THE TERM C C ((M-J) * (N-J)) / ((J+1) * (L-M-N+J+1)) C C THE PROBABILITY THAT THE RANDOM VARIABLE X IS EQUAL TO I IS C CALCULATED CONVENTIONALLY AS C C C(M,I) * C(L-M,N-I) / C(L,N) C C WHERE C(A,B) IS THE NUMBER OF COMBINATIONS OF A OBJECTS C TAKEN B AT A TIME AND I IS THE MINIMUM VALUE OF J SUCH THAT C BOTH BINOMIAL COEFFICIENTS ARE DEFINED IN THE NUMERATOR. C C THE CONTRIBUTIONS OF THE REMAINING VALUES OF J ARE C ACCUMULATED BY MAKING USE OF THE PROPERTY DESCRIBED ABOVE C FOR J=I+1,K. C C C VARIABLES C C C BC1 --- TEMPORARY VARIABLE USED TO HOLD THE LOG OF ONE OF THE C BINOMIAL COEFFICIENTS USED TO CALCULATE THE FIRST C TERM OF THE CUMULATIVE HYPERGEOMETRIC. C BC1 IS DOUBLE PRECISION C C BC2 --- TEMPORARY VARIABLE USED TO HOLD THE LOG OF ONE OF THE C BINOMIAL COEFFICIENTS USED TO CALCULATE THE FIRST C TERM OF THE CUMULATIVE HYPERGEOMETRIC. C BC2 IS DOUBLE PRECISION C C BC3 --- TEMPORARY VARIABLE USED TO HOLD THE LOG OF ONE OF THE C BINOMIAL COEFFICIENTS USED TO CALCULATE THE FIRST C TERM OF THE CUMULATIVE HYPERGEOMETRIC. C BC3 IS DOUBLE PRECISION C C T --- TEMPORARY VARIABLE USED TO HOLD CUMULATIVE C HYPERGEOMETRIC AT EACH TERM. C T IS DOUBLE PRECISION C C J --- VALUE OF HYPERGEOMETRIC RANDOM VARIABLE AT EACH C ITERATION. C J IS INTEGER C C I --- MINIMUM VALUE OF J SUCH THAT BOTH BINOMIAL C COEFFICIENTS ARE DEFINED IN THE NUMERATOR OF THE C EXPRESSION FOR THE FIRST TERM IN THE HYPERGEOMETRIC C CALCULATION. I IS THE MAXIMUM OF 0 AND M-L+N. C I IS INTEGER C C MAXJ --- MAXIMUM VALUE OF J IN CALCULATION. IF MAXJ IS C GREATER THAN THE MINIMUM OF M AND N, THIS FUNCTION C RETURNS A PROBABILITY OF 1. C MAXJ IS INTEGER C C LMN --- L - M - N C USED TO AVOID RECALCULATION IN LOOP. C LMN IS INTEGER C C C FUNCTIONS AND SUBROUTINES C C C LBICO --- FUNCTION THAT CALCULATES THE LOG OF THE BINOMIAL C COEFFICIENT FOR N OUTCOMES AND S SUCCESSES. C LBICO IS DOUBLE PRECISION C C********************************************************************** C .. SCALAR ARGUMENTS .. INTEGER K,L,M,N C .. C .. LOCAL SCALARS .. DOUBLE PRECISION BC1,BC2,BC3,CUM,T INTEGER I,I99993,I99996,J,LMN,MAXJ C .. C .. EXTERNAL FUNCTIONS .. DOUBLE PRECISION LBICO EXTERNAL LBICO C .. C .. INTRINSIC FUNCTIONS .. INTRINSIC DEXP,MAX,MIN,DBLE C .. C .. EXECUTABLE STATEMENTS .. C C CALCULATE MINIMUM AND MAXIMUM VALUES WHICH RANDOM VARIABLE C MAY TAKE. C I = MAX(0,M-L+N) MAXJ = MIN(N,M) C C IF K IS OUT OF VALID BOUNDS FOR J, RETURN 0 OR 1. C OTHERWISE, PROCEED WITH CALCULATION. C IF (.NOT. (K.LT.I)) GO TO 10 STOP 'K LESS THAN I IN CUMHYP' GO TO 40 10 IF (.NOT. (K.GT.MAXJ)) GO TO 20 STOP 'K GREATER THAN MAXJ IN CUMHYP' GO TO 40 C CALCULATE-HYPERGEOMETRIC 20 ASSIGN 30 TO I99996 GO TO 50 30 CONTINUE 40 RETURN C C C STOP '*** EXECUTION FLOWING INTO FLECS PROCEDURES ***' C TO CALCULATE-HYPERGEOMETRIC C********************************************************************** C C CALCULATE I'TH TERM BY METHOD OF BINOMIAL COEFFICIENTS, C THEN CALCULATE CONTRIBUTIONS OF REMAINING TERMS RECURSIVELY C BY MULTIPLICATION. C C********************************************************************** 50 BC1 = LBICO(M,I) BC2 = LBICO(L-M,N-I) BC3 = LBICO(L,N) C T = DEXP(BC1+ (BC2-BC3)) CUM = T C IF (T.LT.0.0D0) GO TO 70 C CALC-CONTRIB-REMAIN-TERMS ASSIGN 60 TO I99993 GO TO 80 C 60 CONTINUE 70 CUMHYP = CUM GO TO I99996 C TO CALC-CONTRIB-REMAIN-TERMS C C C********************************************************************** C C CALCULATE CONTRIBUTIONS OF REMAINING TERMS RECURSIVELY C BY MULTIPLICATION. C C********************************************************************** 80 LMN = L - M - N C J = I + 1 GO TO 100 90 J = J + 1 100 IF (J.GT.K) GO TO 110 T = T*DBLE(((M-J+1)* (N-J+1)))/DBLE((J* (LMN+J))) CUM = CUM + T GO TO 90 110 GO TO I99993 C END C-------------------------------------------------------------- DOUBLE PRECISION FUNCTION CUMNBN(S,XN,PR) C********************************************************************** C C DOUBLE PRECISION FUNCTION CUMNNBN(S,XN,PR) C CUMULATIVE NEGATIVE BINOMIAL DISTRIBUTION C C C FUNCTION C C C RETURNS THE PROBABILITY THAT IT THERE WILL BE S OR FEWER FAILURES C BEFORE THERE ARE XN SUCCESSES, WITH EACH BINOMIAL TRIAL HAVING C A PROBABILITY OF SUCCESS PR. C C PROB(# FAILURES = S | XN SUCCESSES, PR) = C ( XN + S - 1 ) C ( ) * PR^XN * (1-PR)^S C ( S ) C C C ARGUMENTS C C C S --> THE NUMBER OF FAILURES C S IS DOUBLE PRECISION C C XN --> THE NUMBER OF SUCCESSES C XN IS DOUBLE PRECISION C C PR --> THE PROBABILITY OF SUCCESS IN EACH BINOMIAL TRIAL. C PR IS DOUBLE PRECISION C C C METHOD C C C FORMULA 26.5.26 OF ABRAMOWITZ AND STEGUN, HANDBOOK OF C MATHEMATICAL FUNCTIONS (1966) IS USED TO REDUCE THE NEGATIVE C BINOMIAL DISTRIBUTION TO THE CUMULATIVE BETA DISTRIBUTION. C C********************************************************************** C .. SCALAR ARGUMENTS .. DOUBLE PRECISION PR,S,XN C .. C .. EXTERNAL FUNCTIONS .. DOUBLE PRECISION DUMBET EXTERNAL DUMBET C .. C .. EXECUTABLE STATEMENTS .. CUMNBN = DUMBET(PR,XN,S+1.D0) RETURN END C-------------------------------------------------------------- DOUBLE PRECISION FUNCTION CUMNOR(X) C C********************************************************************** C C DOUBLE PRECISION FUNCTION CUMNOR(X) C C C FUNCTION C C C COMPUTES THE CUMULATIVE OF THE NORMAL DISTRIBUTION, I.E., C THE INTEGRAL FROM -INFINITY TO X OF C (1/SQRT(2*PI)) EXP(-U*U/2) DU C C C METHOD C C C CALLS THE DOUBLE PRECISION VERSION OF THIS ROUTINE DUMNOR. C C C ARGUMENTS C C C X --> VALUE AT WHICH CUMULATIVE NORMAL IS EVALUATED C DOUBLE PRECISION X C C********************************************************************** C .. C .. SCALAR ARGUMENTS .. DOUBLE PRECISION X C .. C .. EXTERNAL FUNCTIONS .. DOUBLE PRECISION DUMNOR EXTERNAL DUMNOR C .. C .. EXECUTABLE STATEMENTS .. CUMNOR = DUMNOR(X) RETURN END C------------------------------------------------------------------- DOUBLE PRECISION FUNCTION CUMPOI(S,XLAM) C C********************************************************************** C C DOUBLE PRECISION FUNCTION CUMPOI(S,XLAM) C CUMULATIVE POISSON DISTRIBUTION C C C FUNCTION C C C RETURNS THE PROBABILITY OF S OR FEWER EVENTS IN A POISSON C DISTRIBUTION WITH MEAN XLAM. C C C ARGUMENTS C C C S --> UPPER LIMIT OF CUMULATION OF THE POISSON. C S IS DOUBLE PRECISION C C XLAM --> MEAN OF THE POISSON DISTRIBUTION. C XLAM IS DOUBLE PRECISION C C C METHOD C C C USES FORMULA 26.4.21 OF ABRAMOWITZ AND STEGUN, HANDBOOK OF C MATHEMATICAL FUNCTIONS TO REDUCE THE CUMULATIVE POISSON TO C THE CUMULATIVE CHI-SQUARE DISTRIBUTION. C C********************************************************************** C .. SCALAR ARGUMENTS .. DOUBLE PRECISION S,XLAM C .. C .. LOCAL SCALARS .. DOUBLE PRECISION CHI,DF C .. C .. EXTERNAL FUNCTIONS .. DOUBLE PRECISION CUMCHI EXTERNAL CUMCHI C .. C .. EXECUTABLE STATEMENTS .. DF = 2.0D0* (S+1.0D0) CHI = 2.0D0*XLAM CUMPOI = 1.0D0 - CUMCHI(CHI,DF) RETURN END C----------------------------------------------------------- DOUBLE PRECISION FUNCTION CUMT(T,DF) C C********************************************************************** C C DOUBLE PRECISION FUNCTION CUMT(T,DF) C CUMULATIVE T-DISTRIBUTION C C C FUNCTION C C C COMPUTES THE INTEGRAL FROM -INFINITY TO T OF THE T-DENSITY. C C C ARGUMENTS C C C T --> UPPER LIMIT OF INTEGRATION OF THE T-DENSITY. C T IS DOUBLE PRECISION C C DF --> DEGREES OF FREEDOM OF THE T-DISTRIBUTION. C DF IS DOUBLE PRECISION C C C METHOD C C C FORMULA 26.5.27 OF ABRAMOWITZ AND STEGUN, HANDBOOK OF C MATHEMATICAL FUNCTIONS IS USED TO REDUCE THE T-DISTRIBUTION C TO AN INCOMPLETE BETA. C C********************************************************************** C .. SCALAR ARGUMENTS .. DOUBLE PRECISION DF,T C .. C .. LOCAL SCALARS .. DOUBLE PRECISION XX DOUBLE PRECISION A C .. C .. EXTERNAL FUNCTIONS .. DOUBLE PRECISION DUMBET EXTERNAL DUMBET C .. C .. INTRINSIC FUNCTIONS .. INTRINSIC DBLE,REAL C .. C .. EXECUTABLE STATEMENTS .. XX = DF/ (DF+T**2D0) A = 1.D0-DUMBET(XX,DF/2.0D0,0.5D0) IF (.NOT. (T.LE.0.0D0)) GO TO 10 CUMT = 0.5D0* (1.0D0-A) GO TO 20 10 CUMT = 0.5D0* (1.0D0+A) 20 RETURN END C------------------------------------------------------- DOUBLE PRECISION FUNCTION CXNCHN(PNONC) C C********************************************************************** C C RENAMES CUMCHN WITHOUT THE EXTRA PARAMETERS SO THAT QSTINV C CAN BE USED ON IT. C C********************************************************************** C C .. SCALAR ARGUMENTS .. DOUBLE PRECISION PNONC,XCS,XDF C .. C .. LOCAL SCALARS .. DOUBLE PRECISION CS,DF C .. C .. EXTERNAL FUNCTIONS .. DOUBLE PRECISION CUMCHN EXTERNAL CUMCHN C .. C .. ENTRY POINTS .. DOUBLE PRECISION CTNCHN C .. C .. SAVE STATEMENT .. SAVE CS,DF C .. C .. EXECUTABLE STATEMENTS .. CXNCHN = CUMCHN(CS,DF,PNONC) RETURN C ENTRY CTNCHN(XCS,XDF) CS = XCS DF = XDF CTNCHN = 0.0D0 RETURN END C---------------------------------------- DOUBLE PRECISION FUNCTION DEVLPL(A,N,X) C C********************************************************************** C C DOUBLE PRECISION FUNCTION DEVLPL(A,N,X) C DOUBLE PRECISION EVALUATE A POLYNOMIAL AT X C C C FUNCTION C C C RETURNS C A(1) + A(2)*X + ... + A(N)*X**(N-1) C C C ARGUMENTS C C C A --> ARRAY OF COEFFICIENTS OF THE POLYNOMIAL. C A IS DOUBLE PRECISION(N) C C N --> LENGTH OF A, ALSO DEGREE OF POLYNOMIAL - 1. C N IS INTEGER C C X --> POINT AT WHICH THE POLYNOMIAL IS TO BE EVALUATED. C X IS DOUBLE PRECISION C C********************************************************************** C C .. SCALAR ARGUMENTS .. DOUBLE PRECISION X INTEGER N C .. C .. ARRAY ARGUMENTS .. DOUBLE PRECISION A(N) C .. C .. LOCAL SCALARS .. DOUBLE PRECISION TERM INTEGER I C .. C .. EXECUTABLE STATEMENTS .. TERM = A(N) DO 10,I = N - 1,1,-1 TERM = A(I) + TERM*X 10 CONTINUE DEVLPL = TERM RETURN END C--------------------------------------------------------------------- DOUBLE PRECISION FUNCTION DLN1MX(X) C********************************************************************** C C DOUBLE PRECISION FUNCTION DLN1MX(X) C DOUBLE PRECISION LN(1-X) C C C FUNCTION C C C RETURNS LN(1-X) FOR SMALL X (GOOD ACCURACY IF X .LE. 0.1). C NOTE THAT THE OBVIOUS CODE OF C LOG(1.0-X) C WON'T WORK FOR SMALL X BECAUSE 1.0-X LOSES ACCURACY C C C ARGUMENTS C C C X --> VALUE FOR WHICH LN(1-X) IS DESIRED. C X IS DOUBLE PRECISION C C C METHOD C C C IF X > 0.1, THE OBVIOUS CODE ABOVE IS USED ELSE C THE TAYLOR SERIES FOR 1-X IS EXPANDED TO 20 TERMS. C C********************************************************************** C .. SCALAR ARGUMENTS .. DOUBLE PRECISION X C .. C .. EXTERNAL FUNCTIONS .. DOUBLE PRECISION DLN1PX EXTERNAL DLN1PX C .. C .. EXECUTABLE STATEMENTS .. DLN1MX = DLN1PX(-X) RETURN END C----------------------------------------------- DOUBLE PRECISION FUNCTION DLN1PX(A) C********************************************************************** C C DOUBLE PRECISION FUNCTION DLN1PX(X) C DOUBLE PRECISION LN(1+X) C C C FUNCTION C C C RETURNS LN(1+X) C NOTE THAT THE OBVIOUS CODE OF C LOG(1.0+X) C WON'T WORK FOR SMALL X BECAUSE 1.0+X LOSES ACCURACY C C C ARGUMENTS C C C X --> VALUE FOR WHICH LN(1-X) IS DESIRED. C X IS DOUBLE PRECISION C C C METHOD C C C RENAMES ALNREL FROM: C DIDINATO, A. R. AND MORRIS, A. H. ALGORITHM 708: SIGNIFICANT C DIGIT COMPUTATION OF THE INCOMPLETE BETA FUNCTION RATIOS. ACM C TRANS. MATH. SOFTW. 18 (1993), 360-373. C C********************************************************************** C----------------------------------------------------------------------- C EVALUATION OF THE FUNCTION LN(1 + A) C----------------------------------------------------------------------- C .. SCALAR ARGUMENTS .. DOUBLE PRECISION A C .. C .. LOCAL SCALARS .. DOUBLE PRECISION P1,P2,P3,Q1,Q2,Q3,T,T2,W,X C .. C .. INTRINSIC FUNCTIONS .. INTRINSIC ABS,DBLE,DLOG C .. C .. DATA STATEMENTS .. DATA P1/-.129418923021993D+01/,P2/.405303492862024D+00/, + P3/-.178874546012214D-01/ DATA Q1/-.162752256355323D+01/,Q2/.747811014037616D+00/, + Q3/-.845104217945565D-01/ C .. C .. EXECUTABLE STATEMENTS .. C-------------------------- IF (ABS(A).GT.0.375D0) GO TO 10 T = A/ (A+2.0D0) T2 = T*T W = (((P3*T2+P2)*T2+P1)*T2+1.0D0)/ (((Q3*T2+Q2)*T2+Q1)*T2+1.0D0) DLN1PX = 2.0D0*T*W RETURN C 10 X = 1.D0 + DBLE(A) DLN1PX = DLOG(X) RETURN END C--------------------------------------------------------------------------- DOUBLE PRECISION FUNCTION DLANOR(X) C********************************************************************** C C DOUBLE PRECISION FUNCTION DLANOR( X ) C DOUBLE PRECISION LOGARITH OF THE ASYMPTOTIC NORMAL C C C FUNCTION C C C COMPUTES THE LOGARITHM OF THE CUMULATIVE NORMAL DISTRIBUTION C FROM ABS( X ) TO INFINITY FOR ABS( X ) >= 5. C C C ARGUMENTS C C C X --> VALUE AT WHICH CUMULATIVE NORMAL TO BE EVALUATED C DOUBLE PRECISION X C C C METHOD C C C 23 TERM EXPANSION OF FORMULA 26.2.12 OF ABRAMOWITZ AND STEGUN. C THE RELATIVE ERROR AT X = 5 IS ABOUT 0.5E-5. C C C NOTE C C C ABS(X) MUST BE >= 5 ELSE THERE IS AN ERROR STOP. C C********************************************************************** C .. PARAMETERS .. DOUBLE PRECISION DLSQPI PARAMETER (DLSQPI=0.91893853320467274177D0) C .. C .. SCALAR ARGUMENTS .. DOUBLE PRECISION X C .. C .. LOCAL SCALARS .. DOUBLE PRECISION APPROX,CORREC,XX,XX2 C .. C .. LOCAL ARRAYS .. DOUBLE PRECISION COEF(12) C .. C .. EXTERNAL FUNCTIONS .. DOUBLE PRECISION DEVLPL,DLN1PX EXTERNAL DEVLPL,DLN1PX C .. C .. INTRINSIC FUNCTIONS .. INTRINSIC ABS,LOG C .. C .. DATA STATEMENTS .. DATA COEF/-1.0D0,3.0D0,-15.0D0,105.0D0,-945.0D0,10395.0D0, + -135135.0D0,2027025.0D0,-34459425.0D0,654729075.0D0, + -13749310575D0,316234143225.0D0/ C .. C .. EXECUTABLE STATEMENTS .. XX = ABS(X) IF (XX.LT.5.0D0) STOP ' ARGUMENT TOO SMALL IN DLANOR' APPROX = -DLSQPI - 0.5*XX*XX - LOG(XX) XX2 = XX*XX CORREC = DEVLPL(COEF,12,1.0D0/XX2)/XX2 CORREC = DLN1PX(CORREC) DLANOR = APPROX + CORREC RETURN END C--------------------------------------------------------------------- DOUBLE PRECISION FUNCTION DLNBET(A0,B0) C********************************************************************** C C DOUBLE PRECISION FUNCTION DLNBET( A, B ) C DOUBLE PRECISION LN OF THE COMPLETE BETA C C C FUNCTION C C C RETURNS THE NATURAL LOG OF THE COMPLETE BETA FUNCTION, C I.E., C C LN( GAMMA(A)*GAMMA(B) / GAMMA(A+B) C C C ARGUMENTS C C C A,B --> THE (SYMMETRIC) ARGUMENTS TO THE COMPLETE BETA C DOUBLE PRECISION A, B C C C METHOD C C C RENAMES BETALN FROM: C DIDINATO, A. R. AND MORRIS, A. H. ALGORITHM 708: SIGNIFICANT C DIGIT COMPUTATION OF THE INCOMPLETE BETA FUNCTION RATIOS. ACM C TRANS. MATH. SOFTW. 18 (1993), 360-373. C C********************************************************************** C----------------------------------------------------------------------- C EVALUATION OF THE LOGARITHM OF THE BETA FUNCTION C----------------------------------------------------------------------- C E = 0.5*LN(2*PI) C-------------------------- C .. SCALAR ARGUMENTS .. DOUBLE PRECISION A0,B0 C .. C .. LOCAL SCALARS .. DOUBLE PRECISION A,B,C,E,H,U,V,W,Z INTEGER I,N C .. C .. EXTERNAL FUNCTIONS .. DOUBLE PRECISION ALGDIV,ALNREL,BCORR,GAMLN,GSUMLN EXTERNAL ALGDIV,ALNREL,BCORR,GAMLN,GSUMLN C .. C .. INTRINSIC FUNCTIONS .. INTRINSIC DLOG,DMAX1,DMIN1 C .. C .. DATA STATEMENTS .. DATA E/.918938533204673D0/ C .. C .. EXECUTABLE STATEMENTS .. C-------------------------- A = DMIN1(A0,B0) B = DMAX1(A0,B0) IF (A.GE.8.0D0) GO TO 100 IF (A.GE.1.0D0) GO TO 20 C----------------------------------------------------------------------- C PROCEDURE WHEN A .LT. 1 C----------------------------------------------------------------------- IF (B.GE.8.0D0) GO TO 10 DLNBET = GAMLN(A) + (GAMLN(B)-GAMLN(A+B)) RETURN 10 DLNBET = GAMLN(A) + ALGDIV(A,B) RETURN C----------------------------------------------------------------------- C PROCEDURE WHEN 1 .LE. A .LT. 8 C----------------------------------------------------------------------- 20 IF (A.GT.2.0D0) GO TO 40 IF (B.GT.2.0D0) GO TO 30 DLNBET = GAMLN(A) + GAMLN(B) - GSUMLN(A,B) RETURN 30 W = 0.0D0 IF (B.LT.8.0D0) GO TO 60 DLNBET = GAMLN(A) + ALGDIV(A,B) RETURN C C REDUCTION OF A WHEN B .LE. 1000 C 40 IF (B.GT.1000.0D0) GO TO 80 N = A - 1.0D0 W = 1.0D0 DO 50 I = 1,N A = A - 1.0D0 H = A/B W = W* (H/ (1.0D0+H)) 50 CONTINUE W = DLOG(W) IF (B.LT.8.0D0) GO TO 60 DLNBET = W + GAMLN(A) + ALGDIV(A,B) RETURN C C REDUCTION OF B WHEN B .LT. 8 C 60 N = B - 1.0D0 Z = 1.0D0 DO 70 I = 1,N B = B - 1.0D0 Z = Z* (B/ (A+B)) 70 CONTINUE DLNBET = W + DLOG(Z) + (GAMLN(A)+ (GAMLN(B)-GSUMLN(A,B))) RETURN C C REDUCTION OF A WHEN B .GT. 1000 C 80 N = A - 1.0D0 W = 1.0D0 DO 90 I = 1,N A = A - 1.0D0 W = W* (A/ (1.0D0+A/B)) 90 CONTINUE DLNBET = (DLOG(W)-N*DLOG(B)) + (GAMLN(A)+ALGDIV(A,B)) RETURN C----------------------------------------------------------------------- C PROCEDURE WHEN A .GE. 8 C----------------------------------------------------------------------- 100 W = BCORR(A,B) H = A/B C = H/ (1.0D0+H) U = - (A-0.5D0)*DLOG(C) V = B*ALNREL(H) IF (U.LE.V) GO TO 110 DLNBET = (((-0.5D0*DLOG(B)+E)+W)-V) - U RETURN 110 DLNBET = (((-0.5D0*DLOG(B)+E)+W)-U) - V RETURN END C------------------------------------------------------------------------------ DOUBLE PRECISION FUNCTION DLNGAM(A) C********************************************************************** C C DOUBLE PRECISION FUNCTION DLNGAM(X) C DOUBLE PRECISION LN OF THE GAMMA FUNCTION C C C FUNCTION C C C RETURNS THE NATURAL LOGARITHM OF GAMMA(X). C C C ARGUMENTS C C C X --> VALUE AT WHICH SCALED LOG GAMMA IS TO BE RETURNED C X IS DOUBLE PRECISION C C C METHOD C C C RENAMES GAMLN FROM: C DIDINATO, A. R. AND MORRIS, A. H. ALGORITHM 708: SIGNIFICANT C DIGIT COMPUTATION OF THE INCOMPLETE BETA FUNCTION RATIOS. ACM C TRANS. MATH. SOFTW. 18 (1993), 360-373. C C********************************************************************** C----------------------------------------------------------------------- C EVALUATION OF LN(GAMMA(A)) FOR POSITIVE A C----------------------------------------------------------------------- C WRITTEN BY ALFRED H. MORRIS C NAVAL SURFACE WARFARE CENTER C DAHLGREN, VIRGINIA C-------------------------- C D = 0.5*(LN(2*PI) - 1) C-------------------------- C .. SCALAR ARGUMENTS .. DOUBLE PRECISION A C .. C .. LOCAL SCALARS .. DOUBLE PRECISION C0,C1,C2,C3,C4,C5,D,T,W INTEGER I,N C .. C .. EXTERNAL FUNCTIONS .. DOUBLE PRECISION GAMLN1 EXTERNAL GAMLN1 C .. C .. INTRINSIC FUNCTIONS .. INTRINSIC DLOG C .. C .. DATA STATEMENTS .. C-------------------------- DATA D/.418938533204673D0/ DATA C0/.833333333333333D-01/,C1/-.277777777760991D-02/, + C2/.793650666825390D-03/,C3/-.595202931351870D-03/, + C4/.837308034031215D-03/,C5/-.165322962780713D-02/ C .. C .. EXECUTABLE STATEMENTS .. C----------------------------------------------------------------------- IF (A.GT.0.8D0) GO TO 10 DLNGAM = GAMLN1(A) - DLOG(A) RETURN 10 IF (A.GT.2.25D0) GO TO 20 T = (A-0.5D0) - 0.5D0 DLNGAM = GAMLN1(T) RETURN C 20 IF (A.GE.10.0D0) GO TO 40 N = A - 1.25D0 T = A W = 1.0D0 DO 30 I = 1,N T = T - 1.0D0 W = T*W 30 CONTINUE DLNGAM = GAMLN1(T-1.0D0) + DLOG(W) RETURN C 40 T = (1.0D0/A)**2 W = (((((C5*T+C4)*T+C3)*T+C2)*T+C1)*T+C0)/A DLNGAM = (D+W) + (A-0.5D0)* (DLOG(A)-1.0D0) END C---------------------------------------------------------- DOUBLE PRECISION FUNCTION DSTREM(Z) IMPLICIT DOUBLE PRECISION (A-H,O-P,R-Z),INTEGER (I-N),LOGICAL (Q) C********************************************************************** C C DOUBLE PRECISION FUNCTION DSTREM( Z ) C DOUBLE PRECISION STERLING REMAINDER C C C FUNCTION C C C RETURNS LOG(GAMMA(Z)) - STERLING(Z) C WHERE STERLING(Z) IS STERLING'S APPROXIMATION C TO LOG(GAMMA(Z)) C C C ARGUMENTS C C C Z --> VALUE AT WHICH STERLING REMAINDER CALCULATED C MUST BE POSITIVE. C DOUBLE PRECISION Z C C C METHOD C C C C IF Z >= 6 USES 9 TERMS OF SERIES IN BERNOULLI NUMBERS C (VALUES CALCULATED USING MAPLE) C OTHERWISE COMPUTES DIFFERENCE EXPLICITLY C C********************************************************************** C .. PARAMETERS .. DOUBLE PRECISION HLN2PI PARAMETER (HLN2PI=0.91893853320467274178D0) INTEGER NCOEF PARAMETER (NCOEF=10) C .. C .. SCALAR ARGUMENTS .. DOUBLE PRECISION Z C .. C .. LOCAL SCALARS .. DOUBLE PRECISION STERL C .. C .. LOCAL ARRAYS .. DOUBLE PRECISION COEF(NCOEF) C .. C .. EXTERNAL FUNCTIONS .. DOUBLE PRECISION DEVLPL,DLNGAM EXTERNAL DEVLPL,DLNGAM C .. C .. INTRINSIC FUNCTIONS .. INTRINSIC LOG C .. C .. DATA STATEMENTS .. DATA COEF/0.0D0,0.0833333333333333333333333333333D0, + -0.00277777777777777777777777777778D0, + 0.000793650793650793650793650793651D0, + -0.000595238095238095238095238095238D0, + 0.000841750841750841750841750841751D0, + -0.00191752691752691752691752691753D0, + 0.00641025641025641025641025641026D0, + -0.0295506535947712418300653594771D0, + 0.179644372368830573164938490016D0/ C .. C .. EXECUTABLE STATEMENTS .. C FOR INFORMATION, HERE ARE THE NEXT 11 COEFFICIENTS OF THE C REMAINDER TERM IN STERLING'S FORMULA C -1.39243221690590111642743221691 C 13.4028640441683919944789510007 C -156.848284626002017306365132452 C 2193.10333333333333333333333333 C -36108.7712537249893571732652192 C 691472.268851313067108395250776 C -0.152382215394074161922833649589D8 C 0.382900751391414141414141414141D9 C -0.108822660357843910890151491655D11 C 0.347320283765002252252252252252D12 C -0.123696021422692744542517103493D14 C IF (Z.LE.0.0D0) STOP 'ZERO OR NEGATIVE ARGUMENT IN DSTREM' IF (.NOT. (Z.GT.6.0D0)) GO TO 10 DSTREM = DEVLPL(COEF,10,1.0D0/Z**2)*Z GO TO 20 10 STERL = HLN2PI + (Z-0.5D0)*LOG(Z) - Z DSTREM = DLNGAM(Z) - STERL 20 RETURN END C------------------------------------------------------------------------ DOUBLE PRECISION FUNCTION DUMBET(XX,XA,XB) C********************************************************************** C C DOUBLE PRECISION FUNCTION DUMBET(X,A,B) C DOUBLE PRECISION CUMULATIVE INCOMPLETE BETA DISTRIBUTION C C C FUNCTION C C C RETURNS THE CDF TO X OF THE INCOMPLETE BETA DISTRIBUTION C WITH PARAMETERS A AND B. THIS IS THE INTEGRAL FROM 0 TO X C OF (1/B(A,B))*F(T)) WHERE F(T) = T**(A-1) * (1-T)**(B-1) C C C ARGUMENTS C C C X --> UPPER LIMIT OF INTEGRATION. C X IS DOUBLE PRECISION C C A --> FIRST PARAMETER OF THE BETA DISTRIBUTION. C A IS DOUBLE PRECISION C C B --> SECOND PARAMETER OF THE BETA DISTRIBUTION. C B IS DOUBLE PRECISION C C C METHOD C C C CALLS THE ROUTINE BRATIO. C C REFERENCES C C DIDONATO, ARMIDO R. AND MORRIS, ALFRED H. JR. (1992) ALGORITHIM C 708 SIGNIFICANT DIGIT COMPUTATION OF THE INCOMPLETE BETA FUNCTION C RATIOS. ACM TOMS, VOL.18, NO. 3, SEPT. 1992, 360-373. C C********************************************************************** C .. SCALAR ARGUMENTS .. DOUBLE PRECISION XX,XA,XB C .. C .. LOCAL SCALARS .. DOUBLE PRECISION YY,W,W1 INTEGER IERR C .. C .. EXTERNAL ROUTINES .. EXTERNAL BRATIO C .. C .. EXECUTABLE STATEMENTS .. IF (.NOT. (XX.LE.0.0D0)) GO TO 10 DUMBET = 0.0D0 RETURN 10 IF (.NOT. (XX.GE.1.0D0)) GO TO 20 DUMBET = 1.0D0 RETURN 20 YY = 1.0D0 - XX C CALL BRATIO ROUTINE CALL BRATIO(XA,XB,XX,YY,W,W1,IERR) DUMBET = W RETURN END C------------------------------------------------------ DOUBLE PRECISION FUNCTION DUMGAM(XX,AA) C C********************************************************************** C C DOUBLE PRECISION FUNCTION DUMGAM(X,A) C DOUBLE PRECISION CUMULATIVE INCOMPLETE GAMMA DISTRIBUTION C C C FUNCTION C C C COMPUTES THE CUMULATIVE OF THE INCOMPLETE GAMMA C DISTRIBUTION, I.E., THE INTEGRAL FROM 0 TO X OF C (1/GAM(A))*EXP(-T)*T**(A-1) DT C WHERE GAM(A) IS THE COMPLETE GAMMA FUNCTION OF A, I.E., C GAM(A) = INTEGRAL FROM 0 TO INFINITY OF C EXP(-T)*T**(A-1) DT C C C ARGUMENTS C C C X --> THE UPPER LIMIT OF INTEGRATION OF THE INCOMPLETE GAMMA. C X IS DOUBLE PRECISION C C A --> THE SHAPE PARAMETER OF THE INCOMPLETE GAMMA. C A IS DOUBLE PRECISION C C C METHOD C C C CALLS THE ROUTINE GRATIO. C C********************************************************************** C C .. C .. SCALAR ARGUMENTS .. DOUBLE PRECISION AA,XX C .. C .. LOCAL SCALARS .. DOUBLE PRECISION ANS,QANS C .. C .. EXTERNAL ROUTINES .. EXTERNAL GRATIO C .. C .. EXECUTABLE STATEMENTS .. IF (.NOT. (XX.LE.0.D0)) GO TO 10 DUMGAM = 0.D0 RETURN C CALL GRATIO ROUTINE 10 CALL GRATIO(AA,XX,ANS,QANS,0) DUMGAM = ANS RETURN END C----------------------------------------------------- DOUBLE PRECISION FUNCTION DUMNOR(X) C********************************************************************** C C DOUBLE PRECISION FUNCTION DUMNOR(X) C C C FUNCTION C C C COMPUTES THE CUMULATIVE OF THE NORMAL DISTRIBUTION, I.E., C THE INTEGRAL FROM -INFINITY TO X OF C (1/SQRT(2*PI)) EXP(-U*U/2) DU C C C METHOD C C C THE RATIONAL FUNCTION APPROXIMATION FROM PAGES 90 - 92 OF C KENNEDY AND GENTLE, STATISTICAL COMPUTING, MARCEL DEKKER, NY C 1980. C C C ARGUMENTS C C C X --> ARGUMENT AT WHICH CUMULATIVE NORMAL IS EVALUATED C DOUBLE PRECISION X C C********************************************************************** C C C PIM12 IS PI**(-1/2) C SQRT2 IS SQRT(2) C C C C .. SCALAR ARGUMENTS .. DOUBLE PRECISION X C .. C .. LOCAL SCALARS .. DOUBLE PRECISION DERF,DERFC,PIM12,SQRT2,Z,Z2,ZM2 LOGICAL QDIRCT C .. C .. LOCAL ARRAYS .. DOUBLE PRECISION XDEN1(4),XDEN2(8),XDEN3(5),XNUM1(4),XNUM2(8), + XNUM3(5) C .. C .. EXTERNAL FUNCTIONS .. DOUBLE PRECISION DEVLPL,DLANOR EXTERNAL DEVLPL,DLANOR C .. C .. INTRINSIC FUNCTIONS .. INTRINSIC ABS,EXP C .. C .. DATA STATEMENTS .. DATA XNUM1/2.4266795523053175D2,2.1979261618294152D1, + 6.9963834886191355D0,-3.5609843701815385D-2/ DATA XDEN1/2.1505887586986120D2,9.1164905404514901D1, + 1.5082797630407787D1,1.0000000000000000D0/ DATA XNUM2/3.004592610201616005D2,4.519189537118729422D2, + 3.393208167343436870D2,1.529892850469404039D2, + 4.316222722205673530D1,7.211758250883093659D0, + 5.641955174789739711D-1,-1.368648573827167067D-7/ DATA XDEN2/3.004592609569832933D2,7.909509253278980272D2, + 9.313540948506096211D2,6.389802644656311665D2, + 2.775854447439876434D2,7.700015293522947295D1, + 1.278272731962942351D1,1.000000000000000000D0/ DATA XNUM3/-2.99610707703542174D-3,-4.94730910623250734D-2, + -2.26956593539686930D-1,-2.78661308609647788D-1, + -2.23192459734184686D-2/ DATA XDEN3/1.06209230528467918D-2,1.91308926107829841D-1, + 1.05167510706793207D0,1.98733201817135256D0, + 1.00000000000000000D0/ DATA PIM12/0.5641895835477562869480795D0/ DATA SQRT2/1.4142135623730950488D0/ C .. C .. EXECUTABLE STATEMENTS .. IF (.NOT. (ABS(X).LT.1.0E-30)) GO TO 10 DUMNOR = 0.5 RETURN GO TO 50 10 IF (.NOT. (X.LE.-38.0)) GO TO 20 DUMNOR = 0.0 RETURN GO TO 50 20 IF (.NOT. (X.LE.-15.0)) GO TO 30 DUMNOR = EXP(DLANOR(X)) RETURN GO TO 50 30 IF (.NOT. (X.GT.6.0)) GO TO 40 DUMNOR = 1.0 RETURN GO TO 50 40 CONTINUE 50 Z = ABS(X/SQRT2) Z2 = Z*Z ZM2 = 1.0D0/Z2 IF (Z.LT.0.5D0) THEN DERF = Z*DEVLPL(XNUM1,4,Z2)/DEVLPL(XDEN1,4,Z2) QDIRCT = .TRUE. ELSE IF (Z.LT.4.0D0) THEN DERFC = EXP(-Z2)*DEVLPL(XNUM2,8,Z)/DEVLPL(XDEN2,8,Z) QDIRCT = .FALSE. ELSE DERFC = (EXP(-Z2)/Z)* (PIM12+ZM2*DEVLPL(XNUM3,5,ZM2)/ + DEVLPL(XDEN3,5,ZM2)) QDIRCT = .FALSE. END IF IF (.NOT. (X.GE.0.0)) GO TO 60 IF (.NOT. (QDIRCT)) DERF = 1.0D0 - DERFC DUMNOR = (1.0D0+DERF)/2.0D0 GO TO 70 60 IF (QDIRCT) DERFC = 1.0D0 - DERF DUMNOR = DERFC/2.0D0 70 RETURN END C-------------------------------------------------------------------- DOUBLE PRECISION FUNCTION ERF(X) C----------------------------------------------------------------------- C EVALUATION OF THE REAL ERROR FUNCTION C----------------------------------------------------------------------- C .. SCALAR ARGUMENTS .. DOUBLE PRECISION X C .. C .. LOCAL SCALARS .. DOUBLE PRECISION AX,BOT,C,T,TOP,X2 C .. C .. LOCAL ARRAYS .. DOUBLE PRECISION A(5),B(3),P(8),Q(8),R(5),S(4) C .. C .. INTRINSIC FUNCTIONS .. INTRINSIC ABS,EXP,SIGN C .. C .. DATA STATEMENTS .. C------------------------- C------------------------- C------------------------- C------------------------- DATA C/.564189583547756D0/ DATA A(1)/.771058495001320D-04/,A(2)/-.133733772997339D-02/, + A(3)/.323076579225834D-01/,A(4)/.479137145607681D-01/, + A(5)/.128379167095513D+00/ DATA B(1)/.301048631703895D-02/,B(2)/.538971687740286D-01/, + B(3)/.375795757275549D+00/ DATA P(1)/-1.36864857382717D-07/,P(2)/5.64195517478974D-01/, + P(3)/7.21175825088309D+00/,P(4)/4.31622272220567D+01/, + P(5)/1.52989285046940D+02/,P(6)/3.39320816734344D+02/, + P(7)/4.51918953711873D+02/,P(8)/3.00459261020162D+02/ DATA Q(1)/1.00000000000000D+00/,Q(2)/1.27827273196294D+01/, + Q(3)/7.70001529352295D+01/,Q(4)/2.77585444743988D+02/, + Q(5)/6.38980264465631D+02/,Q(6)/9.31354094850610D+02/, + Q(7)/7.90950925327898D+02/,Q(8)/3.00459260956983D+02/ DATA R(1)/2.10144126479064D+00/,R(2)/2.62370141675169D+01/, + R(3)/2.13688200555087D+01/,R(4)/4.65807828718470D+00/, + R(5)/2.82094791773523D-01/ DATA S(1)/9.41537750555460D+01/,S(2)/1.87114811799590D+02/, + S(3)/9.90191814623914D+01/,S(4)/1.80124575948747D+01/ C .. C .. EXECUTABLE STATEMENTS .. C------------------------- AX = ABS(X) IF (AX.GT.0.5D0) GO TO 10 T = X*X TOP = ((((A(1)*T+A(2))*T+A(3))*T+A(4))*T+A(5)) + 1.0D0 BOT = ((B(1)*T+B(2))*T+B(3))*T + 1.0D0 ERF = X* (TOP/BOT) RETURN C 10 IF (AX.GT.4.0D0) GO TO 20 TOP = ((((((P(1)*AX+P(2))*AX+P(3))*AX+P(4))*AX+P(5))*AX+P(6))*AX+ + P(7))*AX + P(8) BOT = ((((((Q(1)*AX+Q(2))*AX+Q(3))*AX+Q(4))*AX+Q(5))*AX+Q(6))*AX+ + Q(7))*AX + Q(8) ERF = 0.5D0 + (0.5D0-EXP(-X*X)*TOP/BOT) IF (X.LT.0.0D0) ERF = -ERF RETURN C 20 IF (AX.GE.5.8D0) GO TO 30 X2 = X*X T = 1.0D0/X2 TOP = (((R(1)*T+R(2))*T+R(3))*T+R(4))*T + R(5) BOT = (((S(1)*T+S(2))*T+S(3))*T+S(4))*T + 1.0D0 ERF = (C-TOP/ (X2*BOT))/AX ERF = 0.5D0 + (0.5D0-EXP(-X2)*ERF) IF (X.LT.0.0D0) ERF = -ERF RETURN C 30 ERF = SIGN(1.0D0,X) RETURN END C-------------------------------------------------------------- DOUBLE PRECISION FUNCTION ERFC1(IND,X) C----------------------------------------------------------------------- C EVALUATION OF THE COMPLEMENTARY ERROR FUNCTION C C ERFC1(IND,X) = ERFC(X) IF IND = 0 C ERFC1(IND,X) = EXP(X*X)*ERFC(X) OTHERWISE C----------------------------------------------------------------------- C .. SCALAR ARGUMENTS .. DOUBLE PRECISION X INTEGER IND C .. C .. LOCAL SCALARS .. DOUBLE PRECISION AX,BOT,C,E,T,TOP,W C .. C .. LOCAL ARRAYS .. DOUBLE PRECISION A(5),B(3),P(8),Q(8),R(5),S(4) C .. C .. EXTERNAL FUNCTIONS .. DOUBLE PRECISION EXPARG EXTERNAL EXPARG C .. C .. INTRINSIC FUNCTIONS .. INTRINSIC ABS,DBLE,EXP C .. C .. DATA STATEMENTS .. C------------------------- C------------------------- C------------------------- C------------------------- DATA C/.564189583547756D0/ DATA A(1)/.771058495001320D-04/,A(2)/-.133733772997339D-02/, + A(3)/.323076579225834D-01/,A(4)/.479137145607681D-01/, + A(5)/.128379167095513D+00/ DATA B(1)/.301048631703895D-02/,B(2)/.538971687740286D-01/, + B(3)/.375795757275549D+00/ DATA P(1)/-1.36864857382717D-07/,P(2)/5.64195517478974D-01/, + P(3)/7.21175825088309D+00/,P(4)/4.31622272220567D+01/, + P(5)/1.52989285046940D+02/,P(6)/3.39320816734344D+02/, + P(7)/4.51918953711873D+02/,P(8)/3.00459261020162D+02/ DATA Q(1)/1.00000000000000D+00/,Q(2)/1.27827273196294D+01/, + Q(3)/7.70001529352295D+01/,Q(4)/2.77585444743988D+02/, + Q(5)/6.38980264465631D+02/,Q(6)/9.31354094850610D+02/, + Q(7)/7.90950925327898D+02/,Q(8)/3.00459260956983D+02/ DATA R(1)/2.10144126479064D+00/,R(2)/2.62370141675169D+01/, + R(3)/2.13688200555087D+01/,R(4)/4.65807828718470D+00/, + R(5)/2.82094791773523D-01/ DATA S(1)/9.41537750555460D+01/,S(2)/1.87114811799590D+02/, + S(3)/9.90191814623914D+01/,S(4)/1.80124575948747D+01/ C .. C .. EXECUTABLE STATEMENTS .. C------------------------- C C ABS(X) .LE. 0.5 C AX = ABS(X) IF (AX.GT.0.5D0) GO TO 10 T = X*X TOP = ((((A(1)*T+A(2))*T+A(3))*T+A(4))*T+A(5)) + 1.0D0 BOT = ((B(1)*T+B(2))*T+B(3))*T + 1.0D0 ERFC1 = 0.5D0 + (0.5D0-X* (TOP/BOT)) IF (IND.NE.0) ERFC1 = EXP(T)*ERFC1 RETURN C C 0.5 .LT. ABS(X) .LE. 4 C 10 IF (AX.GT.4.0D0) GO TO 20 TOP = ((((((P(1)*AX+P(2))*AX+P(3))*AX+P(4))*AX+P(5))*AX+P(6))*AX+ + P(7))*AX + P(8) BOT = ((((((Q(1)*AX+Q(2))*AX+Q(3))*AX+Q(4))*AX+Q(5))*AX+Q(6))*AX+ + Q(7))*AX + Q(8) ERFC1 = TOP/BOT GO TO 40 C C ABS(X) .GT. 4 C 20 IF (X.LE.-5.6D0) GO TO 60 IF (IND.NE.0) GO TO 30 IF (X.GT.100.0D0) GO TO 70 IF (X*X.GT.-EXPARG(1)) GO TO 70 C 30 T = (1.0D0/X)**2 TOP = (((R(1)*T+R(2))*T+R(3))*T+R(4))*T + R(5) BOT = (((S(1)*T+S(2))*T+S(3))*T+S(4))*T + 1.0D0 ERFC1 = (C-T*TOP/BOT)/AX C C FINAL ASSEMBLY C 40 IF (IND.EQ.0) GO TO 50 IF (X.LT.0.0D0) ERFC1 = 2.0D0*EXP(X*X) - ERFC1 RETURN 50 W = DBLE(X)*DBLE(X) T = W E = W - DBLE(T) ERFC1 = ((0.5D0+ (0.5D0-E))*EXP(-T))*ERFC1 IF (X.LT.0.0D0) ERFC1 = 2.0D0 - ERFC1 RETURN C C LIMIT VALUE FOR LARGE NEGATIVE X C 60 ERFC1 = 2.0D0 IF (IND.NE.0) ERFC1 = 2.0D0*EXP(X*X) RETURN C C LIMIT VALUE FOR LARGE POSITIVE X C WHEN IND = 0 C 70 ERFC1 = 0.0D0 RETURN END C------------------------------------------------------------------ DOUBLE PRECISION FUNCTION ESUM(MU,X) C----------------------------------------------------------------------- C EVALUATION OF EXP(MU + X) C----------------------------------------------------------------------- C .. SCALAR ARGUMENTS .. DOUBLE PRECISION X INTEGER MU C .. C .. LOCAL SCALARS .. DOUBLE PRECISION W C .. C .. INTRINSIC FUNCTIONS .. INTRINSIC EXP C .. C .. EXECUTABLE STATEMENTS .. IF (X.GT.0.0D0) GO TO 10 C IF (MU.LT.0) GO TO 20 W = MU + X IF (W.GT.0.0D0) GO TO 20 ESUM = EXP(W) RETURN C 10 IF (MU.GT.0) GO TO 20 W = MU + X IF (W.LT.0.0D0) GO TO 20 ESUM = EXP(W) RETURN C 20 W = MU ESUM = EXP(W)*EXP(X) RETURN END C------------------------------------------------ DOUBLE PRECISION FUNCTION EXPARG(L) C-------------------------------------------------------------------- C IF L = 0 THEN EXPARG(L) = THE LARGEST POSITIVE W FOR WHICH C EXP(W) CAN BE COMPUTED. C C IF L IS NONZERO THEN EXPARG(L) = THE LARGEST NEGATIVE W FOR C WHICH THE COMPUTED VALUE OF EXP(W) IS NONZERO. C C NOTE... ONLY AN APPROXIMATE VALUE FOR EXPARG(L) IS NEEDED. C-------------------------------------------------------------------- C .. SCALAR ARGUMENTS .. INTEGER L C .. C .. LOCAL SCALARS .. DOUBLE PRECISION LNB INTEGER B,M C .. C .. EXTERNAL FUNCTIONS .. INTEGER IPMPAR EXTERNAL IPMPAR C .. C .. INTRINSIC FUNCTIONS .. INTRINSIC DBLE,DLOG C .. C .. EXECUTABLE STATEMENTS .. C B = IPMPAR(4) IF (B.NE.2) GO TO 10 LNB = .69314718055995D0 GO TO 40 10 IF (B.NE.8) GO TO 20 LNB = 2.0794415416798D0 GO TO 40 20 IF (B.NE.16) GO TO 30 LNB = 2.7725887222398D0 GO TO 40 30 LNB = DLOG(DBLE(B)) C 40 IF (L.EQ.0) GO TO 50 M = IPMPAR(9) - 1 EXPARG = 0.99999D0* (M*LNB) RETURN 50 M = IPMPAR(10) EXPARG = 0.99999D0* (M*LNB) RETURN END C----------------------------------------------------- DOUBLE PRECISION FUNCTION FPSER(A,B,X,EPS) C----------------------------------------------------------------------- C C EVALUATION OF I (A,B) C X C C FOR B .LT. MIN(EPS,EPS*A) AND X .LE. 0.5. C C----------------------------------------------------------------------- C C SET FPSER = X**A C C .. SCALAR ARGUMENTS .. DOUBLE PRECISION A,B,EPS,X C .. C .. LOCAL SCALARS .. DOUBLE PRECISION AN,C,S,T,TOL C .. C .. EXTERNAL FUNCTIONS .. DOUBLE PRECISION EXPARG EXTERNAL EXPARG C .. C .. INTRINSIC FUNCTIONS .. INTRINSIC ABS,DLOG,EXP C .. C .. EXECUTABLE STATEMENTS .. FPSER = 1.0D0 IF (A.LE.1.D-3*EPS) GO TO 10 FPSER = 0.0D0 T = A*DLOG(X) IF (T.LT.EXPARG(1)) RETURN FPSER = EXP(T) C C NOTE THAT 1/B(A,B) = B C 10 FPSER = (B/A)*FPSER TOL = EPS/A AN = A + 1.0D0 T = X S = T/AN 20 AN = AN + 1.0D0 T = X*T C = T/AN S = S + C IF (ABS(C).GT.TOL) GO TO 20 C FPSER = FPSER* (1.0D0+A*S) RETURN END C------------------------------------------------- DOUBLE PRECISION FUNCTION GAM1(A) C ------------------------------------------------------------------ C COMPUTATION OF 1/GAMMA(A+1) - 1 FOR -0.5 .LE. A .LE. 1.5 C ------------------------------------------------------------------ C .. SCALAR ARGUMENTS .. DOUBLE PRECISION A C .. C .. LOCAL SCALARS .. DOUBLE PRECISION BOT,D,S1,S2,T,TOP,W C .. C .. LOCAL ARRAYS .. DOUBLE PRECISION P(7),Q(5),R(9) C .. C .. DATA STATEMENTS .. C ------------------- C ------------------- C ------------------- C ------------------- DATA P(1)/.577215664901533D+00/,P(2)/-.409078193005776D+00/, + P(3)/-.230975380857675D+00/,P(4)/.597275330452234D-01/, + P(5)/.766968181649490D-02/,P(6)/-.514889771323592D-02/, + P(7)/.589597428611429D-03/ DATA Q(1)/.100000000000000D+01/,Q(2)/.427569613095214D+00/, + Q(3)/.158451672430138D+00/,Q(4)/.261132021441447D-01/, + Q(5)/.423244297896961D-02/ DATA R(1)/-.422784335098468D+00/,R(2)/-.771330383816272D+00/, + R(3)/-.244757765222226D+00/,R(4)/.118378989872749D+00/, + R(5)/.930357293360349D-03/,R(6)/-.118290993445146D-01/, + R(7)/.223047661158249D-02/,R(8)/.266505979058923D-03/, + R(9)/-.132674909766242D-03/ DATA S1/.273076135303957D+00/,S2/.559398236957378D-01/ C .. C .. EXECUTABLE STATEMENTS .. C ------------------- T = A D = A - 0.5D0 IF (D.GT.0.0D0) T = D - 0.5D0 IF (T) 40,10,20 C 10 GAM1 = 0.0D0 RETURN C 20 TOP = (((((P(7)*T+P(6))*T+P(5))*T+P(4))*T+P(3))*T+P(2))*T + P(1) BOT = (((Q(5)*T+Q(4))*T+Q(3))*T+Q(2))*T + 1.0D0 W = TOP/BOT IF (D.GT.0.0D0) GO TO 30 GAM1 = A*W RETURN 30 GAM1 = (T/A)* ((W-0.5D0)-0.5D0) RETURN C 40 TOP = (((((((R(9)*T+R(8))*T+R(7))*T+R(6))*T+R(5))*T+R(4))*T+R(3))* + T+R(2))*T + R(1) BOT = (S2*T+S1)*T + 1.0D0 W = TOP/BOT IF (D.GT.0.0D0) GO TO 50 GAM1 = A* ((W+0.5D0)+0.5D0) RETURN 50 GAM1 = T*W/A RETURN END C----------------------------------------------------------------------------- SUBROUTINE GAMINV(A,X,X0,P,Q,IERR) C ---------------------------------------------------------------------- C INVERSE INCOMPLETE GAMMA RATIO FUNCTION C C GIVEN POSITIVE A, AND NONEGATIVE P AND Q WHERE P + Q = 1. C THEN X IS COMPUTED WHERE P(A,X) = P AND Q(A,X) = Q. SCHRODER C ITERATION IS EMPLOYED. THE ROUTINE ATTEMPTS TO COMPUTE X C TO 10 SIGNIFICANT DIGITS IF THIS IS POSSIBLE FOR THE C PARTICULAR COMPUTER ARITHMETIC BEING USED. C C ------------ C C X IS A VARIABLE. IF P = 0 THEN X IS ASSIGNED THE VALUE 0, C AND IF Q = 0 THEN X IS SET TO THE LARGEST FLOATING POINT C NUMBER AVAILABLE. OTHERWISE, GAMINV ATTEMPTS TO OBTAIN C A SOLUTION FOR P(A,X) = P AND Q(A,X) = Q. IF THE ROUTINE C IS SUCCESSFUL THEN THE SOLUTION IS STORED IN X. C C X0 IS AN OPTIONAL INITIAL APPROXIMATION FOR X. IF THE USER C DOES NOT WISH TO SUPPLY AN INITIAL APPROXIMATION, THEN SET C X0 .LE. 0. C C IERR IS A VARIABLE THAT REPORTS THE STATUS OF THE RESULTS. C WHEN THE ROUTINE TERMINATES, IERR HAS ONE OF THE FOLLOWING C VALUES ... C C IERR = 0 THE SOLUTION WAS OBTAINED. ITERATION WAS C NOT USED. C IERR.GT.0 THE SOLUTION WAS OBTAINED. IERR ITERATIONS C WERE PERFORMED. C IERR = -2 (INPUT ERROR) A .LE. 0 C IERR = -3 NO SOLUTION WAS OBTAINED. THE RATIO Q/A C IS TOO LARGE. C IERR = -4 (INPUT ERROR) P + Q .NE. 1 C IERR = -6 20 ITERATIONS WERE PERFORMED. THE MOST C RECENT VALUE OBTAINED FOR X IS GIVEN. C THIS CANNOT OCCUR IF X0 .LE. 0. C IERR = -7 ITERATION FAILED. NO VALUE IS GIVEN FOR X. C THIS MAY OCCUR WHEN X IS APPROXIMATELY 0. C IERR = -8 A VALUE FOR X HAS BEEN OBTAINED, BUT THE C ROUTINE IS NOT CERTAIN OF ITS ACCURACY. C ITERATION CANNOT BE PERFORMED IN THIS C CASE. IF X0 .LE. 0, THIS CAN OCCUR ONLY C WHEN P OR Q IS APPROXIMATELY 0. IF X0 IS C POSITIVE THEN THIS CAN OCCUR WHEN A IS C EXCEEDINGLY CLOSE TO X AND A IS EXTREMELY C LARGE (SAY A .GE. 1.E20). C ---------------------------------------------------------------------- C WRITTEN BY ALFRED H. MORRIS, JR. C NAVAL SURFACE WEAPONS CENTER C DAHLGREN, VIRGINIA C ------------------- C .. SCALAR ARGUMENTS .. DOUBLE PRECISION A,P,Q,X,X0 INTEGER IERR C .. C .. LOCAL SCALARS .. DOUBLE PRECISION A0,A1,A2,A3,AM1,AMAX,AP1,AP2,AP3,APN,B,B1,B2,B3, + B4,C,C1,C2,C3,C4,C5,D,E,E2,EPS,G,H,LN10,PN,QG,QN, + R,RTA,S,S2,SUM,T,TOL,U,W,XMAX,XMIN,XN,Y,Z INTEGER IOP C .. C .. LOCAL ARRAYS .. DOUBLE PRECISION AMIN(2),BMIN(2),DMIN(2),EMIN(2),EPS0(2) C .. C .. EXTERNAL FUNCTIONS .. DOUBLE PRECISION ALNREL,GAMLN,GAMLN1,GAMMA,RCOMP,SPMPAR EXTERNAL ALNREL,GAMLN,GAMLN1,GAMMA,RCOMP,SPMPAR C .. C .. EXTERNAL SUBROUTINES .. EXTERNAL GRATIO C .. C .. INTRINSIC FUNCTIONS .. INTRINSIC ABS,DBLE,DLOG,DMAX1,EXP,SQRT C .. C .. DATA STATEMENTS .. C ------------------- C LN10 = LN(10) C C = EULER CONSTANT C ------------------- C ------------------- C ------------------- C ------------------- DATA LN10/2.302585D0/ DATA C/.577215664901533D0/ DATA A0/3.31125922108741D0/,A1/11.6616720288968D0/, + A2/4.28342155967104D0/,A3/.213623493715853D0/ DATA B1/6.61053765625462D0/,B2/6.40691597760039D0/, + B3/1.27364489782223D0/,B4/.036117081018842D0/ DATA EPS0(1)/1.D-10/,EPS0(2)/1.D-08/ DATA AMIN(1)/500.0D0/,AMIN(2)/100.0D0/ DATA BMIN(1)/1.D-28/,BMIN(2)/1.D-13/ DATA DMIN(1)/1.D-06/,DMIN(2)/1.D-04/ DATA EMIN(1)/2.D-03/,EMIN(2)/6.D-03/ DATA TOL/1.D-5/ C .. C .. EXECUTABLE STATEMENTS .. C ------------------- C ****** E, XMIN, AND XMAX ARE MACHINE DEPENDENT CONSTANTS. C E IS THE SMALLEST NUMBER FOR WHICH 1.0 + E .GT. 1.0. C XMIN IS THE SMALLEST POSITIVE NUMBER AND XMAX IS THE C LARGEST POSITIVE NUMBER. C E = SPMPAR(1) XMIN = SPMPAR(2) XMAX = SPMPAR(3) C ------------------- X = 0.0D0 IF (A.LE.0.0D0) GO TO 300 T = DBLE(P) + DBLE(Q) - 1.D0 IF (ABS(T).GT.E) GO TO 320 C IERR = 0 IF (P.EQ.0.0D0) RETURN IF (Q.EQ.0.0D0) GO TO 270 IF (A.EQ.1.0D0) GO TO 280 C E2 = 2.0D0*E AMAX = 0.4D-10/ (E*E) IOP = 1 IF (E.GT.1.D-10) IOP = 2 EPS = EPS0(IOP) XN = X0 IF (X0.GT.0.0D0) GO TO 160 C C SELECTION OF THE INITIAL APPROXIMATION XN OF X C WHEN A .LT. 1 C IF (A.GT.1.0D0) GO TO 80 G = GAMMA(A+1.0D0) QG = Q*G IF (QG.EQ.0.0D0) GO TO 360 B = QG/A IF (QG.GT.0.6D0*A) GO TO 40 IF (A.GE.0.30D0 .OR. B.LT.0.35D0) GO TO 10 T = EXP(- (B+C)) U = T*EXP(T) XN = T*EXP(U) GO TO 160 C 10 IF (B.GE.0.45D0) GO TO 40 IF (B.EQ.0.0D0) GO TO 360 Y = -DLOG(B) S = 0.5D0 + (0.5D0-A) Z = DLOG(Y) T = Y - S*Z IF (B.LT.0.15D0) GO TO 20 XN = Y - S*DLOG(T) - DLOG(1.0D0+S/ (T+1.0D0)) GO TO 220 20 IF (B.LE.0.01D0) GO TO 30 U = ((T+2.0D0* (3.0D0-A))*T+ (2.0D0-A)* (3.0D0-A))/ + ((T+ (5.0D0-A))*T+2.0D0) XN = Y - S*DLOG(T) - DLOG(U) GO TO 220 30 C1 = -S*Z C2 = -S* (1.0D0+C1) C3 = S* ((0.5D0*C1+ (2.0D0-A))*C1+ (2.5D0-1.5D0*A)) C4 = -S* (((C1/3.0D0+ (2.5D0-1.5D0*A))*C1+ ((A-6.0D0)*A+7.0D0))* + C1+ ((11.0D0*A-46)*A+47.0D0)/6.0D0) C5 = -S* ((((-C1/4.0D0+ (11.0D0*A-17.0D0)/6.0D0)*C1+ ((-3.0D0*A+ + 13.0D0)*A-13.0D0))*C1+0.5D0* (((2.0D0*A-25.0D0)*A+72.0D0)*A- + 61.0D0))*C1+ (((25.0D0*A-195.0D0)*A+477.0D0)*A-379.0D0)/ + 12.0D0) XN = ((((C5/Y+C4)/Y+C3)/Y+C2)/Y+C1) + Y IF (A.GT.1.0D0) GO TO 220 IF (B.GT.BMIN(IOP)) GO TO 220 X = XN RETURN C 40 IF (B*Q.GT.1.D-8) GO TO 50 XN = EXP(- (Q/A+C)) GO TO 70 50 IF (P.LE.0.9D0) GO TO 60 XN = EXP((ALNREL(-Q)+GAMLN1(A))/A) GO TO 70 60 XN = EXP(DLOG(P*G)/A) 70 IF (XN.EQ.0.0D0) GO TO 310 T = 0.5D0 + (0.5D0-XN/ (A+1.0D0)) XN = XN/T GO TO 160 C C SELECTION OF THE INITIAL APPROXIMATION XN OF X C WHEN A .GT. 1 C 80 IF (Q.LE.0.5D0) GO TO 90 W = DLOG(P) GO TO 100 90 W = DLOG(Q) 100 T = SQRT(-2.0D0*W) S = T - (((A3*T+A2)*T+A1)*T+A0)/ ((((B4*T+B3)*T+B2)*T+B1)*T+1.0D0) IF (Q.GT.0.5D0) S = -S C RTA = SQRT(A) S2 = S*S XN = A + S*RTA + (S2-1.0D0)/3.0D0 + S* (S2-7.0D0)/ (36.0D0*RTA) - + ((3.0D0*S2+7.0D0)*S2-16.0D0)/ (810.0D0*A) + + S* ((9.0D0*S2+256.0D0)*S2-433.0D0)/ (38880.0D0*A*RTA) XN = DMAX1(XN,0.0D0) IF (A.LT.AMIN(IOP)) GO TO 110 X = XN D = 0.5D0 + (0.5D0-X/A) IF (ABS(D).LE.DMIN(IOP)) RETURN C 110 IF (P.LE.0.5D0) GO TO 130 IF (XN.LT.3.0D0*A) GO TO 220 Y = - (W+GAMLN(A)) D = DMAX1(2.0D0,A* (A-1.0D0)) IF (Y.LT.LN10*D) GO TO 120 S = 1.0D0 - A Z = DLOG(Y) GO TO 30 120 T = A - 1.0D0 XN = Y + T*DLOG(XN) - ALNREL(-T/ (XN+1.0D0)) XN = Y + T*DLOG(XN) - ALNREL(-T/ (XN+1.0D0)) GO TO 220 C 130 AP1 = A + 1.0D0 IF (XN.GT.0.70D0*AP1) GO TO 170 W = W + GAMLN(AP1) IF (XN.GT.0.15D0*AP1) GO TO 140 AP2 = A + 2.0D0 AP3 = A + 3.0D0 X = EXP((W+X)/A) X = EXP((W+X-DLOG(1.0D0+ (X/AP1)* (1.0D0+X/AP2)))/A) X = EXP((W+X-DLOG(1.0D0+ (X/AP1)* (1.0D0+X/AP2)))/A) X = EXP((W+X-DLOG(1.0D0+ (X/AP1)* (1.0D0+ (X/AP2)* (1.0D0+ + X/AP3))))/A) XN = X IF (XN.GT.1.D-2*AP1) GO TO 140 IF (XN.LE.EMIN(IOP)*AP1) RETURN GO TO 170 C 140 APN = AP1 T = XN/APN SUM = 1.0D0 + T 150 APN = APN + 1.0D0 T = T* (XN/APN) SUM = SUM + T IF (T.GT.1.D-4) GO TO 150 T = W - DLOG(SUM) XN = EXP((XN+T)/A) XN = XN* (1.0D0- (A*DLOG(XN)-XN-T)/ (A-XN)) GO TO 170 C C SCHRODER ITERATION USING P C 160 IF (P.GT.0.5D0) GO TO 220 170 IF (P.LE.1.D10*XMIN) GO TO 350 AM1 = (A-0.5D0) - 0.5D0 180 IF (A.LE.AMAX) GO TO 190 D = 0.5D0 + (0.5D0-XN/A) IF (ABS(D).LE.E2) GO TO 350 C 190 IF (IERR.GE.20) GO TO 330 IERR = IERR + 1 CALL GRATIO(A,XN,PN,QN,0) IF (PN.EQ.0.0D0 .OR. QN.EQ.0.0D0) GO TO 350 R = RCOMP(A,XN) IF (R.EQ.0.0D0) GO TO 350 T = (PN-P)/R W = 0.5D0* (AM1-XN) IF (ABS(T).LE.0.1D0 .AND. ABS(W*T).LE.0.1D0) GO TO 200 X = XN* (1.0D0-T) IF (X.LE.0.0D0) GO TO 340 D = ABS(T) GO TO 210 C 200 H = T* (1.0D0+W*T) X = XN* (1.0D0-H) IF (X.LE.0.0D0) GO TO 340 IF (ABS(W).GE.1.0D0 .AND. ABS(W)*T*T.LE.EPS) RETURN D = ABS(H) 210 XN = X IF (D.GT.TOL) GO TO 180 IF (D.LE.EPS) RETURN IF (ABS(P-PN).LE.TOL*P) RETURN GO TO 180 C C SCHRODER ITERATION USING Q C 220 IF (Q.LE.1.D10*XMIN) GO TO 350 AM1 = (A-0.5D0) - 0.5D0 230 IF (A.LE.AMAX) GO TO 240 D = 0.5D0 + (0.5D0-XN/A) IF (ABS(D).LE.E2) GO TO 350 C 240 IF (IERR.GE.20) GO TO 330 IERR = IERR + 1 CALL GRATIO(A,XN,PN,QN,0) IF (PN.EQ.0.0D0 .OR. QN.EQ.0.0D0) GO TO 350 R = RCOMP(A,XN) IF (R.EQ.0.0D0) GO TO 350 T = (Q-QN)/R W = 0.5D0* (AM1-XN) IF (ABS(T).LE.0.1D0 .AND. ABS(W*T).LE.0.1D0) GO TO 250 X = XN* (1.0D0-T) IF (X.LE.0.0D0) GO TO 340 D = ABS(T) GO TO 260 C 250 H = T* (1.0D0+W*T) X = XN* (1.0D0-H) IF (X.LE.0.0D0) GO TO 340 IF (ABS(W).GE.1.0D0 .AND. ABS(W)*T*T.LE.EPS) RETURN D = ABS(H) 260 XN = X IF (D.GT.TOL) GO TO 230 IF (D.LE.EPS) RETURN IF (ABS(Q-QN).LE.TOL*Q) RETURN GO TO 230 C C SPECIAL CASES C 270 X = XMAX RETURN C 280 IF (Q.LT.0.9D0) GO TO 290 X = -ALNREL(-P) RETURN 290 X = -DLOG(Q) RETURN C C ERROR RETURN C 300 IERR = -2 RETURN C 310 IERR = -3 RETURN C 320 IERR = -4 RETURN C 330 IERR = -6 RETURN C 340 IERR = -7 RETURN C 350 X = XN IERR = -8 RETURN C 360 X = XMAX IERR = -8 RETURN END C----------------------------------------------------------------------- DOUBLE PRECISION FUNCTION GAMLN(A) C----------------------------------------------------------------------- C EVALUATION OF LN(GAMMA(A)) FOR POSITIVE A C----------------------------------------------------------------------- C WRITTEN BY ALFRED H. MORRIS C NAVAL SURFACE WARFARE CENTER C DAHLGREN, VIRGINIA C-------------------------- C D = 0.5*(LN(2*PI) - 1) C-------------------------- C .. SCALAR ARGUMENTS .. DOUBLE PRECISION A C .. C .. LOCAL SCALARS .. DOUBLE PRECISION C0,C1,C2,C3,C4,C5,D,T,W INTEGER I,N C .. C .. EXTERNAL FUNCTIONS .. DOUBLE PRECISION GAMLN1 EXTERNAL GAMLN1 C .. C .. INTRINSIC FUNCTIONS .. INTRINSIC DLOG C .. C .. DATA STATEMENTS .. C-------------------------- DATA D/.418938533204673D0/ DATA C0/.833333333333333D-01/,C1/-.277777777760991D-02/, + C2/.793650666825390D-03/,C3/-.595202931351870D-03/, + C4/.837308034031215D-03/,C5/-.165322962780713D-02/ C .. C .. EXECUTABLE STATEMENTS .. C----------------------------------------------------------------------- IF (A.GT.0.8D0) GO TO 10 GAMLN = GAMLN1(A) - DLOG(A) RETURN 10 IF (A.GT.2.25D0) GO TO 20 T = (A-0.5D0) - 0.5D0 GAMLN = GAMLN1(T) RETURN C 20 IF (A.GE.10.0D0) GO TO 40 N = A - 1.25D0 T = A W = 1.0D0 DO 30 I = 1,N T = T - 1.0D0 W = T*W 30 CONTINUE GAMLN = GAMLN1(T-1.0D0) + DLOG(W) RETURN C 40 T = (1.0D0/A)**2 W = (((((C5*T+C4)*T+C3)*T+C2)*T+C1)*T+C0)/A GAMLN = (D+W) + (A-0.5D0)* (DLOG(A)-1.0D0) END C--------------------------------------------------------------------- DOUBLE PRECISION FUNCTION GAMLN1(A) C----------------------------------------------------------------------- C EVALUATION OF LN(GAMMA(1 + A)) FOR -0.2 .LE. A .LE. 1.25 C----------------------------------------------------------------------- C .. SCALAR ARGUMENTS .. DOUBLE PRECISION A C .. C .. LOCAL SCALARS .. DOUBLE PRECISION P0,P1,P2,P3,P4,P5,P6,Q1,Q2,Q3,Q4,Q5,Q6,R0,R1,R2, + R3,R4,R5,S1,S2,S3,S4,S5,W,X C .. C .. DATA STATEMENTS .. C---------------------- DATA P0/.577215664901533D+00/,P1/.844203922187225D+00/, + P2/-.168860593646662D+00/,P3/-.780427615533591D+00/, + P4/-.402055799310489D+00/,P5/-.673562214325671D-01/, + P6/-.271935708322958D-02/ DATA Q1/.288743195473681D+01/,Q2/.312755088914843D+01/, + Q3/.156875193295039D+01/,Q4/.361951990101499D+00/, + Q5/.325038868253937D-01/,Q6/.667465618796164D-03/ DATA R0/.422784335098467D+00/,R1/.848044614534529D+00/, + R2/.565221050691933D+00/,R3/.156513060486551D+00/, + R4/.170502484022650D-01/,R5/.497958207639485D-03/ DATA S1/.124313399877507D+01/,S2/.548042109832463D+00/, + S3/.101552187439830D+00/,S4/.713309612391000D-02/, + S5/.116165475989616D-03/ C .. C .. EXECUTABLE STATEMENTS .. C---------------------- IF (A.GE.0.6D0) GO TO 10 W = ((((((P6*A+P5)*A+P4)*A+P3)*A+P2)*A+P1)*A+P0)/ + ((((((Q6*A+Q5)*A+Q4)*A+Q3)*A+Q2)*A+Q1)*A+1.0D0) GAMLN1 = -A*W RETURN C 10 X = (A-0.5D0) - 0.5D0 W = (((((R5*X+R4)*X+R3)*X+R2)*X+R1)*X+R0)/ + (((((S5*X+S4)*X+S3)*X+S2)*X+S1)*X+1.0D0) GAMLN1 = X*W RETURN END C----------------------------------------------------------- DOUBLE PRECISION FUNCTION GAMMA(A) C----------------------------------------------------------------------- C C EVALUATION OF THE GAMMA FUNCTION FOR REAL ARGUMENTS C C ----------- C C GAMMA(A) IS ASSIGNED THE VALUE 0 WHEN THE GAMMA FUNCTION CANNOT C BE COMPUTED. C C----------------------------------------------------------------------- C WRITTEN BY ALFRED H. MORRIS, JR. C NAVAL SURFACE WEAPONS CENTER C DAHLGREN, VIRGINIA C----------------------------------------------------------------------- C .. SCALAR ARGUMENTS .. DOUBLE PRECISION A C .. C .. LOCAL SCALARS .. DOUBLE PRECISION BOT,D,G,LNX,PI,R1,R2,R3,R4,R5,S,T,TOP,W,X,Z INTEGER I,J,M,N C .. C .. LOCAL ARRAYS .. DOUBLE PRECISION P(7),Q(7) C .. C .. EXTERNAL FUNCTIONS .. DOUBLE PRECISION EXPARG,SPMPAR EXTERNAL EXPARG,SPMPAR C .. C .. INTRINSIC FUNCTIONS .. INTRINSIC ABS,DBLE,DLOG,EXP,INT,MOD,SIN C .. C .. DATA STATEMENTS .. C-------------------------- C D = 0.5*(LN(2*PI) - 1) C-------------------------- C-------------------------- C-------------------------- DATA PI/3.1415926535898D0/ DATA D/.41893853320467274178D0/ DATA P(1)/.539637273585445D-03/,P(2)/.261939260042690D-02/, + P(3)/.204493667594920D-01/,P(4)/.730981088720487D-01/, + P(5)/.279648642639792D+00/,P(6)/.553413866010467D+00/, + P(7)/1.0D0/ DATA Q(1)/-.832979206704073D-03/,Q(2)/.470059485860584D-02/, + Q(3)/.225211131035340D-01/,Q(4)/-.170458969313360D+00/, + Q(5)/-.567902761974940D-01/,Q(6)/.113062953091122D+01/, + Q(7)/1.0D0/ DATA R1/.820756370353826D-03/,R2/-.595156336428591D-03/, + R3/.793650663183693D-03/,R4/-.277777777770481D-02/, + R5/.833333333333333D-01/ C .. C .. EXECUTABLE STATEMENTS .. C-------------------------- GAMMA = 0.0D0 X = A IF (ABS(A).GE.15.0D0) GO TO 110 C----------------------------------------------------------------------- C EVALUATION OF GAMMA(A) FOR ABS(A) .LT. 15 C----------------------------------------------------------------------- T = 1.0D0 M = INT(A) - 1 C C LET T BE THE PRODUCT OF A-J WHEN A .GE. 2 C IF (M) 40,30,10 10 DO 20 J = 1,M X = X - 1.0D0 T = X*T 20 CONTINUE 30 X = X - 1.0D0 GO TO 80 C C LET T BE THE PRODUCT OF A+J WHEN A .LT. 1 C 40 T = A IF (A.GT.0.0D0) GO TO 70 M = -M - 1 IF (M.EQ.0) GO TO 60 DO 50 J = 1,M X = X + 1.0D0 T = X*T 50 CONTINUE 60 X = (X+0.5D0) + 0.5D0 T = X*T IF (T.EQ.0.0D0) RETURN C 70 CONTINUE C C THE FOLLOWING CODE CHECKS IF 1/T CAN OVERFLOW. THIS C CODE MAY BE OMITTED IF DESIRED. C IF (ABS(T).GE.1.D-30) GO TO 80 IF (ABS(T)*SPMPAR(3).LE.1.0001D0) RETURN GAMMA = 1.0D0/T RETURN C C COMPUTE GAMMA(1 + X) FOR 0 .LE. X .LT. 1 C 80 TOP = P(1) BOT = Q(1) DO 90 I = 2,7 TOP = P(I) + X*TOP BOT = Q(I) + X*BOT 90 CONTINUE GAMMA = TOP/BOT C C TERMINATION C IF (A.LT.1.0D0) GO TO 100 GAMMA = GAMMA*T RETURN 100 GAMMA = GAMMA/T RETURN C----------------------------------------------------------------------- C EVALUATION OF GAMMA(A) FOR ABS(A) .GE. 15 C----------------------------------------------------------------------- 110 IF (ABS(A).GE.1.D3) RETURN IF (A.GT.0.0D0) GO TO 120 X = -A N = X T = X - N IF (T.GT.0.9D0) T = 1.0D0 - T S = SIN(PI*T)/PI IF (MOD(N,2).EQ.0) S = -S IF (S.EQ.0.0D0) RETURN C C COMPUTE THE MODIFIED ASYMPTOTIC SUM C 120 T = 1.0D0/ (X*X) G = ((((R1*T+R2)*T+R3)*T+R4)*T+R5)/X C C ONE MAY REPLACE THE NEXT STATEMENT WITH LNX = ALOG(X) C BUT LESS ACCURACY WILL NORMALLY BE OBTAINED. C LNX = DLOG(X) C C FINAL ASSEMBLY C Z = X G = (D+G) + (Z-0.5D0)* (LNX-1.D0) W = G T = G - DBLE(W) IF (W.GT.0.99999D0*EXPARG(0)) RETURN GAMMA = EXP(W)* (1.0D0+T) IF (A.LT.0.0D0) GAMMA = (1.0D0/ (GAMMA*S))/X RETURN END C---------------------------------------------------- SUBROUTINE GRAT1(A,X,R,P,Q,EPS) C .. SCALAR ARGUMENTS .. DOUBLE PRECISION A,EPS,P,Q,R,X C .. C .. LOCAL SCALARS .. DOUBLE PRECISION A2N,A2NM1,AM0,AN,AN0,B2N,B2NM1,C,CMA,G,H,J,L,SUM, + T,TOL,W,Z C .. C .. EXTERNAL FUNCTIONS .. DOUBLE PRECISION ERF,ERFC1,GAM1,REXP EXTERNAL ERF,ERFC1,GAM1,REXP C .. C .. INTRINSIC FUNCTIONS .. INTRINSIC ABS,DLOG,EXP,SQRT C .. C .. EXECUTABLE STATEMENTS .. C----------------------------------------------------------------------- C EVALUATION OF THE INCOMPLETE GAMMA RATIO FUNCTIONS C P(A,X) AND Q(A,X) C C IT IS ASSUMED THAT A .LE. 1. EPS IS THE TOLERANCE TO BE USED. C THE INPUT ARGUMENT R HAS THE VALUE E**(-X)*X**A/GAMMA(A). C----------------------------------------------------------------------- IF (A*X.EQ.0.0D0) GO TO 120 IF (A.EQ.0.5D0) GO TO 100 IF (X.LT.1.1D0) GO TO 10 GO TO 60 C C TAYLOR SERIES FOR P(A,X)/X**A C 10 AN = 3.0D0 C = X SUM = X/ (A+3.0D0) TOL = 0.1D0*EPS/ (A+1.0D0) 20 AN = AN + 1.0D0 C = -C* (X/AN) T = C/ (A+AN) SUM = SUM + T IF (ABS(T).GT.TOL) GO TO 20 J = A*X* ((SUM/6.0D0-0.5D0/ (A+2.0D0))*X+1.0D0/ (A+1.0D0)) C Z = A*DLOG(X) H = GAM1(A) G = 1.0D0 + H IF (X.LT.0.25D0) GO TO 30 IF (A.LT.X/2.59D0) GO TO 50 GO TO 40 30 IF (Z.GT.-.13394D0) GO TO 50 C 40 W = EXP(Z) P = W*G* (0.5D0+ (0.5D0-J)) Q = 0.5D0 + (0.5D0-P) RETURN C 50 L = REXP(Z) W = 0.5D0 + (0.5D0+L) Q = (W*J-L)*G - H IF (Q.LT.0.0D0) GO TO 90 P = 0.5D0 + (0.5D0-Q) RETURN C C CONTINUED FRACTION EXPANSION C 60 A2NM1 = 1.0D0 A2N = 1.0D0 B2NM1 = X B2N = X + (1.0D0-A) C = 1.0D0 70 A2NM1 = X*A2N + C*A2NM1 B2NM1 = X*B2N + C*B2NM1 AM0 = A2NM1/B2NM1 C = C + 1.0D0 CMA = C - A A2N = A2NM1 + CMA*A2N B2N = B2NM1 + CMA*B2N AN0 = A2N/B2N IF (ABS(AN0-AM0).GE.EPS*AN0) GO TO 70 Q = R*AN0 P = 0.5D0 + (0.5D0-Q) RETURN C C SPECIAL CASES C 80 P = 0.0D0 Q = 1.0D0 RETURN C 90 P = 1.0D0 Q = 0.0D0 RETURN C 100 IF (X.GE.0.25D0) GO TO 110 P = ERF(SQRT(X)) Q = 0.5D0 + (0.5D0-P) RETURN 110 Q = ERFC1(0,SQRT(X)) P = 0.5D0 + (0.5D0-Q) RETURN C 120 IF (X.LE.A) GO TO 80 GO TO 90 END C-------------------------------------------------------- SUBROUTINE GRATIO(A,X,ANS,QANS,IND) C ---------------------------------------------------------------------- C EVALUATION OF THE INCOMPLETE GAMMA RATIO FUNCTIONS C P(A,X) AND Q(A,X) C C ---------- C C IT IS ASSUMED THAT A AND X ARE NONNEGATIVE, WHERE A AND X C ARE NOT BOTH 0. C C ANS AND QANS ARE VARIABLES. GRATIO ASSIGNS ANS THE VALUE C P(A,X) AND QANS THE VALUE Q(A,X). IND MAY BE ANY INTEGER. C IF IND = 0 THEN THE USER IS REQUESTING AS MUCH ACCURACY AS C POSSIBLE (UP TO 14 SIGNIFICANT DIGITS). OTHERWISE, IF C IND = 1 THEN ACCURACY IS REQUESTED TO WITHIN 1 UNIT OF THE C 6-TH SIGNIFICANT DIGIT, AND IF IND .NE. 0,1 THEN ACCURACY C IS REQUESTED TO WITHIN 1 UNIT OF THE 3RD SIGNIFICANT DIGIT. C C ERROR RETURN ... C ANS IS ASSIGNED THE VALUE 2 WHEN A OR X IS NEGATIVE, C WHEN A*X = 0, OR WHEN P(A,X) AND Q(A,X) ARE INDETERMINANT. C P(A,X) AND Q(A,X) ARE COMPUTATIONALLY INDETERMINANT WHEN C X IS EXCEEDINGLY CLOSE TO A AND A IS EXTREMELY LARGE. C ---------------------------------------------------------------------- C WRITTEN BY ALFRED H. MORRIS, JR. C NAVAL SURFACE WEAPONS CENTER C DAHLGREN, VIRGINIA C -------------------- C .. SCALAR ARGUMENTS .. DOUBLE PRECISION A,ANS,QANS,X INTEGER IND C .. C .. LOCAL SCALARS .. DOUBLE PRECISION A2N,A2NM1,ACC,ALOG10,AM0,AMN,AN,AN0,APN,B2N, + B2NM1,C,C0,C1,C2,C3,C4,C5,C6,CMA,D10,D20,D30,D40, + D50,D60,D70,E,E0,G,H,J,L,R,RT2PIN,RTA,RTPI,RTX,S, + SUM,T,T1,THIRD,TOL,TWOA,U,W,X0,Y,Z INTEGER I,IOP,M,MAX,N C .. C .. LOCAL ARRAYS .. DOUBLE PRECISION ACC0(3),BIG(3),D0(13),D1(12),D2(10),D3(8),D4(6), + D5(4),D6(2),E00(3),WK(20),X00(3) C .. C .. EXTERNAL FUNCTIONS .. DOUBLE PRECISION ERF,ERFC1,GAM1,GAMMA,REXP,RLOG,SPMPAR EXTERNAL ERF,ERFC1,GAM1,GAMMA,REXP,RLOG,SPMPAR C .. C .. INTRINSIC FUNCTIONS .. INTRINSIC ABS,DBLE,DLOG,DMAX1,EXP,INT,SQRT C .. C .. DATA STATEMENTS .. C -------------------- C -------------------- C ALOG10 = LN(10) C RT2PIN = 1/SQRT(2*PI) C RTPI = SQRT(PI) C -------------------- C -------------------- C -------------------- C -------------------- C -------------------- C -------------------- C -------------------- C -------------------- C -------------------- DATA ACC0(1)/5.D-15/,ACC0(2)/5.D-7/,ACC0(3)/5.D-4/ DATA BIG(1)/20.0D0/,BIG(2)/14.0D0/,BIG(3)/10.0D0/ DATA E00(1)/.25D-3/,E00(2)/.25D-1/,E00(3)/.14D0/ DATA X00(1)/31.0D0/,X00(2)/17.0D0/,X00(3)/9.7D0/ DATA ALOG10/2.30258509299405D0/ DATA RT2PIN/.398942280401433D0/ DATA RTPI/1.77245385090552D0/ DATA THIRD/.333333333333333D0/ DATA D0(1)/.833333333333333D-01/,D0(2)/-.148148148148148D-01/, + D0(3)/.115740740740741D-02/,D0(4)/.352733686067019D-03/, + D0(5)/-.178755144032922D-03/,D0(6)/.391926317852244D-04/, + D0(7)/-.218544851067999D-05/,D0(8)/-.185406221071516D-05/, + D0(9)/.829671134095309D-06/,D0(10)/-.176659527368261D-06/, + D0(11)/.670785354340150D-08/,D0(12)/.102618097842403D-07/, + D0(13)/-.438203601845335D-08/ DATA D10/-.185185185185185D-02/,D1(1)/-.347222222222222D-02/, + D1(2)/.264550264550265D-02/,D1(3)/-.990226337448560D-03/, + D1(4)/.205761316872428D-03/,D1(5)/-.401877572016461D-06/, + D1(6)/-.180985503344900D-04/,D1(7)/.764916091608111D-05/, + D1(8)/-.161209008945634D-05/,D1(9)/.464712780280743D-08/, + D1(10)/.137863344691572D-06/,D1(11)/-.575254560351770D-07/, + D1(12)/.119516285997781D-07/ DATA D20/.413359788359788D-02/,D2(1)/-.268132716049383D-02/, + D2(2)/.771604938271605D-03/,D2(3)/.200938786008230D-05/, + D2(4)/-.107366532263652D-03/,D2(5)/.529234488291201D-04/, + D2(6)/-.127606351886187D-04/,D2(7)/.342357873409614D-07/, + D2(8)/.137219573090629D-05/,D2(9)/-.629899213838006D-06/, + D2(10)/.142806142060642D-06/ DATA D30/.649434156378601D-03/,D3(1)/.229472093621399D-03/, + D3(2)/-.469189494395256D-03/,D3(3)/.267720632062839D-03/, + D3(4)/-.756180167188398D-04/,D3(5)/-.239650511386730D-06/, + D3(6)/.110826541153473D-04/,D3(7)/-.567495282699160D-05/, + D3(8)/.142309007324359D-05/ DATA D40/-.861888290916712D-03/,D4(1)/.784039221720067D-03/, + D4(2)/-.299072480303190D-03/,D4(3)/-.146384525788434D-05/, + D4(4)/.664149821546512D-04/,D4(5)/-.396836504717943D-04/, + D4(6)/.113757269706784D-04/ DATA D50/-.336798553366358D-03/,D5(1)/-.697281375836586D-04/, + D5(2)/.277275324495939D-03/,D5(3)/-.199325705161888D-03/, + D5(4)/.679778047793721D-04/ DATA D60/.531307936463992D-03/,D6(1)/-.592166437353694D-03/, + D6(2)/.270878209671804D-03/ DATA D70/.344367606892378D-03/ C .. C .. EXECUTABLE STATEMENTS .. C -------------------- C ****** E IS A MACHINE DEPENDENT CONSTANT. E IS THE SMALLEST C FLOATING POINT NUMBER FOR WHICH 1.0 + E .GT. 1.0 . C E = SPMPAR(1) C C -------------------- IF (A.LT.0.0D0 .OR. X.LT.0.0D0) GO TO 430 IF (A.EQ.0.0D0 .AND. X.EQ.0.0D0) GO TO 430 IF (A*X.EQ.0.0D0) GO TO 420 C IOP = IND + 1 IF (IOP.NE.1 .AND. IOP.NE.2) IOP = 3 ACC = DMAX1(ACC0(IOP),E) E0 = E00(IOP) X0 = X00(IOP) C C SELECT THE APPROPRIATE ALGORITHM C IF (A.GE.1.0D0) GO TO 10 IF (A.EQ.0.5D0) GO TO 390 IF (X.LT.1.1D0) GO TO 160 T1 = A*DLOG(X) - X U = A*EXP(T1) IF (U.EQ.0.0D0) GO TO 380 R = U* (1.0D0+GAM1(A)) GO TO 250 C 10 IF (A.GE.BIG(IOP)) GO TO 30 IF (A.GT.X .OR. X.GE.X0) GO TO 20 TWOA = A + A M = INT(TWOA) IF (TWOA.NE.DBLE(M)) GO TO 20 I = M/2 IF (A.EQ.DBLE(I)) GO TO 210 GO TO 220 20 T1 = A*DLOG(X) - X R = EXP(T1)/GAMMA(A) GO TO 40 C 30 L = X/A IF (L.EQ.0.0D0) GO TO 370 S = 0.5D0 + (0.5D0-L) Z = RLOG(L) IF (Z.GE.700.0D0/A) GO TO 410 Y = A*Z RTA = SQRT(A) IF (ABS(S).LE.E0/RTA) GO TO 330 IF (ABS(S).LE.0.4D0) GO TO 270 C T = (1.0D0/A)**2 T1 = (((0.75D0*T-1.0D0)*T+3.5D0)*T-105.0D0)/ (A*1260.0D0) T1 = T1 - Y R = RT2PIN*RTA*EXP(T1) C 40 IF (R.EQ.0.0D0) GO TO 420 IF (X.LE.DMAX1(A,ALOG10)) GO TO 50 IF (X.LT.X0) GO TO 250 GO TO 100 C C TAYLOR SERIES FOR P/R C 50 APN = A + 1.0D0 T = X/APN WK(1) = T DO 60 N = 2,20 APN = APN + 1.0D0 T = T* (X/APN) IF (T.LE.1.D-3) GO TO 70 WK(N) = T 60 CONTINUE N = 20 C 70 SUM = T TOL = 0.5D0*ACC 80 APN = APN + 1.0D0 T = T* (X/APN) SUM = SUM + T IF (T.GT.TOL) GO TO 80 C MAX = N - 1 DO 90 M = 1,MAX N = N - 1 SUM = SUM + WK(N) 90 CONTINUE ANS = (R/A)* (1.0D0+SUM) QANS = 0.5D0 + (0.5D0-ANS) RETURN C C ASYMPTOTIC EXPANSION C 100 AMN = A - 1.0D0 T = AMN/X WK(1) = T DO 110 N = 2,20 AMN = AMN - 1.0D0 T = T* (AMN/X) IF (ABS(T).LE.1.D-3) GO TO 120 WK(N) = T 110 CONTINUE N = 20 C 120 SUM = T 130 IF (ABS(T).LE.ACC) GO TO 140 AMN = AMN - 1.0D0 T = T* (AMN/X) SUM = SUM + T GO TO 130 C 140 MAX = N - 1 DO 150 M = 1,MAX N = N - 1 SUM = SUM + WK(N) 150 CONTINUE QANS = (R/X)* (1.0D0+SUM) ANS = 0.5D0 + (0.5D0-QANS) RETURN C C TAYLOR SERIES FOR P(A,X)/X**A C 160 AN = 3.0D0 C = X SUM = X/ (A+3.0D0) TOL = 3.0D0*ACC/ (A+1.0D0) 170 AN = AN + 1.0D0 C = -C* (X/AN) T = C/ (A+AN) SUM = SUM + T IF (ABS(T).GT.TOL) GO TO 170 J = A*X* ((SUM/6.0D0-0.5D0/ (A+2.0D0))*X+1.0D0/ (A+1.0D0)) C Z = A*DLOG(X) H = GAM1(A) G = 1.0D0 + H IF (X.LT.0.25D0) GO TO 180 IF (A.LT.X/2.59D0) GO TO 200 GO TO 190 180 IF (Z.GT.-.13394D0) GO TO 200 C 190 W = EXP(Z) ANS = W*G* (0.5D0+ (0.5D0-J)) QANS = 0.5D0 + (0.5D0-ANS) RETURN C 200 L = REXP(Z) W = 0.5D0 + (0.5D0+L) QANS = (W*J-L)*G - H IF (QANS.LT.0.0D0) GO TO 380 ANS = 0.5D0 + (0.5D0-QANS) RETURN C C FINITE SUMS FOR Q WHEN A .GE. 1 C AND 2*A IS AN INTEGER C 210 SUM = EXP(-X) T = SUM N = 1 C = 0.0D0 GO TO 230 C 220 RTX = SQRT(X) SUM = ERFC1(0,RTX) T = EXP(-X)/ (RTPI*RTX) N = 0 C = -0.5D0 C 230 IF (N.EQ.I) GO TO 240 N = N + 1 C = C + 1.0D0 T = (X*T)/C SUM = SUM + T GO TO 230 240 QANS = SUM ANS = 0.5D0 + (0.5D0-QANS) RETURN C C CONTINUED FRACTION EXPANSION C 250 TOL = DMAX1(5.0D0*E,ACC) A2NM1 = 1.0D0 A2N = 1.0D0 B2NM1 = X B2N = X + (1.0D0-A) C = 1.0D0 260 A2NM1 = X*A2N + C*A2NM1 B2NM1 = X*B2N + C*B2NM1 AM0 = A2NM1/B2NM1 C = C + 1.0D0 CMA = C - A A2N = A2NM1 + CMA*A2N B2N = B2NM1 + CMA*B2N AN0 = A2N/B2N IF (ABS(AN0-AM0).GE.TOL*AN0) GO TO 260 C QANS = R*AN0 ANS = 0.5D0 + (0.5D0-QANS) RETURN C C GENERAL TEMME EXPANSION C 270 IF (ABS(S).LE.2.0D0*E .AND. A*E*E.GT.3.28D-3) GO TO 430 C = EXP(-Y) W = 0.5D0*ERFC1(1,SQRT(Y)) U = 1.0D0/A Z = SQRT(Z+Z) IF (L.LT.1.0D0) Z = -Z IF (IOP-2) 280,290,300 C 280 IF (ABS(S).LE.1.D-3) GO TO 340 C0 = ((((((((((((D0(13)*Z+D0(12))*Z+D0(11))*Z+D0(10))*Z+D0(9))*Z+ + D0(8))*Z+D0(7))*Z+D0(6))*Z+D0(5))*Z+D0(4))*Z+D0(3))*Z+D0(2))* + Z+D0(1))*Z - THIRD C1 = (((((((((((D1(12)*Z+D1(11))*Z+D1(10))*Z+D1(9))*Z+D1(8))*Z+ + D1(7))*Z+D1(6))*Z+D1(5))*Z+D1(4))*Z+D1(3))*Z+D1(2))*Z+D1(1))* + Z + D10 C2 = (((((((((D2(10)*Z+D2(9))*Z+D2(8))*Z+D2(7))*Z+D2(6))*Z+ + D2(5))*Z+D2(4))*Z+D2(3))*Z+D2(2))*Z+D2(1))*Z + D20 C3 = (((((((D3(8)*Z+D3(7))*Z+D3(6))*Z+D3(5))*Z+D3(4))*Z+D3(3))*Z+ + D3(2))*Z+D3(1))*Z + D30 C4 = (((((D4(6)*Z+D4(5))*Z+D4(4))*Z+D4(3))*Z+D4(2))*Z+D4(1))*Z + + D40 C5 = (((D5(4)*Z+D5(3))*Z+D5(2))*Z+D5(1))*Z + D50 C6 = (D6(2)*Z+D6(1))*Z + D60 T = ((((((D70*U+C6)*U+C5)*U+C4)*U+C3)*U+C2)*U+C1)*U + C0 GO TO 310 C 290 C0 = (((((D0(6)*Z+D0(5))*Z+D0(4))*Z+D0(3))*Z+D0(2))*Z+D0(1))*Z - + THIRD C1 = (((D1(4)*Z+D1(3))*Z+D1(2))*Z+D1(1))*Z + D10 C2 = D2(1)*Z + D20 T = (C2*U+C1)*U + C0 GO TO 310 C 300 T = ((D0(3)*Z+D0(2))*Z+D0(1))*Z - THIRD C 310 IF (L.LT.1.0D0) GO TO 320 QANS = C* (W+RT2PIN*T/RTA) ANS = 0.5D0 + (0.5D0-QANS) RETURN 320 ANS = C* (W-RT2PIN*T/RTA) QANS = 0.5D0 + (0.5D0-ANS) RETURN C C TEMME EXPANSION FOR L = 1 C 330 IF (A*E*E.GT.3.28D-3) GO TO 430 C = 0.5D0 + (0.5D0-Y) W = (0.5D0-SQRT(Y)* (0.5D0+ (0.5D0-Y/3.0D0))/RTPI)/C U = 1.0D0/A Z = SQRT(Z+Z) IF (L.LT.1.0D0) Z = -Z IF (IOP-2) 340,350,360 C 340 C0 = ((((((D0(7)*Z+D0(6))*Z+D0(5))*Z+D0(4))*Z+D0(3))*Z+D0(2))*Z+ + D0(1))*Z - THIRD C1 = (((((D1(6)*Z+D1(5))*Z+D1(4))*Z+D1(3))*Z+D1(2))*Z+D1(1))*Z + + D10 C2 = ((((D2(5)*Z+D2(4))*Z+D2(3))*Z+D2(2))*Z+D2(1))*Z + D20 C3 = (((D3(4)*Z+D3(3))*Z+D3(2))*Z+D3(1))*Z + D30 C4 = (D4(2)*Z+D4(1))*Z + D40 C5 = (D5(2)*Z+D5(1))*Z + D50 C6 = D6(1)*Z + D60 T = ((((((D70*U+C6)*U+C5)*U+C4)*U+C3)*U+C2)*U+C1)*U + C0 GO TO 310 C 350 C0 = (D0(2)*Z+D0(1))*Z - THIRD C1 = D1(1)*Z + D10 T = (D20*U+C1)*U + C0 GO TO 310 C 360 T = D0(1)*Z - THIRD GO TO 310 C C SPECIAL CASES C 370 ANS = 0.0D0 QANS = 1.0D0 RETURN C 380 ANS = 1.0D0 QANS = 0.0D0 RETURN C 390 IF (X.GE.0.25D0) GO TO 400 ANS = ERF(SQRT(X)) QANS = 0.5D0 + (0.5D0-ANS) RETURN 400 QANS = ERFC1(0,SQRT(X)) ANS = 0.5D0 + (0.5D0-QANS) RETURN C 410 IF (ABS(S).LE.2.0D0*E) GO TO 430 420 IF (X.LE.A) GO TO 370 GO TO 380 C C ERROR RETURN C 430 ANS = 2.0D0 RETURN END C--------------------------------------------------------------------- DOUBLE PRECISION FUNCTION GSUMLN(A,B) C----------------------------------------------------------------------- C EVALUATION OF THE FUNCTION LN(GAMMA(A + B)) C FOR 1 .LE. A .LE. 2 AND 1 .LE. B .LE. 2 C----------------------------------------------------------------------- C .. SCALAR ARGUMENTS .. DOUBLE PRECISION A,B C .. C .. LOCAL SCALARS .. DOUBLE PRECISION X C .. C .. EXTERNAL FUNCTIONS .. DOUBLE PRECISION ALNREL,GAMLN1 EXTERNAL ALNREL,GAMLN1 C .. C .. INTRINSIC FUNCTIONS .. INTRINSIC DBLE,DLOG C .. C .. EXECUTABLE STATEMENTS .. X = DBLE(A) + DBLE(B) - 2.D0 IF (X.GT.0.25D0) GO TO 10 GSUMLN = GAMLN1(1.0D0+X) RETURN 10 IF (X.GT.1.25D0) GO TO 20 GSUMLN = GAMLN1(X) + ALNREL(X) RETURN 20 GSUMLN = GAMLN1(X-1.0D0) + DLOG(X* (1.0D0+X)) RETURN END C------------------------------------------------------------------ DOUBLE PRECISION FUNCTION INVNOR(P) C C********************************************************************** C C DOUBLE PRECISION FUNCTION INVNOR(P) C NORMAL DISTRIBUTION INVERSE C C C FUNCTION C C C RETURNS X SUCH THAT CUMNOR(X) = P, I.E., THE INTEGRAL FROM - C INFINITY TO X OF (1/SQRT(2*PI)) EXP(-U*U) DU IS P C C C ARGUMENTS C C C P --> THE PROBABILITY WHOSE NORMAL DEVIATE IS SOUGHT. C P IS DOUBLE PRECISION C C C METHOD C C C THE RATIONAL FUNCTION ON PAGE 95 OF KENNEDY AND GENTLE, C STATISTICAL COMPUTING, MARCEL DEKKER, NY , 1980. C C C NOTE C C C IF P .LT. 1.0E-20 THEN INVNOR RETURNS -10.0. C IF P .GE. 1.0 THEN INVNOR RETURNS 10.0. C C********************************************************************** C C .. SCALAR ARGUMENTS .. DOUBLE PRECISION P C .. C .. LOCAL SCALARS .. DOUBLE PRECISION SIGN,Y,Z C .. C .. LOCAL ARRAYS .. DOUBLE PRECISION XDEN(5),XNUM(5) C .. C .. EXTERNAL FUNCTIONS .. DOUBLE PRECISION DEVLPL EXTERNAL DEVLPL C .. C .. INTRINSIC FUNCTIONS .. INTRINSIC DBLE,DLOG,DSQRT C .. C .. DATA STATEMENTS .. DATA XNUM/-0.322232431088D0,-1.000000000000D0,-0.342242088547D0, + -0.204231210245D-1,-0.453642210148D-4/ DATA XDEN/0.993484626060D-1,0.588581570495D0,0.531103462366D0, + 0.103537752850D0,0.38560700634D-2/ C .. C .. EXECUTABLE STATEMENTS .. IF (.NOT. (P.LT.1.0D-20)) GO TO 10 INVNOR = -10.0D0 RETURN 10 IF (.NOT. (P.GE.1.0D0)) GO TO 20 INVNOR = 10.0D0 RETURN 20 IF (.NOT. (P.LE.0.5D0)) GO TO 30 SIGN = -1.0D0 Z = P GO TO 40 30 SIGN = 1.0D0 Z = 1.0D0 - P 40 Y = DSQRT(-2.0D0*DLOG(Z)) INVNOR = Y + DEVLPL(XNUM,5,Y)/DEVLPL(XDEN,5,Y) INVNOR = SIGN*INVNOR RETURN END C---------------------------------------------------------------- SUBROUTINE INVR(STATUS,X,FX,QLEFT,QHI) C********************************************************************** C C SUBROUTINE INVR(STATUS, X, FX, QLEFT, QHI) C BOUNDS THE ZERO OF THE FUNCTION AND INVOKES ZROR C REVERSE COMMUNICATION C C C FUNCTION C C C BOUNDS THE FUNCTION AND INVOKES ZROR TO PERFORM THE ZERO C FINDING. STINVR MUST HAVE BEEN CALLED BEFORE THIS ROUTINE C IN ORDER TO SET ITS PARAMETERS. C C C ARGUMENTS C C C STATUS <--> AT THE BEGINNING OF A ZERO FINDING PROBLEM, STATUS C SHOULD BE SET TO 0 AND INVR INVOKED. (THE VALUE C OF PARAMETERS OTHER THAN X WILL BE IGNORED ON THIS CAL C C WHEN INVR NEEDS THE FUNCTION EVALUATED, IT WILL SET C STATUS TO 1 AND RETURN. THE VALUE OF THE FUNCTION C SHOULD BE SET IN FX AND INVR AGAIN CALLED WITHOUT C CHANGING ANY OF ITS OTHER PARAMETERS. C C WHEN INVR HAS FINISHED WITHOUT ERROR, IT WILL RETURN C WITH STATUS 0. IN THAT CASE X IS APPROXIMATELY A ROOT C OF F(X). C C IF INVR CANNOT BOUND THE FUNCTION, IT RETURNS STATUS C -1 AND SETS QLEFT AND QHI. C INTEGER STATUS C C X <-- THE VALUE OF X AT WHICH F(X) IS TO BE EVALUATED. C DOUBLE PRECISION X C C FX --> THE VALUE OF F(X) CALCULATED WHEN INVR RETURNS WITH C STATUS = 1. C DOUBLE PRECISION FX C C QLEFT <-- DEFINED ONLY IF QMFINV RETURNS .FALSE. IN THAT C CASE IT IS .TRUE. IF THE STEPPING SEARCH TERMINATED C UNSUCESSFULLY AT SMALL. IF IT IS .FALSE. THE SEARCH C TERMINATED UNSUCESSFULLY AT BIG. C QLEFT IS LOGICAL C C QHI <-- DEFINED ONLY IF QMFINV RETURNS .FALSE. IN THAT C CASE IT IS .TRUE. IF F(X) .GT. Y AT THE TERMINATION C OF THE SEARCH AND .FALSE. IF F(X) .LT. Y AT THE C TERMINATION OF THE SEARCH. C QHI IS LOGICAL C C********************************************************************** C .. SCALAR ARGUMENTS .. DOUBLE PRECISION FX,X,ZABSST,ZABSTO,ZBIG,ZRELST,ZRELTO,ZSMALL, +ZSTPMU INTEGER STATUS LOGICAL QHI,QLEFT C .. C .. LOCAL SCALARS .. DOUBLE PRECISION ABSSTP,ABSTOL,BIG,FBIG,FSMALL,RELSTP,RELTOL, +SMALL,STEP,STPMUL,XHI,XLB,XLO,XSAVE,XUB,YY,ZX,ZY,ZZ INTEGER I99999 LOGICAL QBDD,QCOND,QDUM1,QDUM2,QINCR,QLIM,QOK,QUP C .. C .. EXTERNAL SUBROUTINES .. EXTERNAL STZROR,ZROR C .. C .. INTRINSIC FUNCTIONS .. INTRINSIC ABS,MAX,MIN C .. C .. STATEMENT FUNCTIONS .. LOGICAL QXMON C .. C .. SAVE STATEMENT .. SAVE C .. C .. STATEMENT FUNCTION DEFINITIONS .. QXMON(ZX,ZY,ZZ) = ZX .LE. ZY .AND. ZY .LE. ZZ C .. C .. EXECUTABLE STATEMENTS .. IF (STATUS.GT.0) GO TO 310 QCOND = .NOT. QXMON(SMALL,X,BIG) IF (QCOND) STOP ' SMALL, X, BIG NOT MONOTONE IN INVR' XSAVE = X C C SEE THAT SMALL AND BIG BOUND THE ZERO AND SET QINCR C X = SMALL C GET-FUNCTION-VALUE ASSIGN 10 TO I99999 GO TO 300 10 FSMALL = FX X = BIG C GET-FUNCTION-VALUE ASSIGN 20 TO I99999 GO TO 300 20 FBIG = FX QINCR = FBIG .GT. FSMALL IF (.NOT. (QINCR)) GO TO 50 IF (FSMALL.LE.0.0D0) GO TO 30 STATUS = -1 QLEFT = .TRUE. QHI = .TRUE. RETURN 30 IF (FBIG.GE.0.0D0) GO TO 40 STATUS = -1 QLEFT = .FALSE. QHI = .FALSE. RETURN 40 GO TO 80 50 IF (FSMALL.GE.0.0D0) GO TO 60 STATUS = -1 QLEFT = .TRUE. QHI = .FALSE. RETURN 60 IF (FBIG.LE.0.0D0) GO TO 70 STATUS = -1 QLEFT = .FALSE. QHI = .TRUE. RETURN 70 CONTINUE 80 X = XSAVE STEP =DMAX1(ABSSTP,RELSTP*DABS(X)) C YY = F(X) - Y C GET-FUNCTION-VALUE ASSIGN 90 TO I99999 GO TO 300 90 YY = FX IF (.NOT. (YY.EQ.0.0D0)) GO TO 100 STATUS = 0 QOK = .TRUE. RETURN 100 QUP = (QINCR .AND. (YY.LT.0.0D0)) .OR. + (.NOT.QINCR .AND. (YY.GT.0.0D0)) C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ C C HANDLE CASE IN WHICH WE MUST STEP HIGHER C C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ IF (.NOT. (QUP)) GO TO 170 XLB = XSAVE XUB = DMIN1(XLB+STEP,BIG) GO TO 120 110 IF (QCOND) GO TO 150 C YY = F(XUB) - Y 120 X = XUB C GET-FUNCTION-VALUE ASSIGN 130 TO I99999 GO TO 300 130 YY = FX QBDD = (QINCR .AND. (YY.GE.0.0D0)) .OR. + (.NOT.QINCR .AND. (YY.LE.0.0D0)) QLIM = XUB .GE. BIG QCOND = QBDD .OR. QLIM IF (QCOND) GO TO 140 STEP = STPMUL*STEP XLB = XUB XUB = DMIN1(XLB+STEP,BIG) 140 GO TO 110 150 IF (.NOT. (QLIM.AND..NOT.QBDD)) GO TO 160 STATUS = -1 QLEFT = .FALSE. QHI = .NOT. QINCR X = BIG RETURN 160 GO TO 240 C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ C C HANDLE CASE IN WHICH WE MUST STEP LOWER C C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 170 XUB = XSAVE XLB = DMAX1(XUB-STEP,SMALL) GO TO 190 180 IF (QCOND) GO TO 220 C YY = F(XLB) - Y 190 X = XLB C GET-FUNCTION-VALUE ASSIGN 200 TO I99999 GO TO 300 200 YY = FX QBDD = (QINCR .AND. (YY.LE.0.0D0)) .OR. + (.NOT.QINCR .AND. (YY.GE.0.0D0)) QLIM = XLB .LE. SMALL QCOND = QBDD .OR. QLIM IF (QCOND) GO TO 210 STEP = STPMUL*STEP XUB = XLB XLB = DMAX1(XUB-STEP,SMALL) 210 GO TO 180 220 IF (.NOT. (QLIM.AND..NOT.QBDD)) GO TO 230 STATUS = -1 QLEFT = .TRUE. QHI = QINCR X = SMALL RETURN 230 CONTINUE 240 CALL STZROR(XLB,XUB,ABSTOL,RELTOL) C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ C C IF WE REACH HERE, XLB AND XUB BOUND THE ZERO OF F. C C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ STATUS = 0 GO TO 260 250 IF (.NOT. (STATUS.EQ.1)) GO TO 290 260 CALL ZROR(STATUS,X,FX,XLO,XHI,QDUM1,QDUM2) IF (.NOT. (STATUS.EQ.1)) GO TO 280 C GET-FUNCTION-VALUE ASSIGN 270 TO I99999 GO TO 300 270 CONTINUE 280 GO TO 250 290 X = XLO C IF (QRZERO(F,Y,XLB,XUB,ABSTOL,RELTOL)) CONTINUE C X = XLB STATUS = 0 RETURN ENTRY STINVR(ZSMALL,ZBIG,ZABSST,ZRELST,ZSTPMU,ZABSTO,ZRELTO) C********************************************************************** C C SUBROUTINE STINVR( SMALL, BIG, ABSSTP, RELSTP, STPMUL, C + ABSTOL, RELTOL ) C SET INVERSE FINDER - REVERSE COMMUNICATION C C C FUNCTION C C C CONCISE DESCRIPTION - GIVEN A MONOTONE FUNCTION F FINDS X C SUCH THAT F(X) = Y. USES REVERSE COMMUNICATION -- SEE INVR. C THIS ROUTINE SETS QUANTITIES NEEDED BY INVR. C C MORE PRECISE DESCRIPTION OF INVR - C C F MUST BE A MONOTONE FUNCTION, THE RESULTS OF QMFINV ARE C OTHERWISE UNDEFINED. QINCR MUST BE .TRUE. IF F IS NON- C DECREASING AND .FALSE. IF F IS NON-INCREASING. C C QMFINV WILL RETURN .TRUE. IF AND ONLY IF F(SMALL) AND C F(BIG) BRACKET Y, I. E., C QINCR IS .TRUE. AND F(SMALL).LE.Y.LE.F(BIG) OR C QINCR IS .FALSE. AND F(BIG).LE.Y.LE.F(SMALL) C C IF QMFINV RETURNS .TRUE., THEN THE X RETURNED SATISFIES C THE FOLLOWING CONDITION. LET C TOL(X) = MAX(ABSTOL,RELTOL*ABS(X)) C THEN IF QINCR IS .TRUE., C F(X-TOL(X)) .LE. Y .LE. F(X+TOL(X)) C AND IF QINCR IS .FALSE. C F(X-TOL(X)) .GE. Y .GE. F(X+TOL(X)) C C C ARGUMENTS C C C SMALL --> THE LEFT ENDPOINT OF THE INTERVAL TO BE C SEARCHED FOR A SOLUTION. C SMALL IS DOUBLE PRECISION C C BIG --> THE RIGHT ENDPOINT OF THE INTERVAL TO BE C SEARCHED FOR A SOLUTION. C BIG IS DOUBLE PRECISION C C ABSSTP, RELSTP --> THE INITIAL STEP SIZE IN THE SEARCH C IS MAX(ABSSTP,RELSTP*ABS(X)). SEE ALGORITHM. C ABSSTP IS DOUBLE PRECISION C RELSTP IS DOUBLE PRECISION C C STPMUL --> WHEN A STEP DOESN'T BOUND THE ZERO, THE STEP C SIZE IS MULTIPLIED BY STPMUL AND ANOTHER STEP C TAKEN. A POPULAR VALUE IS 2.0 C DOUBLE PRECISION STPMUL C C ABSTOL, RELTOL --> TWO NUMBERS THAT DETERMINE THE ACCURACY C OF THE SOLUTION. SEE FUNCTION FOR A PRECISE DEFINITION. C ABSTOL IS DOUBLE PRECISION C RELTOL IS DOUBLE PRECISION C C C METHOD C C C COMPARES F(X) WITH Y FOR THE INPUT VALUE OF X THEN USES QINCR C TO DETERMINE WHETHER TO STEP LEFT OR RIGHT TO BOUND THE C DESIRED X. THE INITIAL STEP SIZE IS C MAX(ABSSTP,RELSTP*ABS(S)) FOR THE INPUT VALUE OF X. C ITERATIVELY STEPS RIGHT OR LEFT UNTIL IT BOUNDS X. C AT EACH STEP WHICH DOESN'T BOUND X, THE STEP SIZE IS DOUBLED. C THE ROUTINE IS CAREFUL NEVER TO STEP BEYOND SMALL OR BIG. IF C IT HASN'T BOUNDED X AT SMALL OR BIG, QMFINV RETURNS .FALSE. C AFTER SETTING QLEFT AND QHI. C C IF X IS SUCCESSFULLY BOUNDED THEN ALGORITHM R OF THE PAPER C 'TWO EFFICIENT ALGORITHMS WITH GUARANTEED CONVERGENCE FOR C FINDING A ZERO OF A FUNCTION' BY J. C. P. BUS AND C T. J. DEKKER IN ACM TRANSACTIONS ON MATHEMATICAL C SOFTWARE, VOLUME 1, NO. 4 PAGE 330 (DEC. '75) IS EMPLOYED C TO FIND THE ZERO OF THE FUNCTION F(X)-Y. THIS IS ROUTINE C QRZERO. C C********************************************************************** SMALL = ZSMALL BIG = ZBIG ABSSTP = ZABSST RELSTP = ZRELST STPMUL = ZSTPMU ABSTOL = ZABSTO RELTOL = ZRELTO RETURN STOP '*** EXECUTION FLOWING INTO FLECS PROCEDURES ***' C TO GET-FUNCTION-VALUE 300 STATUS = 1 RETURN 310 CONTINUE GO TO I99999 END C----------------------------------------------------------- INTEGER FUNCTION IPMPAR(I) C----------------------------------------------------------------------- C C IPMPAR PROVIDES THE INTEGER MACHINE CONSTANTS FOR THE COMPUTER C THAT IS USED. IT IS ASSUMED THAT THE ARGUMENT I IS AN INTEGER C HAVING ONE OF THE VALUES 1-10. IPMPAR(I) HAS THE VALUE ... C C INTEGERS. C C ASSUME INTEGERS ARE REPRESENTED IN THE N-DIGIT, BASE-A FORM C C SIGN ( X(N-1)*A**(N-1) + ... + X(1)*A + X(0) ) C C WHERE 0 .LE. X(I) .LT. A FOR I=0,...,N-1. C C IPMPAR(1) = A, THE BASE. C C IPMPAR(2) = N, THE NUMBER OF BASE-A DIGITS. C C IPMPAR(3) = A**N - 1, THE LARGEST MAGNITUDE. C C FLOATING-POINT NUMBERS. C C IT IS ASSUMED THAT THE SINGLE AND DOUBLE PRECISION FLOATING C POINT ARITHMETICS HAVE THE SAME BASE, SAY B, AND THAT THE C NONZERO NUMBERS ARE REPRESENTED IN THE FORM C C SIGN (B**E) * (X(1)/B + ... + X(M)/B**M) C C WHERE X(I) = 0,1,...,B-1 FOR I=1,...,M, C X(1) .GE. 1, AND EMIN .LE. E .LE. EMAX. C C IPMPAR(4) = B, THE BASE. C C SINGLE-PRECISION C C IPMPAR(5) = M, THE NUMBER OF BASE-B DIGITS. C C IPMPAR(6) = EMIN, THE SMALLEST EXPONENT E. C C IPMPAR(7) = EMAX, THE LARGEST EXPONENT E. C C DOUBLE-PRECISION C C IPMPAR(8) = M, THE NUMBER OF BASE-B DIGITS. C C IPMPAR(9) = EMIN, THE SMALLEST EXPONENT E. C C IPMPAR(10) = EMAX, THE LARGEST EXPONENT E. C C----------------------------------------------------------------------- C C TO DEFINE THIS FUNCTION FOR THE COMPUTER BEING USED, ACTIVATE C THE DATA STATMENTS FOR THE COMPUTER BY REMOVING THE C FROM C COLUMN 1. (ALL THE OTHER DATA STATEMENTS SHOULD HAVE C IN C COLUMN 1.) C C----------------------------------------------------------------------- C C IPMPAR IS AN ADAPTATION OF THE FUNCTION I1MACH, WRITTEN BY C P.A. FOX, A.D. HALL, AND N.L. SCHRYER (BELL LABORATORIES). C IPMPAR WAS FORMED BY A.H. MORRIS (NSWC). THE CONSTANTS ARE C FROM BELL LABORATORIES, NSWC, AND OTHER SOURCES. C C----------------------------------------------------------------------- C .. SCALAR ARGUMENTS .. INTEGER I C .. C .. LOCAL ARRAYS .. INTEGER IMACH(10) C .. C .. DATA STATEMENTS .. C C MACHINE CONSTANTS FOR AMDAHL MACHINES. C C DATA IMACH( 1) / 2 / C DATA IMACH( 2) / 31 / C DATA IMACH( 3) / 2147483647 / C DATA IMACH( 4) / 16 / C DATA IMACH( 5) / 6 / C DATA IMACH( 6) / -64 / C DATA IMACH( 7) / 63 / C DATA IMACH( 8) / 14 / C DATA IMACH( 9) / -64 / C DATA IMACH(10) / 63 / C C MACHINE CONSTANTS FOR THE AT&T 3B SERIES, AT&T C PC 7300, AND AT&T 6300. C C DATA IMACH( 1) / 2 / C DATA IMACH( 2) / 31 / C DATA IMACH( 3) / 2147483647 / C DATA IMACH( 4) / 2 / C DATA IMACH( 5) / 24 / C DATA IMACH( 6) / -125 / C DATA IMACH( 7) / 128 / C DATA IMACH( 8) / 53 / C DATA IMACH( 9) / -1021 / C DATA IMACH(10) / 1024 / C C MACHINE CONSTANTS FOR THE BURROUGHS 1700 SYSTEM. C C DATA IMACH( 1) / 2 / C DATA IMACH( 2) / 33 / C DATA IMACH( 3) / 8589934591 / C DATA IMACH( 4) / 2 / C DATA IMACH( 5) / 24 / C DATA IMACH( 6) / -256 / C DATA IMACH( 7) / 255 / C DATA IMACH( 8) / 60 / C DATA IMACH( 9) / -256 / C DATA IMACH(10) / 255 / C C MACHINE CONSTANTS FOR THE BURROUGHS 5700 SYSTEM. C C DATA IMACH( 1) / 2 / C DATA IMACH( 2) / 39 / C DATA IMACH( 3) / 549755813887 / C DATA IMACH( 4) / 8 / C DATA IMACH( 5) / 13 / C DATA IMACH( 6) / -50 / C DATA IMACH( 7) / 76 / C DATA IMACH( 8) / 26 / C DATA IMACH( 9) / -50 / C DATA IMACH(10) / 76 / C C MACHINE CONSTANTS FOR THE BURROUGHS 6700/7700 SYSTEMS. C C DATA IMACH( 1) / 2 / C DATA IMACH( 2) / 39 / C DATA IMACH( 3) / 549755813887 / C DATA IMACH( 4) / 8 / C DATA IMACH( 5) / 13 / C DATA IMACH( 6) / -50 / C DATA IMACH( 7) / 76 / C DATA IMACH( 8) / 26 / C DATA IMACH( 9) / -32754 / C DATA IMACH(10) / 32780 / C C MACHINE CONSTANTS FOR THE CDC 6000/7000 SERIES C 60 BIT ARITHMETIC, AND THE CDC CYBER 995 64 BIT C ARITHMETIC (NOS OPERATING SYSTEM). C C DATA IMACH( 1) / 2 / C DATA IMACH( 2) / 48 / C DATA IMACH( 3) / 281474976710655 / C DATA IMACH( 4) / 2 / C DATA IMACH( 5) / 48 / C DATA IMACH( 6) / -974 / C DATA IMACH( 7) / 1070 / C DATA IMACH( 8) / 95 / C DATA IMACH( 9) / -926 / C DATA IMACH(10) / 1070 / C C MACHINE CONSTANTS FOR THE CDC CYBER 995 64 BIT C ARITHMETIC (NOS/VE OPERATING SYSTEM). C C DATA IMACH( 1) / 2 / C DATA IMACH( 2) / 63 / C DATA IMACH( 3) / 9223372036854775807 / C DATA IMACH( 4) / 2 / C DATA IMACH( 5) / 48 / C DATA IMACH( 6) / -4096 / C DATA IMACH( 7) / 4095 / C DATA IMACH( 8) / 96 / C DATA IMACH( 9) / -4096 / C DATA IMACH(10) / 4095 / C C MACHINE CONSTANTS FOR THE CRAY 1, XMP, 2, AND 3. C C DATA IMACH( 1) / 2 / C DATA IMACH( 2) / 63 / C DATA IMACH( 3) / 9223372036854775807 / C DATA IMACH( 4) / 2 / C DATA IMACH( 5) / 47 / C DATA IMACH( 6) / -8189 / C DATA IMACH( 7) / 8190 / C DATA IMACH( 8) / 94 / C DATA IMACH( 9) / -8099 / C DATA IMACH(10) / 8190 / C C MACHINE CONSTANTS FOR THE DATA GENERAL ECLIPSE S/200. C C DATA IMACH( 1) / 2 / C DATA IMACH( 2) / 15 / C DATA IMACH( 3) / 32767 / C DATA IMACH( 4) / 16 / C DATA IMACH( 5) / 6 / C DATA IMACH( 6) / -64 / C DATA IMACH( 7) / 63 / C DATA IMACH( 8) / 14 / C DATA IMACH( 9) / -64 / C DATA IMACH(10) / 63 / C C MACHINE CONSTANTS FOR THE HARRIS 220. C C DATA IMACH( 1) / 2 / C DATA IMACH( 2) / 23 / C DATA IMACH( 3) / 8388607 / C DATA IMACH( 4) / 2 / C DATA IMACH( 5) / 23 / C DATA IMACH( 6) / -127 / C DATA IMACH( 7) / 127 / C DATA IMACH( 8) / 38 / C DATA IMACH( 9) / -127 / C DATA IMACH(10) / 127 / C C MACHINE CONSTANTS FOR THE HONEYWELL 600/6000 C AND DPS 8/70 SERIES. C C DATA IMACH( 1) / 2 / C DATA IMACH( 2) / 35 / C DATA IMACH( 3) / 34359738367 / C DATA IMACH( 4) / 2 / C DATA IMACH( 5) / 27 / C DATA IMACH( 6) / -127 / C DATA IMACH( 7) / 127 / C DATA IMACH( 8) / 63 / C DATA IMACH( 9) / -127 / C DATA IMACH(10) / 127 / C C MACHINE CONSTANTS FOR THE HP 2100 C 3 WORD DOUBLE PRECISION OPTION WITH FTN4 C C DATA IMACH( 1) / 2 / C DATA IMACH( 2) / 15 / C DATA IMACH( 3) / 32767 / C DATA IMACH( 4) / 2 / C DATA IMACH( 5) / 23 / C DATA IMACH( 6) / -128 / C DATA IMACH( 7) / 127 / C DATA IMACH( 8) / 39 / C DATA IMACH( 9) / -128 / C DATA IMACH(10) / 127 / C C MACHINE CONSTANTS FOR THE HP 2100 C 4 WORD DOUBLE PRECISION OPTION WITH FTN4 C C DATA IMACH( 1) / 2 / C DATA IMACH( 2) / 15 / C DATA IMACH( 3) / 32767 / C DATA IMACH( 4) / 2 / C DATA IMACH( 5) / 23 / C DATA IMACH( 6) / -128 / C DATA IMACH( 7) / 127 / C DATA IMACH( 8) / 55 / C DATA IMACH( 9) / -128 / C DATA IMACH(10) / 127 / C C MACHINE CONSTANTS FOR THE HP 9000. C C DATA IMACH( 1) / 2 / C DATA IMACH( 2) / 31 / C DATA IMACH( 3) / 2147483647 / C DATA IMACH( 4) / 2 / C DATA IMACH( 5) / 24 / C DATA IMACH( 6) / -126 / C DATA IMACH( 7) / 128 / C DATA IMACH( 8) / 53 / C DATA IMACH( 9) / -1021 / C DATA IMACH(10) / 1024 / C C MACHINE CONSTANTS FOR THE IBM 360/370 SERIES, C THE ICL 2900, THE ITEL AS/6, THE XEROX SIGMA C 5/7/9 AND THE SEL SYSTEMS 85/86. C C DATA IMACH( 1) / 2 / C DATA IMACH( 2) / 31 / C DATA IMACH( 3) / 2147483647 / C DATA IMACH( 4) / 16 / C DATA IMACH( 5) / 6 / C DATA IMACH( 6) / -64 / C DATA IMACH( 7) / 63 / C DATA IMACH( 8) / 14 / C DATA IMACH( 9) / -64 / C DATA IMACH(10) / 63 / C C MACHINE CONSTANTS FOR THE IBM PC. C C DATA IMACH(1)/2/ C DATA IMACH(2)/31/ C DATA IMACH(3)/2147483647/ C DATA IMACH(4)/2/ C DATA IMACH(5)/24/ C DATA IMACH(6)/-125/ C DATA IMACH(7)/128/ C DATA IMACH(8)/53/ C DATA IMACH(9)/-1021/ C DATA IMACH(10)/1024/ C C MACHINE CONSTANTS FOR THE MACINTOSH II - ABSOFT C MACFORTRAN II. C C DATA IMACH( 1) / 2 / C DATA IMACH( 2) / 31 / C DATA IMACH( 3) / 2147483647 / C DATA IMACH( 4) / 2 / C DATA IMACH( 5) / 24 / C DATA IMACH( 6) / -125 / C DATA IMACH( 7) / 128 / C DATA IMACH( 8) / 53 / C DATA IMACH( 9) / -1021 / C DATA IMACH(10) / 1024 / C C MACHINE CONSTANTS FOR THE MICROVAX - VMS FORTRAN. C C DATA IMACH( 1) / 2 / C DATA IMACH( 2) / 31 / C DATA IMACH( 3) / 2147483647 / C DATA IMACH( 4) / 2 / C DATA IMACH( 5) / 24 / C DATA IMACH( 6) / -127 / C DATA IMACH( 7) / 127 / C DATA IMACH( 8) / 56 / C DATA IMACH( 9) / -127 / C DATA IMACH(10) / 127 / C C MACHINE CONSTANTS FOR THE PDP-10 (KA PROCESSOR). C C DATA IMACH( 1) / 2 / C DATA IMACH( 2) / 35 / C DATA IMACH( 3) / 34359738367 / C DATA IMACH( 4) / 2 / C DATA IMACH( 5) / 27 / C DATA IMACH( 6) / -128 / C DATA IMACH( 7) / 127 / C DATA IMACH( 8) / 54 / C DATA IMACH( 9) / -101 / C DATA IMACH(10) / 127 / C C MACHINE CONSTANTS FOR THE PDP-10 (KI PROCESSOR). C C DATA IMACH( 1) / 2 / C DATA IMACH( 2) / 35 / C DATA IMACH( 3) / 34359738367 / C DATA IMACH( 4) / 2 / C DATA IMACH( 5) / 27 / C DATA IMACH( 6) / -128 / C DATA IMACH( 7) / 127 / C DATA IMACH( 8) / 62 / C DATA IMACH( 9) / -128 / C DATA IMACH(10) / 127 / C C MACHINE CONSTANTS FOR THE PDP-11 FORTRAN SUPPORTING C 32-BIT INTEGER ARITHMETIC. C C DATA IMACH( 1) / 2 / C DATA IMACH( 2) / 31 / C DATA IMACH( 3) / 2147483647 / C DATA IMACH( 4) / 2 / C DATA IMACH( 5) / 24 / C DATA IMACH( 6) / -127 / C DATA IMACH( 7) / 127 / C DATA IMACH( 8) / 56 / C DATA IMACH( 9) / -127 / C DATA IMACH(10) / 127 / C C MACHINE CONSTANTS FOR THE SEQUENT BALANCE 8000. C C DATA IMACH( 1) / 2 / C DATA IMACH( 2) / 31 / C DATA IMACH( 3) / 2147483647 / C DATA IMACH( 4) / 2 / C DATA IMACH( 5) / 24 / C DATA IMACH( 6) / -125 / C DATA IMACH( 7) / 128 / C DATA IMACH( 8) / 53 / C DATA IMACH( 9) / -1021 / C DATA IMACH(10) / 1024 / C C MACHINE CONSTANTS FOR THE SILICON GRAPHICS IRIS-4D C SERIES (MIPS R3000 PROCESSOR). C C DATA IMACH( 1) / 2 / C DATA IMACH( 2) / 31 / C DATA IMACH( 3) / 2147483647 / C DATA IMACH( 4) / 2 / C DATA IMACH( 5) / 24 / C DATA IMACH( 6) / -125 / C DATA IMACH( 7) / 128 / C DATA IMACH( 8) / 53 / C DATA IMACH( 9) / -1021 / C DATA IMACH(10) / 1024 / C C MACHINE CONSTANTS FOR IEEE ARITHMETIC MACHINES, SUCH AS THE AT&T C 3B SERIES, MOTOROLA 68000 BASED MACHINES (E.G. SUN 3 AND AT&T C PC 7300), AND 8087 BASED MICROS (E.G. IBM PC AND AT&T 6300). C DATA IMACH( 1) / 2 / DATA IMACH( 2) / 31 / DATA IMACH( 3) / 2147483647 / DATA IMACH( 4) / 2 / DATA IMACH( 5) / 24 / DATA IMACH( 6) / -125 / DATA IMACH( 7) / 128 / DATA IMACH( 8) / 53 / DATA IMACH( 9) / -1021 / DATA IMACH(10) / 1024 / C C MACHINE CONSTANTS FOR THE UNIVAC 1100 SERIES. C C DATA IMACH( 1) / 2 / C DATA IMACH( 2) / 35 / C DATA IMACH( 3) / 34359738367 / C DATA IMACH( 4) / 2 / C DATA IMACH( 5) / 27 / C DATA IMACH( 6) / -128 / C DATA IMACH( 7) / 127 / C DATA IMACH( 8) / 60 / C DATA IMACH( 9) /-1024 / C DATA IMACH(10) / 1023 / C C MACHINE CONSTANTS FOR THE VAX 11/780. C C DATA IMACH( 1) / 2 / C DATA IMACH( 2) / 31 / C DATA IMACH( 3) / 2147483647 / C DATA IMACH( 4) / 2 / C DATA IMACH( 5) / 24 / C DATA IMACH( 6) / -127 / C DATA IMACH( 7) / 127 / C DATA IMACH( 8) / 56 / C DATA IMACH( 9) / -127 / C DATA IMACH(10) / 127 / C IPMPAR = IMACH(I) RETURN END C------------------------------------------------------------- DOUBLE PRECISION FUNCTION LBICO(N,S) IMPLICIT DOUBLE PRECISION (A-H,O-P,R-Z),INTEGER (I-N),LOGICAL (Q) C********************************************************************** C C DOUBLE PRECISION FUNCTION LBICO(N,S) C CALCULATE LOG OF BINOMIAL COEFFICIENT C C C FUNCTION C C C CALCULATES THE LOG OF THE BINOMIAL COEFFICIENT OF S OBJECTS C IN A SET OF SIZE N. C C C ARGUMENTS C C C N --> SAMPLE SIZE. C N IS INTEGER C C S --> NUMBER OF OBJECTS IN SUBSET. C S IS INTEGER C C C METHOD C C C THE LOG OF THE BINOMIAL COEFFICIENT IS CALCULATED DIRECTLY, C USING ROUTINE XLFACT, WHICH CALCULATES THE FACTORIAL OF ITS C ONE INTEGER ARGUMENT USING A TABLE OF LOGARITHMS. C C********************************************************************** INTEGER N,S DOUBLE PRECISION NLFACT,SLFACT,NMSLFC DOUBLE PRECISION XLFACT C C *** MAIN BODY OF PROGRAM C C C IF S IS EQUAL TO ZERO OR N, RETURN ZERO. C IF S IS OUT OF VALID BOUNDS, STOP EXECUTION. C OTHERWISE, PROCEED WITH CALCULATION. C IF (.NOT. (S.EQ.0)) GO TO 10 LBICO = 0.0 GO TO 70 10 IF (.NOT. (S.EQ.N)) GO TO 20 LBICO = 0.0 GO TO 70 20 IF (.NOT. (S.LT.0)) GO TO 30 STOP 'S LESS THAN 0 IN LBICO' GO TO 70 30 IF (.NOT. (S.GT.N)) GO TO 40 STOP 'S GREATER THAN 0 IN LBICO' GO TO 70 40 IF (.NOT. (N.LE.0)) GO TO 50 STOP 'N LESS THAN OR EQUAL TO 0 IN LBICO' GO TO 70 C CALCULATE-LOG-BINOMIAL-COEFFICIENT 50 ASSIGN 60 TO I99993 GO TO 80 60 CONTINUE 70 RETURN C C C STOP '*** EXECUTION FLOWING INTO FLECS PROCEDURES ***' C TO CALCULATE-LOG-BINOMIAL-COEFFICIENT C********************************************************************** C C CALCULATE LOG OF BINOMIAL COEFFICIENT FROM LOGARITHMS OF C N!, S! AND (N-S)! C C********************************************************************** C C CALCULATE LOGARITHMS OF N!, S! AND (N-S)! * 80 NLFACT = XLFACT(N) SLFACT = XLFACT(S) NMSLFC = XLFACT(N-S) C C CALCULATE LOG OF BINOMIAL COEFFICIENT. C LBICO = NLFACT- (SLFACT+NMSLFC) C GO TO I99993 C END C--------------------------------------------------------------- DOUBLE PRECISION FUNCTION PSI(XX) C--------------------------------------------------------------------- C C EVALUATION OF THE DIGAMMA FUNCTION C C ----------- C C PSI(XX) IS ASSIGNED THE VALUE 0 WHEN THE DIGAMMA FUNCTION CANNOT C BE COMPUTED. C C THE MAIN COMPUTATION INVOLVES EVALUATION OF RATIONAL CHEBYSHEV C APPROXIMATIONS PUBLISHED IN MATH. COMP. 27, 123-127(1973) BY C CODY, STRECOK AND THACHER. C C--------------------------------------------------------------------- C PSI WAS WRITTEN AT ARGONNE NATIONAL LABORATORY FOR THE FUNPACK C PACKAGE OF SPECIAL FUNCTION SUBROUTINES. PSI WAS MODIFIED BY C A.H. MORRIS (NSWC). C--------------------------------------------------------------------- C .. SCALAR ARGUMENTS .. DOUBLE PRECISION XX C .. C .. LOCAL SCALARS .. DOUBLE PRECISION AUG,DEN,DX0,PIOV4,SGN,UPPER,W,X,XMAX1,XMX0, + XSMALL,Z INTEGER I,M,N,NQ C .. C .. LOCAL ARRAYS .. DOUBLE PRECISION P1(7),P2(4),Q1(6),Q2(4) C .. C .. EXTERNAL FUNCTIONS .. DOUBLE PRECISION SPMPAR INTEGER IPMPAR EXTERNAL SPMPAR,IPMPAR C .. C .. INTRINSIC FUNCTIONS .. INTRINSIC ABS,COS,DBLE,DLOG,DMIN1,INT,SIN C .. C .. DATA STATEMENTS .. C--------------------------------------------------------------------- C C PIOV4 = PI/4 C DX0 = ZERO OF PSI TO EXTENDED PRECISION C C--------------------------------------------------------------------- C--------------------------------------------------------------------- C C COEFFICIENTS FOR RATIONAL APPROXIMATION OF C PSI(X) / (X - X0), 0.5 .LE. X .LE. 3.0 C C--------------------------------------------------------------------- C--------------------------------------------------------------------- C C COEFFICIENTS FOR RATIONAL APPROXIMATION OF C PSI(X) - LN(X) + 1 / (2*X), X .GT. 3.0 C C--------------------------------------------------------------------- DATA PIOV4/.785398163397448D0/ DATA DX0/1.461632144968362341262659542325721325D0/ DATA P1(1)/.895385022981970D-02/,P1(2)/.477762828042627D+01/, + P1(3)/.142441585084029D+03/,P1(4)/.118645200713425D+04/, + P1(5)/.363351846806499D+04/,P1(6)/.413810161269013D+04/, + P1(7)/.130560269827897D+04/ DATA Q1(1)/.448452573429826D+02/,Q1(2)/.520752771467162D+03/, + Q1(3)/.221000799247830D+04/,Q1(4)/.364127349079381D+04/, + Q1(5)/.190831076596300D+04/,Q1(6)/.691091682714533D-05/ DATA P2(1)/-.212940445131011D+01/,P2(2)/-.701677227766759D+01/, + P2(3)/-.448616543918019D+01/,P2(4)/-.648157123766197D+00/ DATA Q2(1)/.322703493791143D+02/,Q2(2)/.892920700481861D+02/, + Q2(3)/.546117738103215D+02/,Q2(4)/.777788548522962D+01/ C .. C .. EXECUTABLE STATEMENTS .. C--------------------------------------------------------------------- C C MACHINE DEPENDENT CONSTANTS ... C C XMAX1 = THE SMALLEST POSITIVE FLOATING POINT CONSTANT C WITH ENTIRELY INTEGER REPRESENTATION. ALSO USED C AS NEGATIVE OF LOWER BOUND ON ACCEPTABLE NEGATIVE C ARGUMENTS AND AS THE POSITIVE ARGUMENT BEYOND WHICH C PSI MAY BE REPRESENTED AS ALOG(X). C C XSMALL = ABSOLUTE ARGUMENT BELOW WHICH PI*COTAN(PI*X) C MAY BE REPRESENTED BY 1/X. C C--------------------------------------------------------------------- XMAX1 = IPMPAR(3) XMAX1 = DMIN1(XMAX1,1.0D0/SPMPAR(1)) XSMALL = 1.D-9 C--------------------------------------------------------------------- X = XX AUG = 0.0D0 IF (X.GE.0.5D0) GO TO 50 C--------------------------------------------------------------------- C X .LT. 0.5, USE REFLECTION FORMULA C PSI(1-X) = PSI(X) + PI * COTAN(PI*X) C--------------------------------------------------------------------- IF (ABS(X).GT.XSMALL) GO TO 10 IF (X.EQ.0.0D0) GO TO 100 C--------------------------------------------------------------------- C 0 .LT. ABS(X) .LE. XSMALL. USE 1/X AS A SUBSTITUTE C FOR PI*COTAN(PI*X) C--------------------------------------------------------------------- AUG = -1.0D0/X GO TO 40 C--------------------------------------------------------------------- C REDUCTION OF ARGUMENT FOR COTAN C--------------------------------------------------------------------- 10 W = -X SGN = PIOV4 IF (W.GT.0.0D0) GO TO 20 W = -W SGN = -SGN C--------------------------------------------------------------------- C MAKE AN ERROR EXIT IF X .LE. -XMAX1 C--------------------------------------------------------------------- 20 IF (W.GE.XMAX1) GO TO 100 NQ = INT(W) W = W - DBLE(NQ) NQ = INT(W*4.0D0) W = 4.0D0* (W-DBLE(NQ)*.25D0) C--------------------------------------------------------------------- C W IS NOW RELATED TO THE FRACTIONAL PART OF 4.0 * X. C ADJUST ARGUMENT TO CORRESPOND TO VALUES IN FIRST C QUADRANT AND DETERMINE SIGN C--------------------------------------------------------------------- N = NQ/2 IF ((N+N).NE.NQ) W = 1.0D0 - W Z = PIOV4*W M = N/2 IF ((M+M).NE.N) SGN = -SGN C--------------------------------------------------------------------- C DETERMINE FINAL VALUE FOR -PI*COTAN(PI*X) C--------------------------------------------------------------------- N = (NQ+1)/2 M = N/2 M = M + M IF (M.NE.N) GO TO 30 C--------------------------------------------------------------------- C CHECK FOR SINGULARITY C--------------------------------------------------------------------- IF (Z.EQ.0.0D0) GO TO 100 C--------------------------------------------------------------------- C USE COS/SIN AS A SUBSTITUTE FOR COTAN, AND C SIN/COS AS A SUBSTITUTE FOR TAN C--------------------------------------------------------------------- AUG = SGN* ((COS(Z)/SIN(Z))*4.0D0) GO TO 40 30 AUG = SGN* ((SIN(Z)/COS(Z))*4.0D0) 40 X = 1.0D0 - X 50 IF (X.GT.3.0D0) GO TO 70 C--------------------------------------------------------------------- C 0.5 .LE. X .LE. 3.0 C--------------------------------------------------------------------- DEN = X UPPER = P1(1)*X C DO 60 I = 1,5 DEN = (DEN+Q1(I))*X UPPER = (UPPER+P1(I+1))*X 60 CONTINUE C DEN = (UPPER+P1(7))/ (DEN+Q1(6)) XMX0 = DBLE(X) - DX0 PSI = DEN*XMX0 + AUG RETURN C--------------------------------------------------------------------- C IF X .GE. XMAX1, PSI = LN(X) C--------------------------------------------------------------------- 70 IF (X.GE.XMAX1) GO TO 90 C--------------------------------------------------------------------- C 3.0 .LT. X .LT. XMAX1 C--------------------------------------------------------------------- W = 1.0D0/ (X*X) DEN = W UPPER = P2(1)*W C DO 80 I = 1,3 DEN = (DEN+Q2(I))*W UPPER = (UPPER+P2(I+1))*W 80 CONTINUE C AUG = UPPER/ (DEN+Q2(4)) - 0.5D0/X + AUG 90 PSI = AUG + DLOG(X) RETURN C--------------------------------------------------------------------- C ERROR RETURN C--------------------------------------------------------------------- 100 PSI = 0.0D0 RETURN END C--------------------------------------------------------------------- LOGICAL FUNCTION QRZERO(F,Y,XLO,XHI,ABSTOL,RELTOL) C********************************************************************** C C LOGICAL FUNCTION QRZERO(F,Y,XLO,XHI,ABSTOL,RELTOL) C QR ZERO OF A BOUNDED MONOTONE FUNCTION F - Y C C C FUNCTION C C C CONCISE DESCRIPTION - GIVEN A FUNCTION F C FIND XLO SUCH THAT F(XLO) = Y. C C MORE PRECISE DESCRIPTION - C C INPUT CONDITION. F IS A REAL FUNCTION OF A SINGLE C REAL ARGUMENT AND XLO AND XHI ARE SUCH THAT C (F(XLO)-Y)*(F(XHI)-Y) .LE. 0.0 C C IF THE INPUT CONDITION IS MET, QRZERO RETURNS .TRUE. C AND OUTPUT VALUES OF XLO AND XHI SATISFY THE FOLLOWING C (F(XLO)-Y)*(F(XHI)-Y) .LE. 0. C ABS(F(XLO)-Y) .LE. ABS(F(XHI)-Y) C ABS(XLO-XHI) .LE. TOL(X) C WHERE C TOL(X) = MAX(ABSTOL,RELTOL*ABS(X)) C C IF THIS ALGORITHM DOES NOT FIND XLO AND XHI SATISFYING C THESE CONDITIONS THEN QRZERO RETURNS .FALSE. THIS C IMPLIES THAT THE INPUT CONDITION WAS NOT MET. C C C ARGUMENTS C C C F --> A FUNCTION OF A SINGLE REAL ARGUMENT. C F MUST BE DECLARED EXTERNAL IN THE CALLING PROGRAM. C F IS EXTERNAL C C Y --> THE EQUATION TO BE SOLVED IS F(X) = Y. C Y IS REAL C C XLO <--> ON INPUT, THE LEFT ENDPOINT OF THE INTERVAL TO BE C SEARCHED FOR A SOLUTION. ON OUTPUT, AN ESTIMATE OF C AN IMPROVED INTERVAL (SEE FUNCTION). C XLO IS DOUBLE PRECISION C C XHI <--> ON INPUT, THE RIGHT ENDPOINT OF THE INTERVAL TO BE C FOR A SOLUTION. ON OUTPUT, AN ESTIMATE OF AN IMPROVED C INTERVAL. C XHI IS DOUBLE PRECISION C C ABSTOL, RELTOL --> TWO NUMBERS THAT DETERMINE THE ACCURACY C OF THE SOLUTION. SEE FUNCTION FOR A C PRECISE DEFINITION. C ABSTOL IS DOUBLE PRECISION C RELTOL IS DOUBLE PRECISION C C C METHOD C C C ALGORITHM R OF THE PAPER 'TWO EFFICIENT ALGORITHMS WITH C GUARANTEED CONVERGENCE FOR FINDING A ZERO OF A FUNCTION' C BY J. C. P. BUS AND T. J. DEKKER IN ACM TRANSACTIONS ON C MATHEMATICAL SOFTWARE, VOLUME 1, NO. 4 PAGE 330 C (DEC. '75) IS EMPLOYED TO FIND THE ZERO OF F(X)-Y. C C********************************************************************** C .. C .. SCALAR ARGUMENTS .. DOUBLE PRECISION Y,XLO,XHI,ABSTOL,RELTOL C .. C .. FUNCTION ARGUMENTS .. DOUBLE PRECISION F EXTERNAL F C .. C .. LOCAL SCALARS .. DOUBLE PRECISION A,B,C,D,FA,FB,FC,FD,FDA,FDB,M,MB,P,Q,TOL,W,ZX INTEGER EXT LOGICAL FIRST C .. C .. INTRINSIC FUNCTIONS .. INTRINSIC DABS,DMAX1,DSIGN C .. C .. STATEMENT FUNCTIONS .. DOUBLE PRECISION FTOL C .. C .. SAVE STATEMENT .. SAVE C .. C .. STATEMENT FUNCTION DEFINITIONS .. FTOL(ZX) = 0.5D0*DMAX1(ABSTOL,RELTOL*DABS(ZX)) C .. C .. EXECUTABLE STATEMENTS .. B = XLO FB = F(XLO) - Y XLO = XHI A = XLO FA = F(XLO) - Y FIRST = .TRUE. 10 C = A FC = FA EXT = 0 20 IF (.NOT. (DABS(FC).LT.DABS(FB))) GO TO 40 IF (.NOT. (C.NE.A)) GO TO 30 D = A FD = FA 30 A = B FA = FB XLO = C B = XLO FB = FC C = A FC = FA 40 TOL = FTOL(XLO) M = (C+B)*.5D0 MB = M - B IF (.NOT. (ABS(MB).GT.TOL)) GO TO 170 IF (.NOT. (EXT.GT.3)) GO TO 50 W = MB GO TO 130 50 TOL = DSIGN(TOL,MB) P = (B-A)*FB IF (.NOT. (FIRST)) GO TO 60 Q = FA - FB FIRST = .FALSE. GO TO 70 60 FDB = (FD-FB)/ (D-B) FDA = (FD-FA)/ (D-A) P = FDA*P Q = FDB*FA - FDA*FB 70 IF (.NOT. (P.LT.0.0D0)) GO TO 80 P = -P Q = -Q 80 IF (EXT.EQ.3) P = P*2.0D0 IF (.NOT. ((P*1.0D0).EQ.0.0D0.OR.P.LE. (Q*TOL))) GO TO 90 W = TOL GO TO 120 90 IF (.NOT. (P.LT. (MB*Q))) GO TO 100 W = P/Q GO TO 110 100 W = MB 110 CONTINUE 120 CONTINUE 130 D = A FD = FA A = B FA = FB B = B + W XLO = B FB = F(XLO) - Y IF (.NOT. ((FC*FB).GE.0.0D0)) GO TO 140 GO TO 10 140 IF (.NOT. (W.EQ.MB)) GO TO 150 EXT = 0 GO TO 160 150 EXT = EXT + 1 160 GO TO 20 170 XHI = C QRZERO = (FC.GE.0.0D0 .AND. FB.LE.0.0D0) .OR. + (FC.LT.0.0D0 .AND. FB.GE.0.0D0) RETURN END C---------------------------------------------------------------------------- LOGICAL FUNCTION QSTINV(F,QINCR,X,Y,SMALL,BIG,ABSSTP,RELSTP, + ABSTOL,RELTOL,QLEFT,QHI) C********************************************************************** C C QSTINV IS AN EXACT RENAMED COPY OF QMFINV. IT IS INTENDED TO BE C USED ONLY FOR INVERTING STATISTICAL FUNCTIONS AND NOT FOR GENERAL C PURPOSES. C C********************************************************************** C C .. SCALAR ARGUMENTS .. DOUBLE PRECISION ABSSTP,ABSTOL,BIG,RELSTP,RELTOL,SMALL,X,Y LOGICAL QHI,QINCR,QLEFT C .. C .. FUNCTION ARGUMENTS .. DOUBLE PRECISION F EXTERNAL F C .. C .. LOCAL SCALARS .. DOUBLE PRECISION STEP,XLB,XUB,YY,ZX,ZY,ZZ LOGICAL QBDD,QCOND,QLIM,QUP C .. C .. EXTERNAL FUNCTIONS .. LOGICAL QSTZRO EXTERNAL QSTZRO C .. C .. INTRINSIC FUNCTIONS .. INTRINSIC DABS,DMAX1,DMIN1 C .. C .. STATEMENT FUNCTIONS .. LOGICAL QXMON C .. C .. SAVE STATEMENT .. SAVE C .. C .. STATEMENT FUNCTION DEFINITIONS .. QXMON(ZX,ZY,ZZ) = ZX .LE. ZY .AND. ZY .LE. ZZ C .. C .. EXECUTABLE STATEMENTS .. QCOND = .NOT. QXMON(SMALL,X,BIG) IF (.NOT. (QCOND)) GO TO 10 STOP ' SMALL, X, BIG NOT MONOTONE IN QSTINV' 10 QSTINV = .TRUE. STEP = DMAX1(ABSSTP,RELSTP*DABS(X)) YY = F(X) - Y IF (YY.EQ.0.0D0) RETURN QUP = (QINCR .AND. (YY.LT.0.0D0)) .OR. + (.NOT.QINCR .AND. (YY.GT.0.0D0)) C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ C C HANDLE CASE IN WHICH WE MUST STEP HIGHER C C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ IF (.NOT. (QUP)) GO TO 70 XLB = X XUB = DMIN1(XLB+STEP,BIG) GO TO 30 20 IF (QCOND) GO TO 60 30 YY = F(XUB) - Y QBDD = (QINCR .AND. (YY.GE.0.0D0)) .OR. + (.NOT.QINCR .AND. (YY.LE.0.0D0)) QLIM = XUB .GE. BIG QCOND = QBDD .OR. QLIM IF (QCOND) GO TO 40 STEP = STEP*2.0D0 XLB = XUB XUB = DMIN1(XLB+STEP,BIG) 40 IF (.NOT. (QLIM.AND..NOT.QBDD)) GO TO 50 QSTINV = .FALSE. QLEFT = .FALSE. QHI = .NOT. QINCR X = BIG RETURN 50 GO TO 20 60 GO TO 130 C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ C C HANDLE CASE IN WHICH WE MUST STEP LOWER C C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 70 XUB = X XLB = DMAX1(XUB-STEP,SMALL) GO TO 90 80 IF (QCOND) GO TO 110 90 YY = F(XLB) - Y QBDD = (QINCR .AND. (YY.LE.0.0D0)) .OR. + (.NOT.QINCR .AND. (YY.GE.0.0D0)) QLIM = XLB .LE. SMALL QCOND = QBDD .OR. QLIM IF (QCOND) GO TO 100 STEP = STEP*2.0D0 XUB = XLB XLB = DMAX1(XUB-STEP,SMALL) 100 GO TO 80 110 IF (.NOT. (QLIM.AND..NOT.QBDD)) GO TO 120 QSTINV = .FALSE. QLEFT = .TRUE. QHI = QINCR X = SMALL RETURN 120 CONTINUE 130 IF (QSTZRO(F,Y,XLB,XUB,ABSTOL,RELTOL)) CONTINUE C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ C C IF WE REACH HERE, XLB AND XUB BOUND THE ZERO OF F. C C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ X = XLB RETURN END C------------------------------------------------------------ LOGICAL FUNCTION QSTZRO(F,Y,XLO,XHI,ABSTOL,RELTOL) C********************************************************************** C C C QSTZRO IS AN EXACT COPY OF QRZERO. QSTZRO IS INTENDED FOR USE ONL C FOR FINDING INVERSES OF CUMULATIVE STATISTICAL FUNCTIONS SO IS NOT C FURTHER DOCUMENTED FOR GENERAL USE. C C********************************************************************** C C .. SCALAR ARGUMENTS .. DOUBLE PRECISION ABSTOL,RELTOL,XHI,XLO,Y C .. C .. FUNCTION ARGUMENTS .. DOUBLE PRECISION F EXTERNAL F C .. C .. LOCAL SCALARS .. DOUBLE PRECISION A,B,C,D,FA,FB,FC,FD,FDA,FDB,M,MB,P,Q,TOL,W,ZX INTEGER EXT LOGICAL FIRST C .. C .. INTRINSIC FUNCTIONS .. INTRINSIC DABS,DMAX1,DSIGN C .. C .. STATEMENT FUNCTIONS .. DOUBLE PRECISION FTOL C .. C .. SAVE STATEMENT .. SAVE C .. C .. STATEMENT FUNCTION DEFINITIONS .. FTOL(ZX) = 0.5D0*DMAX1(ABSTOL,RELTOL*DABS(ZX)) C .. C .. EXECUTABLE STATEMENTS .. B = XLO FB = F(XLO) - Y XLO = XHI A = XLO FA = F(XLO) - Y FIRST = .TRUE. 10 C = A FC = FA EXT = 0 20 IF (.NOT. (ABS(FC).LT.ABS(FB))) GO TO 40 IF (.NOT. (C.NE.A)) GO TO 30 D = A FD = FA 30 A = B FA = FB XLO = C B = XLO FB = FC C = A FC = FA 40 TOL = FTOL(XLO) M = (C+B)*.5D0 MB = M - B IF (.NOT. (ABS(MB).GT.TOL)) GO TO 170 IF (.NOT. (EXT.GT.3)) GO TO 50 W = MB GO TO 130 50 TOL = SIGN(TOL,MB) P = (B-A)*FB IF (.NOT. (FIRST)) GO TO 60 Q = FA - FB FIRST = .FALSE. GO TO 70 60 FDB = (FD-FB)/ (D-B) FDA = (FD-FA)/ (D-A) P = FDA*P Q = FDB*FA - FDA*FB 70 IF (.NOT. (P.LT.0.0D0)) GO TO 80 P = -P Q = -Q 80 IF (EXT.EQ.3) P = P*2.0D0 IF (.NOT. ((P*1.0D0).EQ.0.0D0.OR.P.LE. (Q*TOL))) GO TO 90 W = TOL GO TO 120 90 IF (.NOT. (P.LT. (MB*Q))) GO TO 100 W = P/Q GO TO 110 100 W = MB 110 CONTINUE 120 CONTINUE 130 D = A FD = FA A = B FA = FB B = B + W XLO = B FB = F(XLO) - Y IF (.NOT. ((FC*FB).GE.0.0D0)) GO TO 140 GO TO 10 140 IF (.NOT. (W.EQ.MB)) GO TO 150 EXT = 0 GO TO 160 150 EXT = EXT + 1 160 GO TO 20 170 XHI = C QSTZRO = (FC.GE.0.0D0 .AND. FB.LE.0.0D0) .OR. + (FC.LT.0.0D0 .AND. FB.GE.0.0D0) RETURN END C----------------------------------------------------------- DOUBLE PRECISION FUNCTION RCOMP(A,X) C ------------------- C EVALUATION OF EXP(-X)*X**A/GAMMA(A) C ------------------- C RT2PIN = 1/SQRT(2*PI) C ------------------- C .. SCALAR ARGUMENTS .. DOUBLE PRECISION A,X C .. C .. LOCAL SCALARS .. DOUBLE PRECISION RT2PIN,T,T1,U C .. C .. EXTERNAL FUNCTIONS .. DOUBLE PRECISION GAM1,GAMMA,RLOG EXTERNAL GAM1,GAMMA,RLOG C .. C .. INTRINSIC FUNCTIONS .. INTRINSIC DLOG,EXP,SQRT C .. C .. DATA STATEMENTS .. DATA RT2PIN/.398942280401433D0/ C .. C .. EXECUTABLE STATEMENTS .. C ------------------- RCOMP = 0.0D0 IF (A.GE.20.0D0) GO TO 20 T = A*DLOG(X) - X IF (A.GE.1.0D0) GO TO 10 RCOMP = (A*EXP(T))* (1.0D0+GAM1(A)) RETURN 10 RCOMP = EXP(T)/GAMMA(A) RETURN C 20 U = X/A IF (U.EQ.0.0D0) RETURN T = (1.0D0/A)**2 T1 = (((0.75D0*T-1.0D0)*T+3.5D0)*T-105.0D0)/ (A*1260.0D0) T1 = T1 - A*RLOG(U) RCOMP = RT2PIN*SQRT(A)*EXP(T1) RETURN END C-------------------------------------------------------------------- DOUBLE PRECISION FUNCTION REVLPL(A,N,X) C C********************************************************************** C C DOUBLE PRECISION FUNCTION REVLPL(A,N,X) C DOUBLE PRECISION EVALUATE A POLYNOMIAL AT X C C C FUNCTION C C C RETURNS C A(1) + A(2)*X + ... + A(N)*X**(N-1) C C C ARGUMENTS C C C A --> ARRAY OF COEFFICIENTS OF THE POLYNOMIAL. C A IS DOUBLE PRECISION(N) C C N --> LENGTH OF A, ALSO DEGREE OF POLYNOMIAL - 1. C N IS INTEGER C C X --> POINT AT WHICH THE POLYNOMIAL IS TO BE EVALUATED. C X IS DOUBLE PRECISION C C********************************************************************** C C .. SCALAR ARGUMENTS .. DOUBLE PRECISION X INTEGER N C .. C .. ARRAY ARGUMENTS .. DOUBLE PRECISION A(N) C .. C .. LOCAL SCALARS .. DOUBLE PRECISION TERM INTEGER I C .. C .. EXECUTABLE STATEMENTS .. TERM = A(N) DO 10,I = N - 1,1,-1 TERM = A(I) + TERM*X 10 CONTINUE REVLPL = TERM RETURN END C------------------------------------------------------------------ DOUBLE PRECISION FUNCTION REXP(X) C----------------------------------------------------------------------- C EVALUATION OF THE FUNCTION EXP(X) - 1 C----------------------------------------------------------------------- C .. SCALAR ARGUMENTS .. DOUBLE PRECISION X C .. C .. LOCAL SCALARS .. DOUBLE PRECISION P1,P2,Q1,Q2,Q3,Q4,W C .. C .. INTRINSIC FUNCTIONS .. INTRINSIC ABS,EXP C .. C .. DATA STATEMENTS .. DATA P1/.914041914819518D-09/,P2/.238082361044469D-01/, + Q1/-.499999999085958D+00/,Q2/.107141568980644D+00/, + Q3/-.119041179760821D-01/,Q4/.595130811860248D-03/ C .. C .. EXECUTABLE STATEMENTS .. C----------------------- IF (ABS(X).GT.0.15D0) GO TO 10 REXP = X* (((P2*X+P1)*X+1.0D0)/ ((((Q4*X+Q3)*X+Q2)*X+Q1)*X+1.0D0)) RETURN C 10 W = EXP(X) IF (X.GT.0.0D0) GO TO 20 REXP = (W-0.5D0) - 0.5D0 RETURN 20 REXP = W* (0.5D0+ (0.5D0-1.0D0/W)) RETURN END C------------------------------------------------------------------------ DOUBLE PRECISION FUNCTION RLOG(X) C ------------------- C COMPUTATION OF X - 1 - LN(X) C ------------------- C .. SCALAR ARGUMENTS .. DOUBLE PRECISION X C .. C .. LOCAL SCALARS .. DOUBLE PRECISION A,B,P0,P1,P2,Q1,Q2,R,T,U,W,W1 C .. C .. INTRINSIC FUNCTIONS .. INTRINSIC DBLE,DLOG C .. C .. DATA STATEMENTS .. C ------------------- DATA A/.566749439387324D-01/ DATA B/.456512608815524D-01/ DATA P0/.333333333333333D+00/,P1/-.224696413112536D+00/, + P2/.620886815375787D-02/ DATA Q1/-.127408923933623D+01/,Q2/.354508718369557D+00/ C .. C .. EXECUTABLE STATEMENTS .. C ------------------- IF (X.LT.0.61D0 .OR. X.GT.1.57D0) GO TO 40 IF (X.LT.0.82D0) GO TO 10 IF (X.GT.1.18D0) GO TO 20 C C ARGUMENT REDUCTION C U = (X-0.5D0) - 0.5D0 W1 = 0.0D0 GO TO 30 C 10 U = DBLE(X) - 0.7D0 U = U/0.7D0 W1 = A - U*0.3D0 GO TO 30 C 20 U = 0.75D0*DBLE(X) - 1.D0 W1 = B + U/3.0D0 C C SERIES EXPANSION C 30 R = U/ (U+2.0D0) T = R*R W = ((P2*T+P1)*T+P0)/ ((Q2*T+Q1)*T+1.0D0) RLOG = 2.0D0*T* (1.0D0/ (1.0D0-R)-R*W) + W1 RETURN C C 40 R = (X-0.5D0) - 0.5D0 RLOG = R - DLOG(X) RETURN END C------------------------------------------------------- DOUBLE PRECISION FUNCTION RLOG1(X) C----------------------------------------------------------------------- C EVALUATION OF THE FUNCTION X - LN(1 + X) C----------------------------------------------------------------------- C .. SCALAR ARGUMENTS .. DOUBLE PRECISION X C .. C .. LOCAL SCALARS .. DOUBLE PRECISION A,B,H,P0,P1,P2,Q1,Q2,R,T,W,W1 C .. C .. INTRINSIC FUNCTIONS .. INTRINSIC DBLE,DLOG C .. C .. DATA STATEMENTS .. C------------------------ DATA A/.566749439387324D-01/ DATA B/.456512608815524D-01/ DATA P0/.333333333333333D+00/,P1/-.224696413112536D+00/, + P2/.620886815375787D-02/ DATA Q1/-.127408923933623D+01/,Q2/.354508718369557D+00/ C .. C .. EXECUTABLE STATEMENTS .. C------------------------ IF (X.LT.-0.39D0 .OR. X.GT.0.57D0) GO TO 40 IF (X.LT.-0.18D0) GO TO 10 IF (X.GT.0.18D0) GO TO 20 C C ARGUMENT REDUCTION C H = X W1 = 0.0D0 GO TO 30 C 10 H = DBLE(X) + 0.3D0 H = H/0.7D0 W1 = A - H*0.3D0 GO TO 30 C 20 H = 0.75D0*DBLE(X) - 0.25D0 W1 = B + H/3.0D0 C C SERIES EXPANSION C 30 R = H/ (H+2.0D0) T = R*R W = ((P2*T+P1)*T+P0)/ ((Q2*T+Q1)*T+1.0D0) RLOG1 = 2.0D0*T* (1.0D0/ (1.0D0-R)-R*W) + W1 RETURN C C 40 W = (X+0.5D0) + 0.5D0 RLOG1 = X - DLOG(W) RETURN END C--------------------------------------------------------- DOUBLE PRECISION FUNCTION SPMPAR(I) C----------------------------------------------------------------------- C C SPMPAR PROVIDES THE SINGLE PRECISION MACHINE CONSTANTS FOR C THE COMPUTER BEING USED. IT IS ASSUMED THAT THE ARGUMENT C I IS AN INTEGER HAVING ONE OF THE VALUES 1, 2, OR 3. IF THE C SINGLE PRECISION ARITHMETIC BEING USED HAS M BASE B DIGITS AND C ITS SMALLEST AND LARGEST EXPONENTS ARE EMIN AND EMAX, THEN C C SPMPAR(1) = B**(1 - M), THE MACHINE PRECISION, C C SPMPAR(2) = B**(EMIN - 1), THE SMALLEST MAGNITUDE, C C SPMPAR(3) = B**EMAX*(1 - B**(-M)), THE LARGEST MAGNITUDE. C C----------------------------------------------------------------------- C WRITTEN BY C ALFRED H. MORRIS, JR. C NAVAL SURFACE WARFARE CENTER C DAHLGREN VIRGINIA C----------------------------------------------------------------------- C .. SCALAR ARGUMENTS .. INTEGER I C .. C .. LOCAL SCALARS .. DOUBLE PRECISION B,BINV,BM1,ONE,W,Z INTEGER EMAX,EMIN,IBETA,M C .. C .. EXTERNAL FUNCTIONS .. INTEGER IPMPAR EXTERNAL IPMPAR C .. C .. INTRINSIC FUNCTIONS .. INTRINSIC DBLE C .. C .. EXECUTABLE STATEMENTS .. C IF (I.GT.1) GO TO 10 B = IPMPAR(4) M = IPMPAR(8) SPMPAR = B** (1-M) RETURN C 10 IF (I.GT.2) GO TO 20 B = IPMPAR(4) EMIN = IPMPAR(9) ONE = DBLE(1) BINV = ONE/B W = B** (EMIN+2) SPMPAR = ((W*BINV)*BINV)*BINV RETURN C 20 IBETA = IPMPAR(4) M = IPMPAR(8) EMAX = IPMPAR(10) C B = IBETA BM1 = IBETA - 1 ONE = DBLE(1) Z = B** (M-1) W = ((Z-ONE)*B+BM1)/ (B*Z) C Z = B** (EMAX-2) SPMPAR = ((W*Z)*B)*B RETURN END C------------------------------------------------------- DOUBLE PRECISION FUNCTION T1(P,DF) C C C .. SCALAR ARGUMENTS .. DOUBLE PRECISION DF,P C .. C .. LOCAL SCALARS .. DOUBLE PRECISION DENPOW,SUM,TERM,X,XP,XX INTEGER I C .. C .. LOCAL ARRAYS .. DOUBLE PRECISION COEF(5,4),DENOM(4) INTEGER IDEG(4) C .. C .. EXTERNAL FUNCTIONS .. DOUBLE PRECISION INVNOR,REVLPL EXTERNAL INVNOR,REVLPL C .. C .. INTRINSIC FUNCTIONS .. INTRINSIC DABS C .. C .. DATA STATEMENTS .. DATA (COEF(I,1),I=1,5)/1.0D0,1.0D0,3*0.0D0/ DATA (COEF(I,2),I=1,5)/3.0D0,16.0D0,5.0D0,2*0.0D0/ DATA (COEF(I,3),I=1,5)/-15.0D0,17.0D0,19.0D0,3.0D0,0.0D0/ DATA (COEF(I,4),I=1,5)/-945.0D0,-1920.0D0,1482.0D0,776.0D0,79.0D0/ DATA IDEG/2,3,4,5/ DATA DENOM/4.0D0,96.0D0,384.0D0,92160.0D0/ C .. C .. EXECUTABLE STATEMENTS .. X = ABS(INVNOR(P)) XX = X*X SUM = X DENPOW = 1.0D0 DO 10,I = 1,4 TERM = REVLPL(COEF(1,I),IDEG(I),XX)*X DENPOW = DENPOW*DF SUM = SUM + TERM/ (DENPOW*DENOM(I)) 10 CONTINUE IF (.NOT. (P.GE.0.5D0)) GO TO 20 XP = SUM GO TO 30 20 XP = -SUM 30 T1 = XP RETURN END C----------------------------------------------------------- DOUBLE PRECISION FUNCTION XLFACT(N) C********************************************************************** C C DOUBLE PRECISION FUNCTION XLFACT(N) C X LOG FACTORIAL C RETURNS THE LOG OF N! C C C ARGUMENTS C C C N --> INTEGER VALUE WHOSE FACTORIAL IS BEING EVALUATED C********************************************************************** C .. SCALAR ARGUMENTS .. INTEGER N C .. C .. EXTERNAL FUNCTIONS .. DOUBLE PRECISION DLNGAM EXTERNAL DLNGAM C .. C .. INTRINSIC FUNCTIONS .. INTRINSIC DBLE C .. C .. EXECUTABLE STATEMENTS .. XLFACT = DLNGAM(DBLE(N+1)) RETURN END C--------------------------------------------------------------------- SUBROUTINE ZROR(STATUS,X,FX,XLO,XHI,QLEFT,QHI) C********************************************************************** C C SUBROUTINE ZROR(STATUS, X, FX, XLO, XHI, QLEFT, QHI) C ZERO OF A FUNCTION -- REVERSE COMMUNICATION C C C FUNCTION C C C PERFORMS THE ZERO FINDING. STZROR MUST HAVE BEEN CALLED BEFORE C THIS ROUTINE IN ORDER TO SET ITS PARAMETERS. C C C ARGUMENTS C C C STATUS <--> AT THE BEGINNING OF A ZERO FINDING PROBLEM, STATUS C SHOULD BE SET TO 0 AND ZROR INVOKED. (THE VALUE C OF OTHER PARAMETERS WILL BE IGNORED ON THIS CALL.) C C WHEN ZROR NEEDS THE FUNCTION EVALUATED, IT WILL SET C STATUS TO 1 AND RETURN. THE VALUE OF THE FUNCTION C SHOULD BE SET IN FX AND ZROR AGAIN CALLED WITHOUT C CHANGING ANY OF ITS OTHER PARAMETERS. C C WHEN ZROR HAS FINISHED WITHOUT ERROR, IT WILL RETURN C WITH STATUS 0. IN THAT CASE (XLO,XHI) BOUND THE ANSWE C C IF ZROR FINDS AN ERROR (WHICH IMPLIES THAT F(XLO)-Y AN C F(XHI)-Y HAVE THE SAME SIGN, IT RETURNS STATUS -1. IN C THIS CASE, XLO AND XHI ARE UNDEFINED. C INTEGER STATUS C C X <-- THE VALUE OF X AT WHICH F(X) IS TO BE EVALUATED. C DOUBLE PRECISIONX C C FX --> THE VALUE OF F(X) CALCULATED WHEN ZROR RETURNS WITH C STATUS = 1. C DOUBLE PRECISIONFX C C XLO <-- WHEN ZROR RETURNS WITH STATUS = 0, XLO BOUNDS THE C INVERVAL IN X CONTAINING THE SOLUTION BELOW. C DOUBLE PRECISIONXLO C C XHI <-- WHEN ZROR RETURNS WITH STATUS = 0, XHI BOUNDS THE C INVERVAL IN X CONTAINING THE SOLUTION ABOVE. C DOUBLE PRECISIONXHI C C QLEFT <-- .TRUE. IF THE STEPPING SEARCH TERMINATED UNSUCESSFULLY C AT XLO. IF IT IS .FALSE. THE SEARCH TERMINATED C UNSUCESSFULLY AT XHI. C QLEFT IS LOGICAL C C QHI <-- .TRUE. IF F(X) .GT. Y AT THE TERMINATION OF THE C SEARCH AND .FALSE. IF F(X) .LT. Y AT THE C TERMINATION OF THE SEARCH. C QHI IS LOGICAL C C********************************************************************** C .. SCALAR ARGUMENTS .. DOUBLE PRECISION FX,X,XHI,XLO,ZABSTL,ZRELTL,ZXHI,ZXLO INTEGER STATUS LOGICAL QHI,QLEFT C .. C .. SAVE STATEMENT .. SAVE C .. C .. LOCAL SCALARS .. DOUBLE PRECISION A,ABSTOL,B,C,D,FA,FB,FC,FD,FDA,FDB,M,MB,P,Q, +RELTOL,TOL,W,XXHI,XXLO,ZX INTEGER EXT,I99999 LOGICAL FIRST,QRZERO C .. C .. INTRINSIC FUNCTIONS .. INTRINSIC DABS,DMAX1,DSIGN C .. C .. STATEMENT FUNCTIONS .. DOUBLE PRECISION FTOL C .. C .. STATEMENT FUNCTION DEFINITIONS .. FTOL(ZX) = 0.5D0*DMAX1(ABSTOL,RELTOL*DABS(ZX)) C .. C .. EXECUTABLE STATEMENTS .. IF (STATUS.GT.0) GO TO 280 XLO = XXLO XHI = XXHI B = XLO X = XLO C GET-FUNCTION-VALUE ASSIGN 10 TO I99999 GO TO 270 10 FB = FX XLO = XHI A = XLO X = XLO C GET-FUNCTION-VALUE ASSIGN 20 TO I99999 GO TO 270 C C CHECK THAT F(ZXLO) < 0 < F(ZXHI) OR C F(ZXLO) > 0 > F(ZXHI) C 20 IF (.NOT. (FB.LT.0.0D0)) GO TO 40 IF (.NOT. (FX.LT.0.0D0)) GO TO 30 STATUS = -1 QLEFT = FX .LT. FB QHI = .FALSE. RETURN 30 CONTINUE 40 IF (.NOT. (FB.GT.0.0D0)) GO TO 60 IF (.NOT. (FX.GT.0.0D0)) GO TO 50 STATUS = -1 QLEFT = FX .GT. FB QHI = .TRUE. RETURN 50 CONTINUE 60 FA = FX C FIRST = .TRUE. 70 C = A FC = FA EXT = 0 80 IF (.NOT. (DABS(FC).LT.DABS(FB))) GO TO 100 IF (.NOT. (C.NE.A)) GO TO 90 D = A FD = FA 90 A = B FA = FB XLO = C B = XLO FB = FC C = A FC = FA 100 TOL = FTOL(XLO) M = (C+B)*.5 MB = M - B IF (.NOT. (DABS(MB).GT.TOL)) GO TO 240 IF (.NOT. (EXT.GT.3)) GO TO 110 W = MB GO TO 190 110 TOL = DSIGN(TOL,MB) P = (B-A)*FB IF (.NOT. (FIRST)) GO TO 120 Q = FA - FB FIRST = .FALSE. GO TO 130 120 FDB = (FD-FB)/ (D-B) FDA = (FD-FA)/ (D-A) P = FDA*P Q = FDB*FA - FDA*FB 130 IF (.NOT. (P.LT.0.0D0)) GO TO 140 P = -P Q = -Q 140 IF (EXT.EQ.3) P = P*2.0D0 IF (.NOT. ((P*1.0).EQ.0.0D0.OR.P.LE. (Q*TOL))) GO TO 150 W = TOL GO TO 180 150 IF (.NOT. (P.LT. (MB*Q))) GO TO 160 W = P/Q GO TO 170 160 W = MB 170 CONTINUE 180 CONTINUE 190 D = A FD = FA A = B FA = FB B = B + W XLO = B X = XLO C GET-FUNCTION-VALUE ASSIGN 200 TO I99999 GO TO 270 200 FB = FX IF (.NOT. ((FC*FB).GE.0.0D0)) GO TO 210 GO TO 70 210 IF (.NOT. (W.EQ.MB)) GO TO 220 EXT = 0 GO TO 230 220 EXT = EXT + 1 230 GO TO 80 240 XHI = C QRZERO = (FC.GE.0.0D0 .AND. FB.LE.0.0D0) .OR. + (FC.LT.0.0D0 .AND. FB.GE.0.0D0) IF (.NOT. (QRZERO)) GO TO 250 STATUS = 0 GO TO 260 250 STATUS = -1 260 RETURN ENTRY STZROR(ZXLO,ZXHI,ZABSTL,ZRELTL) C********************************************************************** C C SUBROUTINE STZROR( XLO, XHI, ABSTOL, RELTOL ) C SET ZERO FINDER - REVERSE COMMUNICATION VERSION C C C FUNCTION C C C C SETS QUANTITIES NEEDED BY ZROR. THE FUNCTION OF ZROR C AND THE QUANTITIES SET IS GIVEN HERE. C C CONCISE DESCRIPTION - GIVEN A FUNCTION F C FIND XLO SUCH THAT F(XLO) = 0. C C MORE PRECISE DESCRIPTION - C C INPUT CONDITION. F IS A REAL FUNCTION OF A SINGLE C REAL ARGUMENT AND XLO AND XHI ARE SUCH THAT C F(XLO)*F(XHI) .LE. 0.0 C C IF THE INPUT CONDITION IS MET, QRZERO RETURNS .TRUE. C AND OUTPUT VALUES OF XLO AND XHI SATISFY THE FOLLOWING C F(XLO)*F(XHI) .LE. 0. C ABS(F(XLO) .LE. ABS(F(XHI) C ABS(XLO-XHI) .LE. TOL(X) C WHERE C TOL(X) = DMAX1(ABSTOL,RELTOL*ABS(X)) C C IF THIS ALGORITHM DOES NOT FIND XLO AND XHI SATISFYING C THESE CONDITIONS THEN QRZERO RETURNS .FALSE. THIS C IMPLIES THAT THE INPUT CONDITION WAS NOT MET. C C C ARGUMENTS C C C XLO --> THE LEFT ENDPOINT OF THE INTERVAL TO BE C SEARCHED FOR A SOLUTION. C XLO IS DOUBLE PRECISION C C XHI --> THE RIGHT ENDPOINT OF THE INTERVAL TO BE C FOR A SOLUTION. C XHI IS DOUBLE PRECISION C C ABSTOL, RELTOL --> TWO NUMBERS THAT DETERMINE THE ACCURACY C OF THE SOLUTION. SEE FUNCTION FOR A C PRECISE DEFINITION. C ABSTOL IS DOUBLE PRECISION C RELTOL IS DOUBLE PRECISION C C C METHOD C C C ALGORITHM R OF THE PAPER 'TWO EFFICIENT ALGORITHMS WITH C GUARANTEED CONVERGENCE FOR FINDING A ZERO OF A FUNCTION' C BY J. C. P. BUS AND T. J. DEKKER IN ACM TRANSACTIONS ON C MATHEMATICAL SOFTWARE, VOLUME 1, NO. 4 PAGE 330 C (DEC. '75) IS EMPLOYED TO FIND THE ZERO OF F(X)-Y. C C********************************************************************** XXLO = ZXLO XXHI = ZXHI ABSTOL = ZABSTL RELTOL = ZRELTL RETURN STOP '*** EXECUTION FLOWING INTO FLECS PROCEDURES ***' C TO GET-FUNCTION-VALUE 270 STATUS = 1 RETURN 280 CONTINUE GO TO I99999 END