MODULE ERFCOREMOD IMPLICIT NONE INTERFACE CALERF MODULE PROCEDURE CALERF END INTERFACE INTERFACE DERF MODULE PROCEDURE DERF END INTERFACE INTERFACE DERFC MODULE PROCEDURE DERFC END INTERFACE INTERFACE DERFCX MODULE PROCEDURE DERFCX END INTERFACE CONTAINS C-------------------------------------------------------------------- C C DERF subprogram computes approximate values for erf(x). C (see comments heading CALERF). C C Author/date: W. J. Cody, January 8, 1985 C C-------------------------------------------------------------------- FUNCTION DERF( X ) RESULT (VALUE) IMPLICIT NONE DOUBLE PRECISION, INTENT(IN) :: X DOUBLE PRECISION :: VALUE INTEGER, PARAMETER :: JINT = 0 CALL CALERF(X,VALUE,JINT) RETURN END FUNCTION DERF C-------------------------------------------------------------------- C C DERFC subprogram computes approximate values for erfc(x). C (see comments heading CALERF). C C Author/date: W. J. Cody, January 8, 1985 C C-------------------------------------------------------------------- FUNCTION DERFC( X ) RESULT (VALUE) IMPLICIT NONE DOUBLE PRECISION, INTENT(IN) :: X DOUBLE PRECISION :: VALUE INTEGER, PARAMETER :: JINT = 1 CALL CALERF(X,VALUE,JINT) RETURN END FUNCTION DERFC C------------------------------------------------------------------ C C DERFCX subprogram computes approximate values for exp(x*x) * erfc(x). C (see comments heading CALERF). C C Author/date: W. J. Cody, March 30, 1987 C C------------------------------------------------------------------ FUNCTION DERFCX( X ) RESULT (VALUE) IMPLICIT NONE DOUBLE PRECISION, INTENT(IN) :: X DOUBLE PRECISION :: VALUE INTEGER, PARAMETER :: JINT = 2 CALL CALERF(X,VALUE,JINT) RETURN END FUNCTION DERFCX SUBROUTINE CALERF(ARG,RESULT,JINT) IMPLICIT NONE C------------------------------------------------------------------ C C CALERF packet evaluates erf(x), erfc(x), and exp(x*x)*erfc(x) C for a real argument x. It contains three FUNCTION type C subprograms: ERF, ERFC, and ERFCX (or DERF, DERFC, and DERFCX), C and one SUBROUTINE type subprogram, CALERF. The calling C statements for the primary entries are: C C Y=ERF(X) (or Y=DERF(X)), C C Y=ERFC(X) (or Y=DERFC(X)), C and C Y=ERFCX(X) (or Y=DERFCX(X)). C C The routine CALERF is intended for internal packet use only, C all computations within the packet being concentrated in this C routine. The function subprograms invoke CALERF with the C statement C C CALL CALERF(ARG,RESULT,JINT) C C where the parameter usage is as follows C C Function Parameters for CALERF C call ARG Result JINT C C ERF(ARG) ANY REAL ARGUMENT ERF(ARG) 0 C ERFC(ARG) ABS(ARG) .LT. XBIG ERFC(ARG) 1 C ERFCX(ARG) XNEG .LT. ARG .LT. XMAX ERFCX(ARG) 2 C C The main computation evaluates near-minimax approximations C from "Rational Chebyshev approximations for the error function" C by W. J. Cody, Math. Comp., 1969, PP. 631-638. This C transportable program uses rational functions that theoretically C approximate erf(x) and erfc(x) to at least 18 significant C decimal digits. The accuracy achieved depends on the arithmetic C system, the compiler, the intrinsic functions, and proper C selection of the machine-dependent constants. C C******************************************************************* C******************************************************************* C C Explanation of machine-dependent constants C C XMIN = the smallest positive floating-point number. C XINF = the largest positive finite floating-point number. C XNEG = the largest negative argument acceptable to ERFCX; C the negative of the solution to the equation C 2*exp(x*x) = XINF. C XSMALL = argument below which erf(x) may be represented by C 2*x/sqrt(pi) and above which x*x will not underflow. C A conservative value is the largest machine number X C such that 1.0 + X = 1.0 to machine precision. C XBIG = largest argument acceptable to ERFC; solution to C the equation: W(x) * (1-0.5/x**2) = XMIN, where C W(x) = exp(-x*x)/[x*sqrt(pi)]. C XHUGE = argument above which 1.0 - 1/(2*x*x) = 1.0 to C machine precision. A conservative value is C 1/[2*sqrt(XSMALL)] C XMAX = largest acceptable argument to ERFCX; the minimum C of XINF and 1/[sqrt(pi)*XMIN]. C C Approximate values for some important machines are: C C XMIN XINF XNEG XSMALL C C C 7600 (S.P.) 3.13E-294 1.26E+322 -27.220 7.11E-15 C CRAY-1 (S.P.) 4.58E-2467 5.45E+2465 -75.345 7.11E-15 C IEEE (IBM/XT, C SUN, etc.) (S.P.) 1.18E-38 3.40E+38 -9.382 5.96E-8 C IEEE (IBM/XT, C SUN, etc.) (D.P.) 2.23D-308 1.79D+308 -26.628 1.11D-16 C IBM 195 (D.P.) 5.40D-79 7.23E+75 -13.190 1.39D-17 C UNIVAC 1108 (D.P.) 2.78D-309 8.98D+307 -26.615 1.73D-18 C VAX D-Format (D.P.) 2.94D-39 1.70D+38 -9.345 1.39D-17 C VAX G-Format (D.P.) 5.56D-309 8.98D+307 -26.615 1.11D-16 C C C XBIG XHUGE XMAX C C C 7600 (S.P.) 25.922 8.39E+6 1.80X+293 C CRAY-1 (S.P.) 75.326 8.39E+6 5.45E+2465 C IEEE (IBM/XT, C SUN, etc.) (S.P.) 9.194 2.90E+3 4.79E+37 C IEEE (IBM/XT, C SUN, etc.) (D.P.) 26.543 6.71D+7 2.53D+307 C IBM 195 (D.P.) 13.306 1.90D+8 7.23E+75 C UNIVAC 1108 (D.P.) 26.582 5.37D+8 8.98D+307 C VAX D-Format (D.P.) 9.269 1.90D+8 1.70D+38 C VAX G-Format (D.P.) 26.569 6.71D+7 8.98D+307 C C******************************************************************* C******************************************************************* C C Error returns C C The program returns ERFC = 0 for ARG .GE. XBIG; C C ERFCX = XINF for ARG .LT. XNEG; C and C ERFCX = 0 for ARG .GE. XMAX. C C C Intrinsic functions required are: C C ABS, AINT, EXP C C C Author: W. J. Cody C Mathematics and Computer Science Division C Argonne National Laboratory C Argonne, IL 60439 C C Latest modification: March 19, 1990 C Updated to F90 by pab 23.03.2003 C C------------------------------------------------------------------ DOUBLE PRECISION, INTENT(IN) :: ARG INTEGER, INTENT(IN) :: JINT DOUBLE PRECISION, INTENT(INOUT):: RESULT ! Local variables INTEGER :: I DOUBLE PRECISION :: DEL,X,XDEN,XNUM,Y,YSQ C------------------------------------------------------------------ C Mathematical constants C------------------------------------------------------------------ DOUBLE PRECISION, PARAMETER :: ZERO = 0.0D0 DOUBLE PRECISION, PARAMETER :: HALF = 0.05D0 DOUBLE PRECISION, PARAMETER :: ONE = 1.0D0 DOUBLE PRECISION, PARAMETER :: TWO = 2.0D0 DOUBLE PRECISION, PARAMETER :: FOUR = 4.0D0 DOUBLE PRECISION, PARAMETER :: SIXTEN = 16.0D0 DOUBLE PRECISION, PARAMETER :: SQRPI = 5.6418958354775628695D-1 DOUBLE PRECISION, PARAMETER :: THRESH = 0.46875D0 C------------------------------------------------------------------ C Machine-dependent constants C------------------------------------------------------------------ DOUBLE PRECISION, PARAMETER :: XNEG = -26.628D0 DOUBLE PRECISION, PARAMETER :: XSMALL = 1.11D-16 DOUBLE PRECISION, PARAMETER :: XBIG = 26.543D0 DOUBLE PRECISION, PARAMETER :: XHUGE = 6.71D7 DOUBLE PRECISION, PARAMETER :: XMAX = 2.53D307 DOUBLE PRECISION, PARAMETER :: XINF = 1.79D308 !--------------------------------------------------------------- ! Coefficents to the rational polynomials !-------------------------------------------------------------- DOUBLE PRECISION, DIMENSION(5) :: A, Q DOUBLE PRECISION, DIMENSION(4) :: B DOUBLE PRECISION, DIMENSION(9) :: C DOUBLE PRECISION, DIMENSION(8) :: D DOUBLE PRECISION, DIMENSION(6) :: P C------------------------------------------------------------------ C Coefficients for approximation to erf in first interval C------------------------------------------------------------------ PARAMETER (A = (/ 3.16112374387056560D00, & 1.13864154151050156D02,3.77485237685302021D02, & 3.20937758913846947D03, 1.85777706184603153D-1/)) PARAMETER ( B = (/2.36012909523441209D01,2.44024637934444173D02, & 1.28261652607737228D03,2.84423683343917062D03/)) C------------------------------------------------------------------ C Coefficients for approximation to erfc in second interval C------------------------------------------------------------------ PARAMETER ( C=(/5.64188496988670089D-1,8.88314979438837594D0, 1 6.61191906371416295D01,2.98635138197400131D02, 2 8.81952221241769090D02,1.71204761263407058D03, 3 2.05107837782607147D03,1.23033935479799725D03, 4 2.15311535474403846D-8/)) PARAMETER ( D =(/1.57449261107098347D01,1.17693950891312499D02, 1 5.37181101862009858D02,1.62138957456669019D03, 2 3.29079923573345963D03,4.36261909014324716D03, 3 3.43936767414372164D03,1.23033935480374942D03/)) C------------------------------------------------------------------ C Coefficients for approximation to erfc in third interval C------------------------------------------------------------------ PARAMETER ( P =(/3.05326634961232344D-1,3.60344899949804439D-1, 1 1.25781726111229246D-1,1.60837851487422766D-2, 2 6.58749161529837803D-4,1.63153871373020978D-2/)) PARAMETER (Q =(/2.56852019228982242D00,1.87295284992346047D00, 1 5.27905102951428412D-1,6.05183413124413191D-2, 2 2.33520497626869185D-3/)) C------------------------------------------------------------------ X = ARG Y = ABS(X) IF (Y .LE. THRESH) THEN C------------------------------------------------------------------ C Evaluate erf for |X| <= 0.46875 C------------------------------------------------------------------ !YSQ = ZERO IF (Y .GT. XSMALL) THEN YSQ = Y * Y XNUM = A(5)*YSQ XDEN = YSQ DO I = 1, 3 XNUM = (XNUM + A(I)) * YSQ XDEN = (XDEN + B(I)) * YSQ END DO RESULT = X * (XNUM + A(4)) / (XDEN + B(4)) ELSE RESULT = X * A(4) / B(4) ENDIF IF (JINT .NE. 0) RESULT = ONE - RESULT IF (JINT .EQ. 2) RESULT = EXP(YSQ) * RESULT GO TO 800 C------------------------------------------------------------------ C Evaluate erfc for 0.46875 <= |X| <= 4.0 C------------------------------------------------------------------ ELSE IF (Y .LE. FOUR) THEN XNUM = C(9)*Y XDEN = Y DO I = 1, 7 XNUM = (XNUM + C(I)) * Y XDEN = (XDEN + D(I)) * Y END DO RESULT = (XNUM + C(8)) / (XDEN + D(8)) IF (JINT .NE. 2) THEN YSQ = AINT(Y*SIXTEN)/SIXTEN DEL = (Y-YSQ)*(Y+YSQ) RESULT = EXP(-YSQ*YSQ) * EXP(-DEL) * RESULT END IF C------------------------------------------------------------------ C Evaluate erfc for |X| > 4.0 C------------------------------------------------------------------ ELSE RESULT = ZERO IF (Y .GE. XBIG) THEN IF ((JINT .NE. 2) .OR. (Y .GE. XMAX)) GO TO 300 IF (Y .GE. XHUGE) THEN RESULT = SQRPI / Y GO TO 300 END IF END IF YSQ = ONE / (Y * Y) XNUM = P(6)*YSQ XDEN = YSQ DO I = 1, 4 XNUM = (XNUM + P(I)) * YSQ XDEN = (XDEN + Q(I)) * YSQ ENDDO RESULT = YSQ *(XNUM + P(5)) / (XDEN + Q(5)) RESULT = (SQRPI - RESULT) / Y IF (JINT .NE. 2) THEN YSQ = AINT(Y*SIXTEN)/SIXTEN DEL = (Y-YSQ)*(Y+YSQ) RESULT = EXP(-YSQ*YSQ) * EXP(-DEL) * RESULT END IF END IF C------------------------------------------------------------------ C Fix up for negative argument, erf, etc. C------------------------------------------------------------------ 300 IF (JINT .EQ. 0) THEN RESULT = (HALF - RESULT) + HALF IF (X .LT. ZERO) RESULT = -RESULT ELSE IF (JINT .EQ. 1) THEN IF (X .LT. ZERO) RESULT = TWO - RESULT ELSE IF (X .LT. ZERO) THEN IF (X .LT. XNEG) THEN RESULT = XINF ELSE YSQ = AINT(X*SIXTEN)/SIXTEN DEL = (X-YSQ)*(X+YSQ) Y = EXP(YSQ*YSQ) * EXP(DEL) RESULT = (Y+Y) - RESULT END IF END IF END IF 800 RETURN END SUBROUTINE CALERF END MODULE ERFCOREMOD MODULE DUNNETMOD SUBROUTINE MVNPRD(A, B, BPD, EPS, N, INF, IERC, HINC, PROB, & BOUND,IFAULT) C C ALGORITHM AS 251.1 APPL.STATIST. (1989), VOL.38, NO.3 C C FOR A MULTIVARIATE NORMAL VECTOR WITH CORRELATION STRUCTURE C DEFINED BY RHO(I,J) = BPD(I) * BPD(J), COMPUTES THE PROBABILITY C THAT THE VECTOR FALLS IN A RECTANGLE IN N-SPACE WITH ERROR C LESS THAN EPS. C INTEGER NN PARAMETER (NN = 50) REAL A(*), B(*), BPD(*), ESTT(22), FV(5), FD(5), F1T(22), * F2T(22), F3T(22), G1T(22), G3T(22), PSUM(22), H(NN), HL(NN), * BB(NN) INTEGER INF(*), INFT(NN), LDIR(22) REAL ZERO, HALF, ONE, TWO, FOUR, SIX, PT1, PT24, ONEP5, * X2880, SMALL, DXMIN, SQRT2, PROB, ERRL, BI, START, * Z, HINC, ADDN, EPS2, EPS1, EPS, ZU, Z2, Z3, Z4, Z5, ZZ, * ERFAC, EL, EL1, BOUND, PART0, PART2, PART3, FUNC0, FUNC2, * FUNCN, WT, CONTRB, DLG, DX, DA, ESTL, ESTR, SUM, EXCESS, ERROR, * PROB1, SAFE INTEGER N, IERC, IFAULT, I, NTM, NMAX, LVL, NR, NDIM REAL ALNORM, PPND7 EXTERNAL ALNORM, PPND7 DATA ZERO, HALF, ONE, TWO, FOUR, SIX /0.0, 0.5, 1.0, 2.0, * 4.0, 6.0/ DATA PT1, PT24, ONEP5, X2880 /0.1, 0.24, 1.5, 2880.0/ DATA SMALL, DXMIN, SQRT2 /1.0E-10, 0.0000001, 1.41421356237310/ C C CHECK FOR INPUT VALUES OUT OF RANGE. C PROB = ZERO BOUND = ZERO IFAULT = 1 IF (N .LT. 1 .OR. N .GT. NN) RETURN DO 10 I = 1, N BI = ABS(BPD(I)) IFAULT = 2 IF (BI .GE. ONE) RETURN IFAULT = 3 IF (INF(I) .LT. 0 .OR. INF(I) .GT. 2) RETURN IFAULT = 4 IF (INF(I) .EQ. 2 .AND. A(I) .LE. B(I)) RETURN 10 CONTINUE IFAULT = 0 PROB = ONE C C CHECK WHETHER ANY BPD(I) = 0. C NDIM = 0 DO 20 I = 1, N IF (BPD(I) .NE. ZERO) THEN NDIM = NDIM + 1 H(NDIM) = A(I) HL(NDIM) = B(I) BB(NDIM) = BPD(I) INFT(NDIM) = INF(I) ELSE C C IF ANY BPD(I) = 0, THE CONTRIBUTION TO PROB FOR THAT C VARIABLE IS COMPUTED FROM A UNIVARIATE NORMAL. C IF (INF(I) .LT. 1) THEN PROB = PROB * (ONE - ALNORM(B(I), .FALSE.)) ELSE IF (INF(I) .EQ. 1) THEN PROB = PROB * ALNORM(A(I), .FALSE.) ELSE PROB = PROB * (ALNORM(A(I), .FALSE.) - * ALNORM(B(I), .FALSE.)) END IF IF (PROB .LE. SMALL) PROB = ZERO END IF 20 CONTINUE IF (NDIM .EQ. 0 .OR. PROB .EQ. ZERO) RETURN C C IF NOT ALL BPD(I) = 0, PROB IS COMPUTED BY SIMPSON'S RULE. C BUT FIRST, INITIALIZE THE VARIABLES. C Z = ZERO IF (HINC .LE. ZERO) HINC = PT24 ADDN = -ONE DO 30 I = 1, NDIM IF (INFT(I) .EQ. 2 .OR. * (INFT(I) .NE. INFT(1) .AND. BB(I) * BB(1) .GT. ZERO) .OR. * (INFT(I) .EQ. INFT(1) .AND. BB(I) * BB(1) .LT. ZERO)) * ADDN = ZERO 30 CONTINUE C C THE VALUE OF ADDN IS TO BE ADDED TO THE PRODUCT EXPRESSIONS IN C THE INTEGRAND TO INSURE THAT THE LIMITING VALUE IS ZERO. C PROB1 = ZERO NTM = 0 NMAX = 400 IF (IERC .EQ. 0) NMAX = NMAX * 2 CALL PFUNC (Z, H, HL, BB, NDIM, INFT, ADDN, SAFE, FUNC0, NTM, * IERC, PART0) EPS2 = EPS * PT1 * HALF C C SET UPPER BOUND ON Z AND APPORTION EPS. C ZU = -PPND7(EPS2, IFAULT) / SQRT2 IF (IFAULT .NE. 0) THEN IFAULT = 6 RETURN END IF NR = IFIX(ZU / HINC) + 1 ERFAC = ONE IF (IERC .NE. 0) ERFAC = X2880 / HINC ** 5 EL = (EPS - EPS2) / FLOAT(NR) * ERFAC EL1 = EL C C START COMPUTATIONS FOR THE INTERVAL (Z, Z + HINC). C 40 ERROR = ZERO LVL = 0 FV(1) = PART0 FD(1) = SAFE START = Z DA = HINC Z3 = START + HALF * DA CALL PFUNC(Z3, H, HL, BB, NDIM, INFT, ADDN, FD(3), FUNCN, NTM, * IERC, FV(3)) Z5 = START + DA CALL PFUNC(Z5, H, HL, BB, NDIM, INFT, ADDN, FD(5), FUNC2, NTM, * IERC, FV(5)) PART2 = FV(5) SAFE = FD(5) WT = DA / SIX CONTRB = WT * (FV(1) + FOUR * FV(3) + FV(5)) DLG = ZERO IF (IERC .NE. 0) THEN CALL WMAX(FD(1), FD(3), FD(5), DLG) IF (DLG .LE. EL) GO TO 90 DX = DA GO TO 60 END IF LVL = 1 LDIR(LVL) = 2 PSUM(LVL) = ZERO C C BISECT INTERVAL. IF IERC = 1, COMPUTE ESTIMATE ON LEFT C HALF; IF IERC = 0, ON BOTH HALVES. C 50 DX = HALF * DA WT = DX / SIX Z2 = START + HALF * DX CALL PFUNC(Z2, H, HL, BB, NDIM, INFT, ADDN, FD(2), FUNCN, NTM, * IERC,FV(2)) ESTL = WT * (FV(1) + FOUR * FV(2) + FV(3)) IF (IERC .EQ. 0) THEN Z4 = START + ONEP5 * DX CALL PFUNC(Z4, H, HL, BB, NDIM, INFT, ADDN, FD(4), FUNCN, * NTM, IERC, FV(4)) ESTR = WT * (FV(3) + FOUR * FV(4) + FV(5)) SUM = ESTL + ESTR DLG = ABS(CONTRB - SUM) EPS1 = EL / TWO ** (LVL - 1) ERRL = DLG ELSE FV(3) = FV(2) FD(3) = FD(2) CALL WMAX(FD(1), FD(3), FD(5), DLG) ERRL = DLG / TWO ** (5 * LVL) SUM = ESTL EPS1 = EL * (TWO ** LVL) ** 4 END IF C C STOP SUBDIVIDING INTERVAL WHEN ACCURACY IS SUFFICIENT, C OR IF INTERVAL TOO NARROW OR SUBDIVIDED TOO OFTEN. C IF (DLG .LE. EPS1 .OR. DLG .LT. SMALL) GO TO 70 IF (IFAULT .EQ. 0 .AND. NTM .GE. NMAX) IFAULT = 5 IF (ABS(DX) .LE. DXMIN .OR. LVL .GT. 21) IFAULT = 7 IF (IFAULT .NE. 0) GO TO 70 C C RAISE LEVEL. STORE INFORMATION FOR RIGHT HALF AND APPLY C SIMPSON'S RULE TO LEFT HALF. C 60 LVL = LVL + 1 LDIR(LVL) = 1 F1T(LVL) = FV(3) F3T(LVL) = FV(5) DA = DX FV(5) = FV(3) IF (IERC .EQ. 0) THEN F2T(LVL) = FV(4) ESTT(LVL) = ESTR CONTRB = ESTL FV(3) = FV(2) ELSE G1T(LVL) = FD(3) G3T(LVL) = FD(5) FD(5) = FD(3) END IF GO TO 50 C C ACCEPT APPROXIMATE VALUE FOR INTERVAL. C RESTORE SAVED INFORMATION TO PROCESS C RIGHT HALF INTERVAL. C 70 ERROR = ERROR + ERRL 80 IF (LDIR(LVL) .EQ. 1) THEN PSUM(LVL) = SUM LDIR(LVL) = 2 IF (IERC .EQ. 0) DX = DX * TWO START = START + DX DA = HINC / TWO ** (LVL - 1) FV(1) = F1T(LVL) IF (IERC .EQ. 0) THEN FV(3) = F2T(LVL) CONTRB = ESTT(LVL) ELSE FV(3) = F3T(LVL) FD(1) = G1T(LVL) FD(5) = G3T(LVL) END IF FV(5) = F3T(LVL) GO TO 50 END IF SUM = SUM + PSUM(LVL) LVL = LVL - 1 IF (LVL .GT. 0) GO TO 80 CONTRB = SUM LVL = 1 DLG = ERROR 90 PROB1 = PROB1 + CONTRB BOUND = BOUND + DLG EXCESS = EL - DLG EL = EL1 IF (EXCESS .GT. ZERO) EL = EL1 + EXCESS IF ((FUNC0 .GT. ZERO .AND. FUNC2 .LE. FUNC0) .OR. * (FUNC0 .LT. ZERO .AND. FUNC2 .GE. FUNC0)) THEN ZZ = -SQRT2 * Z5 PART3 = ABS(FUNC2) * ALNORM(ZZ, .FALSE.) + BOUND / ERFAC IF (PART3 .LE. EPS .OR. NTM .GE. NMAX .OR. Z5 .GE. ZU) GOTO 100 END IF Z = Z5 PART0 = PART2 FUNC0 = FUNC2 IF (Z .LT. ZU .AND. NTM .LT. NMAX) GO TO 40 100 PROB = (PROB1 - ADDN * HALF) * PROB BOUND = PART3 IF (NTM .GE. NMAX .AND. IFAULT .EQ. 0) IFAULT = 5 IF (BOUND .GT. EPS .AND. IFAULT .EQ. 0) IFAULT = 8 RETURN END SUBROUTINE PFUNC(Z, A, B, BPD, N, INF, ADDN, DERIV, FUNCN, NTM, * IERC, RESULT) C C ALGORITHM AS 251.2 APPL.STATIST. (1989), VOL.38, NO.3 C C C COMPUTE FUNCTION IN INTEGRAND AND ITS 4TH DERIVATIVE. C INTEGER NN PARAMETER (NN = 50) REAL A(*), B(*), BPD(*), FOU(NN), FOU1(4, NN), TMP(4), GOU(NN), * GOU1(4, NN), FF(4), GF(4), TERM(4), GERM(4) INTEGER INF(*) REAL ZERO, ONE, TWO, THREE, FOUR, SIX, EIGHT, TWELVE, SIXTN, * SMALL, Z, U, U1, U2, BI, HI, HLI, BP, ADDN, DERIV, FUNCN, * RESULT, RSLT1, RSLT2, DEN, SQRT2, SQRTPI, PHI, PHI1, PHI2, * PHI3, PHI4, FRM, GRM INTEGER N, NTM, IERC, INFI, I, J, K, M, L, IK REAL ALNORM EXTERNAL ALNORM DATA ZERO, ONE, TWO, THREE, FOUR, SIX, EIGHT, TWELVE, SIXTN, * SMALL /0.0, 1.0, 2.0, 3.0, 4.0, 6.0, 8.0, 12.0, 16.0, 0.1E-12/ DATA SQRT2, SQRTPI /1.41421356237310, 1.77245385090552/ DERIV = ZERO NTM = NTM + 1 RSLT1 = ONE RSLT2 = ONE BI = ONE HI = A(1) + ONE HLI = B(1) + ONE INFI = -1 DO 60 I = 1, N IF (BPD(I) .EQ. BI .AND. A(I) .EQ. HI .AND. B(I) .EQ. HLI .AND. * INF(I) .EQ. INFI) THEN FOU(I) = FOU(I - 1) GOU(I) = GOU(I - 1) DO 10 IK = 1, 4 FOU1(IK, I) = FOU1(IK, I - 1) GOU1(IK, I) = GOU1(IK, I - 1) 10 CONTINUE ELSE BI = BPD(I) HI = A(I) HLI = B(I) INFI = INF(I) IF (BI .EQ. ZERO) THEN IF (INFI .LT. 1) THEN FOU(I) = ONE - ALNORM(HLI, .FALSE.) ELSE IF (INFI .EQ. 1) THEN FOU(I) = ALNORM(HI, .FALSE.) ELSE FOU(I) = ALNORM(HI, .FALSE.) - ALNORM(HLI, .FALSE.) END IF GOU(I) = FOU(I) DO 20 IK = 1, 4 FOU1(IK, I) = ZERO GOU1(IK, I) = ZERO 20 CONTINUE ELSE DEN = SQRT(ONE - BI * BI) BP = BI * SQRT2 / DEN IF (INFI .LT. 1) THEN U = -HLI / DEN + Z * BP FOU(I) = ALNORM(U, .FALSE.) CALL ASSIGN (U, BP, FOU1(1, I)) BP = -BP U = -HLI / DEN + Z * BP GOU(I) = ALNORM(U, .FALSE.) CALL ASSIGN (U, BP, GOU1(1, I)) ELSE IF (INFI .EQ. 1) THEN U = HI / DEN + Z * BP GOU(I) = ALNORM(U, .FALSE.) CALL ASSIGN (U, BP, GOU1(1, I)) BP = -BP U = HI / DEN + Z * BP FOU(I) = ALNORM(U, .FALSE.) CALL ASSIGN (U, BP, FOU1(1, I)) ELSE U2 = -HLI / DEN + Z * BP CALL ASSIGN (U2, BP, FOU1(1, I)) BP = -BP U1 = HI / DEN + Z * BP CALL ASSIGN (U1, BP, TMP(1)) FOU(I) = ALNORM(U1, .FALSE.) + ALNORM(U2, .FALSE.) - ONE DO 30 IK = 1, 4 FOU1(IK, I) = FOU1(IK, I) + TMP(IK) 30 CONTINUE IF (-HLI .EQ. HI) THEN GOU(I) = FOU(I) DO 40 IK = 1, 4 GOU1(IK, I) = FOU1(IK, I) 40 CONTINUE ELSE U2 = -HLI / DEN + Z * BP CALL ASSIGN (U2, BP, GOU1(1, I)) BP = -BP U1 = HI / DEN + Z * BP GOU(I) = ALNORM(U1, .FALSE.) + ALNORM(U2, .FALSE.)-ONE CALL ASSIGN (U1, BP, TMP(1)) DO 50 IK = 1, 4 GOU1(IK, I) = GOU1(IK, I) + TMP(IK) 50 CONTINUE END IF END IF END IF END IF RSLT1 = RSLT1 * FOU(I) RSLT2 = RSLT2 * GOU(I) IF (RSLT1 .LE. SMALL) RSLT1 = ZERO IF (RSLT2 .LE. SMALL) RSLT2 = ZERO 60 CONTINUE FUNCN = RSLT1 + RSLT2 + ADDN RESULT = FUNCN * EXP(-Z * Z) / SQRTPI C C IF 4TH DERIVATIVE IS NOT WANTED, STOP HERE. C OTHERWISE, PROCEED TO COMPUTE 4TH DERIVATIVE. C IF (IERC .EQ. 0) RETURN DO 70 IK = 1, 4 FF(IK) = ZERO GF(IK) = ZERO 70 CONTINUE DO 100 I = 1, N FRM = ONE GRM = ONE DO 80 J = 1, N IF (J .EQ. 1) GO TO 80 FRM = FRM * FOU(J) GRM = GRM * GOU(J) IF (FRM .LE. SMALL) FRM = ZERO IF (GRM .LE. SMALL) GRM = ZERO 80 CONTINUE DO 90 IK = 1, 4 FF(IK) = FF(IK) + FRM * FOU1(IK, I) GF(IK) = GF(IK) + GRM * GOU1(IK, I) 90 CONTINUE 100 CONTINUE IF (N .LE. 2) GO TO 230 DO 130 I = 1, N DO 120 J = I + 1, N TERM(2) = FOU1(1, I) * FOU1(1, J) GERM(2) = GOU1(1, I) * GOU1(1, J) TERM(3) = FOU1(2, I) * FOU1(1, J) GERM(3) = GOU1(2, I) * GOU1(1, J) TERM(4) = FOU1(3, I) * FOU1(1, J) GERM(4) = GOU1(3, I) * GOU1(1, J) TERM(1) = FOU1(2, I) * FOU1(2, J) GERM(1) = GOU1(2, I) * GOU1(2, J) DO 110 K = 1, N IF (K .EQ. I .OR. K .EQ. J) GO TO 110 CALL TOOSML (1, TERM, FOU(K)) CALL TOOSML (1, GERM, GOU(K)) 110 CONTINUE FF(2) = FF(2) + TWO * TERM(2) FF(3) = FF(3) + TWO * TERM(3) * THREE FF(4) = FF(4) + TWO * (TERM(4) * FOUR + TERM(1) * THREE) GF(2) = GF(2) + TWO * GERM(2) GF(3) = GF(3) + TWO * GERM(3) * THREE GF(4) = GF(4) + TWO * (GERM(4) * FOUR + GERM(1) * THREE) 120 CONTINUE 130 CONTINUE DO 170 I = 1, N DO 160 J = I + 1, N DO 150 K = J + 1, N TERM(3) = FOU1(1, I) * FOU1(1, J) * FOU1(1, K) TERM(4) = FOU1(2, I) * FOU1(1, J) * FOU1(1, K) GERM(3) = GOU1(1, I) * GOU1(1, J) * GOU1(1, K) GERM(4) = GOU1(2, I) * GOU1(1, J) * GOU1(1, K) IF (N .GT. 3) THEN DO 140 M = 1, N IF (M .EQ. I .OR. M .EQ. J .OR. M .EQ. K) GO TO 140 CALL TOOSML (3, TERM, FOU(M)) CALL TOOSML (3, GERM, GOU(M)) 140 CONTINUE END IF FF(3) = FF(3) + SIX * TERM(3) FF(4) = FF(4) + SIX * TERM(4) * SIX GF(3) = GF(3) + SIX * GERM(3) GF(4) = GF(4) + SIX * GERM(4) * SIX 150 CONTINUE 160 CONTINUE 170 CONTINUE IF (N .LE. 3) GO TO 230 DO 220 I = 1, N DO 210 J = I + 1, N DO 200 K = J + 1, N DO 190 M = K + 1, N TERM(4) = FOU1(1, I) * FOU1(1, J) * FOU1(1, K) * FOU1(1, M) GERM(4) = GOU1(1, I) * GOU1(1, J) * GOU1(1, K) * GOU1(1, M) IF (N .GT. 4) THEN DO 180 L = 1, N IF (L .EQ. I .OR. L .EQ. J .OR. L .EQ. K .OR. L .EQ. M)GOTO 180 CALL TOOSML (4, TERM, FOU(L)) CALL TOOSML (4, GERM, GOU(L)) 180 CONTINUE END IF FF(4) = FF(4) + FOUR * SIX * TERM(4) GF(4) = GF(4) + FOUR * SIX * GERM(4) 190 CONTINUE 200 CONTINUE 210 CONTINUE 220 CONTINUE C 230 CONTINUE PHI = EXP(-Z * Z) / SQRTPI PHI1 = -TWO * Z * PHI PHI2 = (FOUR * Z ** 2 - TWO) * PHI PHI3 = (-EIGHT * Z ** 3 + TWELVE * Z) * PHI PHI4 = (SIXTN * Z ** 2 * (Z ** 2 - THREE) + TWELVE) * PHI DERIV = PHI * (FF(4) + GF(4)) + FOUR * PHI1 * (FF(3) + GF(3)) * + SIX * PHI2 * (FF(2) + GF(2)) + FOUR * PHI3 * (FF(1) + GF(1)) * + PHI4 * FUNCN RETURN END SUBROUTINE ASSIGN (U, BP, FF) C C ALGORITHM AS 251.3 APPL.STATIST. (1989), VOL.38, NO.3 C C C COMPUTE DERIVATIVES OF NORMAL CDF'S. C REAL FF(4) REAL U, U2, BP, HALF, ONE, THREE, SQ2PI, T1, T2, T3 INTEGER I DATA HALF, ONE, THREE, SQ2PI /0.5, 1.0, 3.0, 2.50662827463100/ DATA ZERO, UMAX, SMALL /0.0, 8.0, 0.1E-07/ IF (ABS(U) .GT. UMAX) THEN DO 10 I = 1, 4 FF(I) = ZERO 10 CONTINUE ELSE U2 = U * U T1 = BP * EXP(-HALF * U2) / SQ2PI T2 = BP * T1 T3 = BP * T2 FF(1) = T1 FF(2) = -U * T2 FF(3) = (U2 - ONE) * T3 FF(4) = (THREE - U2) * U * BP * T3 DO 20 I = 1, 4 IF(ABS(FF(I)) .LT. SMALL) FF(I) = ZERO 20 CONTINUE END IF RETURN END SUBROUTINE WMAX(W1, W2, W3, DLG) C C ALGORITHM AS 251.4 APPL.STATIST. (1989), VOL.38, NO.3 C C C LARGEST ABSOLUTE VALUE OF QUADRATIC FUNCTION FITTED C TO THREE POINTS. C REAL W1, W2, W3, DLG, QUAD, QLIM, QMIN, ONE, TWO, B2C DATA ONE, TWO, QMIN /1.0, 2.0, 0.00001/ DLG = MAX( ABS(W1), ABS(W3) ) QUAD = W1 - W2 * TWO + W3 QLIM = MAX( ABS(W1 - W3) / TWO , QMIN) IF (ABS(QUAD) .LE. QLIM) RETURN B2C = (W1 - W3) / QUAD / TWO IF (ABS(B2C) .GE. ONE) RETURN DLG = MAX( DLG, ABS(W2 - B2C * QUAD * B2C / TWO) ) RETURN END SUBROUTINE TOOSML (N, FF, F) C C ALGORITHM AS 251.5 APPL.STATIST. (1989), VOL.38, NO.3 C C C MULTIPLY FF(I) BY F FOR I = N TO 4. SET TO ZERO IF TOO SMALL. C REAL FF(4), F, ZERO, SMALL INTEGER N, I DATA ZERO, SMALL /0.0, 0.1E-12/ DO 10 I = N, 4 FF(I) = FF(I) * F IF (ABS(FF(I)) .LE. SMALL) FF(I) = ZERO 10 CONTINUE RETURN END REAL FUNCTION ALNORM(X, UPPER) C C ALGORITHM AS 66 APPL. STATIST. (1973) VOL.22, P.424 C C EVALUATES THE TAIL AREA OF THE STANDARDIZED NORMAL CURVE C FROM X TO INFINITY IF UPPER IS .TRUE. OR C FROM MINUS INFINITY TO X IF UPPER IS .FALSE. C REAL LTONE, UTZERO, ZERO, HALF, ONE, CON, A1, A2, A3, $ A4, A5, A6, A7, B1, B2, B3, B4, B5, B6, B7, B8, B9, $ B10, B11, B12, X, Y, Z, ZEXP LOGICAL UPPER, UP C C LTONE AND UTZERO MUST BE SET TO SUIT THE PARTICULAR COMPUTER C (SEE INTRODUCTORY TEXT) C DATA LTONE, UTZERO /7.0, 18.66/ DATA ZERO, HALF, ONE, CON /0.0, 0.5, 1.0, 1.28/ DATA A1, A2, A3, $ A4, A5, A6, $ A7 $ /0.398942280444, 0.399903438504, 5.75885480458, $ 29.8213557808, 2.62433121679, 48.6959930692, $ 5.92885724438/ DATA B1, B2, B3, $ B4, B5, B6, $ B7, B8, B9, $ B10, B11, B12 $ /0.398942280385, 3.8052E-8, 1.00000615302, $ 3.98064794E-4, 1.98615381364, 0.151679116635, $ 5.29330324926, 4.8385912808, 15.1508972451, $ 0.742380924027, 30.789933034, 3.99019417011/ C ZEXP(Z) = EXP(Z) C UP = UPPER Z = X IF (Z .GE. ZERO) GOTO 10 UP = .NOT. UP Z = -Z 10 IF (Z .LE. LTONE .OR. UP .AND. Z .LE. UTZERO) GOTO 20 ALNORM = ZERO GOTO 40 20 Y = HALF * Z * Z IF (Z .GT. CON) GOTO 30 C ALNORM = HALF - Z * (A1 - A2 * Y / (Y + A3 - A4 / (Y + A5 + $ A6 / (Y + A7)))) GOTO 40 C 30 ALNORM = B1 * ZEXP(-Y) / (Z - B2 + B3 / (Z + B4 + B5 / (Z - $ B6 + B7 / (Z + B8 - B9 / (Z + B10 + B11 / (Z + B12)))))) C 40 IF (.NOT. UP) ALNORM = ONE - ALNORM RETURN END REAL FUNCTION PPND7 (P, IFAULT) C C ALGORITHM AS241 APPL. STATIST. (1988) VOL. 37, NO. 3 C C PRODUCES THE NORMAL DEVIATE Z CORRESPONDING TO A GIVEN LOWER C TAIL AREA OF P; Z IS ACCURATE TO ABOUT 1 PART IN 10**7. C C THE HASH SUMS BELOW ARE THE SUMS OF THE MANTISSAS OF THE C COEFFICIENTS. THEY ARE INCLUDED FOR USE IN CHECKING C TRANSCRIPTION. C INTEGER IFAULT REAL ZERO, ONE, HALF, SPLIT1, SPLIT2, CONST1, CONST2, * A0, A1, A2, A3, B1, B2, B3, C0, C1, C2, C3, D1, D2, * E0, E1, E2, E3, F1, F2, P, Q, R PARAMETER (ZERO = 0.0E0, ONE = 1.0E0, HALF = 0.5E0, * SPLIT1 = 0.425E0, SPLIT2 = 5.0E0, * CONST1 = 0.180625E0, CONST2 = 1.6E0) C C COEFFICIENTS FOR P CLOSE TO 1/2 PARAMETER (A0 = 3.38713 27179E0, * A1 = 5.04342 71938E1, * A2 = 1.59291 13202E2, * A3 = 5.91093 74720E1, * B1 = 1.78951 69469E1, * B2 = 7.87577 57664E1, * B3 = 6.71875 63600E1) C HASH SUM AB 32.31845 77772 C C COEFFICIENTS FOR P NEITHER CLOSE TO 1/2 NOR 0 OR 1 PARAMETER (C0 = 1.42343 72777E0, * C1 = 2.75681 53900E0, * C2 = 1.30672 84816E0, * C3 = 1.70238 21103E-1, * D1 = 7.37001 64250E-1, * D2 = 1.20211 32975E-1) C HASH SUM CD 15.76149 29821 C C COEFFICIENTS FOR P NEAR 0 OR 1 PARAMETER (E0 = 6.65790 51150E0, * E1 = 3.08122 63860E0, * E2 = 4.28682 94337E-1, * E3 = 1.73372 03997E-2, * F1 = 2.41978 94225E-1, * F2 = 1.22582 02635E-2) C HASH SUM EF 19.40529 10204 C IFAULT = 0 Q = P - HALF IF (ABS(Q) .LE. SPLIT1) THEN R = CONST1 - Q * Q PPND7 = Q * (((A3 * R + A2) * R + A1) * R + A0) / * (((B3 * R + B2) * R + B1) * R + ONE) RETURN ELSE IF (Q .LT. 0) THEN R = P ELSE R = ONE - P ENDIF IF (R .LE. ZERO) THEN IFAULT = 1 PPND7 = ZERO RETURN ENDIF R = SQRT(-LOG(R)) IF (R .LE. SPLIT2) THEN R = R - CONST2 PPND7 = (((C3 * R + C2) * R + C1) * R + C0) / * ((D2 * R + D1) * R + ONE) ELSE R = R - SPLIT2 PPND7 = (((E3 * R + E2) * R + E1) * R + E0) / * ((F2 * R + F1) * R + ONE) ENDIF IF (Q .LT. 0) PPND7 = -PPND7 RETURN ENDIF END SUBROUTINE SIMPSN (NDF,A,B,BPD,ERRB,N,INF,D,IERC,HNC,PROB, * BND,IFLT) C C STUDENTIZES A MULTIVARIATE INTEGRAL USING SIMPSON'S RULE. C DIMENSION A(*),B(*),BPD(*),INF(*),D(*), * FV(5),F1T(30),F2T(30),F3T(30), * LDIR(30),PSUM(30),ESTT(30),ERRR(30),GV(5),G1T(30),G2T(30), * G3T(30),GSUM(30) DATA ZERO,HALF,ONE,ONEP5,TWO,FOUR,SIX,DXMIN /0.0,0.5,1.0,1.5, * 2.0,4.0,6.0,0.000004/ PROB = ZERO BOUNDA = ZERO BOUNDG = ZERO IFLAG = 0 IER = 0 START = -ONE DAX = ONE ERB2 = ERRB * HALF EPS1 = ERB2 * HALF CALL FUN (ZERO,NDF,A,B,BPD,ERB2,N,INF,D,F0,G0,IERC,HNC,IER) 10 FV(1) = ZERO GV(1) = ZERO ERROR = ZERO DA = DAX LVL = 1 Z3 = START + HALF*DA CALL FUN(Z3,NDF,A,B,BPD,ERB2,N,INF,D,FV(3),GV(3),IERC,HNC,IER) FV(5) = F0 GV(5) = G0 WT = ABS(DA) / SIX CONTRB = WT * (FV(1) + FOUR * FV(3) + FV(5)) CONTRG = WT * (GV(1) + FOUR * GV(3) + GV(5)) LDIR(LVL) = 2 PSUM(LVL) = ZERO GSUM(LVL) = ZERO C C BISECT INTERVAL; COMPUTE ESTIMATES FOR EACH HALF. C 20 DX = HALF * DA WT = ABS(DX) / SIX Z2 = START + HALF * DX CALL FUN(Z2,NDF,A,B,BPD,ERB2,N,INF,D,FV(2),GV(2),IERC,HNC,IER) Z4 = START + ONEP5 * DX CALL FUN(Z4,NDF,A,B,BPD,ERB2,N,INF,D,FV(4),GV(4),IERC,HNC,IER) ESTL = WT * (FV(1) + FOUR * FV(2) + FV(3)) ESTR = WT * (FV(3) + FOUR * FV(4) + FV(5)) ESTGL = WT * (GV(1) + FOUR * GV(2) + GV(3)) ESTGR = WT * (GV(3) + FOUR * GV(4) + GV(5)) SUM = ESTL + ESTR SUMG = ESTGL + ESTGR DLG = ABS(CONTRB - SUM) ERRL = DLG C C STOP BISECTING WHEN ACCURACY SUFFICIENT, OR IF C INTERVAL TOO NARROW OR BISECTED TOO OFTEN. C 30 IF (DLG .LE. EPS1) GO TO 50 IF (ABS(DX) .LE. DXMIN .OR. LVL .GE. 30) GO TO 40 C C RAISE LEVEL. STORE INFORMATION FOR RIGHT HALF C AND APPLY SIMPSON'S RULE TO LEFT HALF. C LVL = LVL + 1 LDIR(LVL) = 1 F1T(LVL) = FV(3) F2T(LVL) = FV(4) F3T(LVL) = FV(5) G1T(LVL) = GV(3) G2T(LVL) = GV(4) G3T(LVL) = GV(5) DA = DX FV(5) = FV(3) FV(3) = FV(2) GV(5) = GV(3) GV(3) = GV(2) ESTT(LVL) = ESTR CONTRB = ESTL CONTRG = ESTGL EPS1 = EPS1 * HALF ERRR(LVL) = EPS1 GO TO 20 C C ACCEPT APPROXIMATE VALUE FOR INTERVAL. C 40 IFLAG = 11 50 ERROR = ERROR + ERRL 60 IF (LDIR(LVL) .EQ. 1) GO TO 70 SUM = SUM + PSUM(LVL) SUMG = SUMG + GSUM(LVL) LVL = LVL - 1 IF (LVL .GT. 0) GO TO 60 CONTRB = SUM CONTRG = SUMG LVL = 1 DLG = ERROR GO TO 80 C C RESTORE SAVED INFORMATION TO PROCESS RIGHT HALF. C 70 PSUM(LVL) = SUM GSUM(LVL) = SUMG LDIR(LVL) = 2 DA = DAX / TWO**(LVL-1) START = START + DX * TWO FV(1) = F1T(LVL) FV(3) = F2T(LVL) FV(5) = F3T(LVL) GV(1) = G1T(LVL) GV(3) = G2T(LVL) GV(5) = G3T(LVL) CONTRB = ESTT(LVL) EXCESS = EPS1 - DLG EPS1 = ERRR(LVL) IF (EXCESS .GT. ZERO) EPS1 = EPS1 + EXCESS GO TO 20 80 PROB = PROB + CONTRB BOUNDG = BOUNDG + CONTRG BOUNDA = BOUNDA + DLG IF (Z4 .LE. ZERO) GO TO 90 IF (IFLT .EQ. 0) IFLT = IER IF (IFLT .EQ. 0) IFLT = IFLAG BOUNDA = BOUNDA + BOUNDG IF (BND .LT. BOUNDA) BND = BOUNDA RETURN 90 EPS1 = ERB2 * HALF EXCESS = EPS1 - BND IF (EXCESS .GT. ZERO) EPS1 = EPS1 + EXCESS START = ONE DAX = -ONE GO TO 10 END FUNCTION SDIST(Y,N) C C COMPUTE Y**(N/2 - 1) EXP(-Y) / GAMMA(N/2) C C (Revised: 1994-01-19) C DATA ZERO, HALF, ONE, X23 / 0.0, 0.5, 1.0, -23.0 / DATA SQRTPI / 1.77245385090552 / SDIST = ZERO IF (Y .LE. ZERO) RETURN JJ = N/2 - 1 JK = 2 * JJ - N + 2 JKP = JJ - JK SDIST = ONE IF (JK .LT. 0) SDIST = SDIST / SQRT(Y) / SQRTPI IF (JKP .EQ. 0) GO TO 20 XN = FLOAT(N) * HALF TEST = ALOG(Y) - Y / FLOAT(JKP) IF ( TEST .LT. X23 ) THEN SDIST = ZERO RETURN ENDIF SDIST = ALOG ( SDIST ) DO 10 J = 1, JKP XN = XN - ONE SDIST = SDIST + TEST - ALOG(XN) 10 CONTINUE IF ( SDIST .LT. X23 ) THEN SDIST = ZERO ELSE SDIST = EXP( SDIST ) ENDIF RETURN 20 SDIST = SDIST * EXP(-Y) RETURN END SUBROUTINE FUN (Z,NDF,H,HL,BPD,ERB2,N,INF,D,F0,G0,IERC * ,HNC,IER) INTEGER NN PARAMETER (NN=50) DIMENSION A(NN),B(NN),H(*),HL(*),BPD(*),INF(*),D(*) DATA ZERO, ONE, TWO, SMALL / 0.0, 1.0, 2.0, 1.0E-08 / F0 = ZERO G0 = ZERO IF (Z .LE. -ONE .OR. Z .GE. ONE) RETURN DF = FLOAT(NDF) ARG = (ONE + Z) / (ONE - Z) TERM = ARG * DF * TWO / (ONE-Z)**2 * SDIST(DF/TWO*ARG*ARG,NDF) IF (TERM .LE. SMALL) RETURN DO 10 I = 1, N A(I) = ARG * H(I) - D(I) B(I) = ARG * HL(I) - D(I) 10 CONTINUE CALL MVNPRD (A,B,BPD,ERB2,N,INF,IERC,HNC,PROB,BND,IFLT) IF (IER .EQ. 0) IER = IFLT G0 = TERM * BND F0 = TERM * PROB RETURN END C * * * * * * * * * * * * * * * * * * * * * * * * * * * * C Charles Dunnett C Dept. of Mathematics and Statistics C McMaster University C Hamilton, Ontario L8S 4K1 C Canada C E-mail: dunnett@mcmaster.ca C Tel.: (905) 525-9140 (Ext. 27104) C * * * * * * * * * * * * * * * * * * * * * * * * * * * * END MODULE DUNNETMOD