C 
C f2py -m mvnprd -h mvnprd.pyf mvnprd.f  only: mvnprd
C  edit mvnprd.pyf with input and output and then
C f2py mvnprd.pyf mvnprd.f -c --fcompiler=gnu95 --compiler=mingw32 -lmsvcr71
C
C f2py --fcompiler=gnu95 --compiler=mingw32 -lmsvcr71 -m mvnprd  -c mvnprd.f 

C Altarnative: compile mvnprd and link to it through mvnprd_interface.f
C 
C gfortran -fPIC -c mvnprd.f
C f2py -m mvnprdmod  -c mvnprd.o mvnprd_interface.f --fcompiler=gnu95 --compiler=mingw32 -lmsvcr71
C
C df -c mvnprd.f
C f2py -m mvnprdmod  -c mvnprd.obj mvnprd_interface.f --fcompiler=compaqv --compiler=mingw32 -lmsvcr71

! This is a MEX-file for MATLAB.
!     and contains a mex-interface to Charles W. Dunnett's programs
!     ,MVNPRD and MVSTUD subroutines for computing multivariate normal
!     or student T probabilities with product correlation structure. The
!     file should compile without errors on (Fortran77) standard Fortran
!     compilers. 
*
* The mex-interface was written by
*     Per Andreas Brodtkorb
*     Norwegian Defence Research Establishment
*     P.O. Box 115m
*     N-3191 Horten
*     Norway
*     Email: Per.Brodtkorb@ffi.no
*
*     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)
*
* MVNPRDMEX Computes multivariate normal  or student T probability
*           with product correlation structure.

*
*  CALL [value,bound,inform] = mvnprdmex(RHO,A,B,D,NDF,abseps,IERC,HNC)
*
*     RHO    REAL, array of coefficients defining the correlation
*            coefficient by:
*                correlation(I,J) =  RHO(I)*RHO(J) for J/=I
*            where 
*                1 < RHO(I) < 1
*     A		 REAL, array of lower integration limits.
*     B		 REAL, array of upper integration limits.
*	       NOTE: any values greater the 37, are considered as
*              infinite values.
*     D         Real array of means
*     NDF       Degrees of freedom, NDF<=0 gives normal probabilities
*     ABSEPS REAL absolute error tolerance.
*     IERC   INTEGER 1 if strict error control based on fourth
*                      derivative
*                    0 if intuitive error control based on halving the
*                      intervals
*     HINC  REAL start interval width of simpson rule
*
* OUTPUT:
*     VALUE  REAL estimated value for the integral
*     BOUND  REAL bound on the error of the approximation
*     INFORM INTEGER, termination status parameter:
*            0, if normal completion with ERROR < EPS;
*            1, if N > 100 or N < 1.
*            2, IF  any abs(rho)>=1      
*            4, if  ANY(B(I)<=A(i))
*            5, if number of terms computed exceeds maximum number of
*                  evaluation points
*            6, if fault accurs in normal subroutines
*            7, if subintervals are too narrow or too many
*            8, if bounds exceeds abseps            
*
*
* MVNPRDMEX calculates multivariate normal or student T probability
* with product correlation structure for rectangular regions.
* The accuracy is up to around single precision, i.e., about 1e-7.
*
* This file was successfully compiled for matlab 5.3
* using Compaq Visual Fortran 6.1, and Windows 2000.
* The example here uses Fortran77 source.
* First, you will need to modify your mexopts.bat file.
* To find it, issue the command prefdir(1) from the Matlab command line,
* the directory it answers with will contain your mexopts.bat file.
* Open it for editing. The first section will look like:
*
*rem ********************************************************************
*rem General parameters
*rem ********************************************************************
*set MATLAB=%MATLAB%
*set DF_ROOT=C:\Program Files\Microsoft Visual Studio
*set VCDir=%DF_ROOT%\VC98
*set MSDevDir=%DF_ROOT%\Common\msdev98
*set DFDir=%DF_ROOT%\DF98
*set PATH=%MSDevDir%\bin;%DFDir%\BIN;%VCDir%\BIN;%PATH%
*set INCLUDE=%DFDir%\INCLUDE;%DFDir%\IMSL\INCLUDE;%INCLUDE%
*set LIB=%DFDir%\LIB;%VCDir%\LIB
*
* then you are ready to compile this file at the matlab prompt using the
* following command:
*  mex -O mvnprdmex.f



C The rest of this file contains:
C 1. A Readme file provided by Charles Dunnett, the author of AS 251.
C 2. The published algorithm AS 251 together with the two other AS algorithms
C    which it calls (AS 66 and AS 241).
C 3. A driver program (MVTIN) for either multivariate normal or t.
C 4. MVSTUD for calculating multivariate t probabilities.
C ***************************************************************************

C Date: Mon, 10 Apr 1995 16:49:10 +0059 (EDT)
C From: "Charles W. Dunnett" <dunnett@mcmail.cis.mcmaster.ca>
C Subject: Readme for AS251 (extended version incl. multivariate t)

C MVTIN is a driver program for computing multivariate normal or t
C probability integrals over arbitrary rectangular regions.  The
C correlation structure is assumed to be of product form, rho_ij =
C b_i x b_j, where -1 < b_i < +1.

C It requires the following:-

C     1.	 MVNPRD  (published as algorithm AS 251 in Applied
C       Statistics (1989), 38: 564-579; see also the correction
C       note in Applied Statistics (1993), 42: 709),

C     2.	 ALNORM and PPND7 (published as algorithms AS 66 and
C       AS 241, respectively, in Applied Statistics, and

C     3.	 MVSTUD   ( which Studentizes MVNPRD).



      SUBROUTINE MVNPRD(A, B, BPD, EPS, N, INF, IERC, HINC, PROB, BOUND,
     *  IFAULT)
      implicit none
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 = 100)
      DOUBLE PRECISION 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)
      DOUBLE PRECISION 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, TSUM, EXCESS, ERROR,
     *  PROB1, SAFE, ONEP5,X2880
      INTEGER N, IERC, IFAULT, I, NTM, NMAX, LVL, NR, NDIM
      DOUBLE PRECISION 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
      NR = NINT(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))
         TSUM = ESTL + ESTR
         DLG = ABS(CONTRB - TSUM)
         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)
         TSUM = 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) = TSUM
         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
      TSUM = TSUM + PSUM(LVL)
      LVL = LVL - 1
      IF (LVL .GT. 0) GO TO 80
      CONTRB = TSUM
      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)
      implicit none
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 = 100)
      DOUBLE PRECISION 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(*)
      DOUBLE PRECISION 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
      DOUBLE PRECISION 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)
      implicit none
C
C        ALGORITHM AS 251.3  APPL.STATIST. (1989), VOL.38, NO.3
C
C
C        COMPUTE DERIVATIVES OF NORMAL CDF'S.
C
      DOUBLE PRECISION FF(4)
      DOUBLE PRECISION U, U2, BP, HALF, ONE, THREE, SQ2PI, T1, T2, T3
     $     ,ZERO, UMAX, SMALL
      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)
      implicit none
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
      DOUBLE PRECISION 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)
      implicit none
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
      DOUBLE PRECISION 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
      DOUBLE PRECISION FUNCTION ALNORM(X, UPPER)
      implicit none
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
      DOUBLE PRECISION 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
      DOUBLE PRECISION FUNCTION PPND7 (P, IFAULT)
      implicit none
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
      DOUBLE PRECISION 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 MVSTUD(NDF,A,B,BPD,ERRB,N,INF,D,IERC,HNC,PROB,
     *    BND,IFLT)
         implicit none
C
C        COMPUTE MULTIVARIATE STUDENT INTEGRAL,
C        USING MVNPRD (DUNNETT, APPL. STAT., 1989)
C        IF RHO(I,J) = BPD(I)*BPD(J).
C
C        IF RHO(I,J) HAS GENERAL STRUCTURE, USE
C        MULNOR (SCHERVISH, APPL. STAT., 1984) AND REPLACE
C        CALL MVNPRD(A,B,BPD,EPS,N,INF,IERC,HNC,PROB,BND,IFLT)
C        BY CALL MULNOR(A,B,SIG,EPS,N,INF,PROB,BND,IFLT).
C
C        AUTHOR: C.W. DUNNETT, MCMASTER UNVERSITY
C
C        BASED ON ADAPTIVE SIMPSON'S RULE ALGORITHM
C        DESCRIBED IN SHAMPINE & ALLEN: "NUMERICAL
C        COMPUTING", (1974), PAGE 240.
C
C        PARAMETERS ARE SAME AS IN ALGORITHM AS 251
C        IN APPL. STAT. (1989), VOL. 38: 564-579
C        WITH THE FOLLOWING ADDITIONS:
C             NDF   INTEGER      INPUT  DEGREES OF FREEDOM
C             D     REAL ARRAY   INPUT  NON-CENTRALITY VECTOR
C        (PUT NDF = 0 FOR INFINITE D.F.)
C
      DOUBLE PRECISION :: HNC,PROB,BND
      INTEGER :: NN, MAXDF,I,IERC,NDF,N,IFLT
      PARAMETER (NN=100, MAXDF = 150)
	integer :: INF(*)
      DOUBLE PRECISION :: A(*),B(*),BPD(*),D(*),F(3),
     &     AA(NN),BB(NN)
      DOUBLE PRECISION :: ERB2, ERRB, AX,BX,XX
      DOUBLE PRECISION,SAVE :: ZERO,HALF,TWO,THREE,FOUR
      INTEGER :: NF
      !DIMENSION A(*),B(*),BPD(*),INF(*),D(*),F(3),AA(NN),BB(NN)
      DATA ZERO,HALF,TWO,THREE,FOUR / 0.0, 0.5, 2.0, 3.0, 4.0 /
      !external float
      DO 10 I = 1, N
         AA(I) = A(I) - D(I)
         BB(I) = B(I) - D(I)
   10 CONTINUE
      IF (NDF .LE. 0) THEN
         CALL MVNPRD(AA,BB,BPD,ERRB,N,INF,IERC,HNC,PROB,BND,IFLT)
         RETURN
      ENDIF
      BND   = ZERO
      IFLT  =  0
      
      ERB2  = ERRB
C
C        CHECK IF D.F. EXCEED MAXDF; IF YES, THEN PROB
C        IS COMPUTED BY QUADRATIC INTERPOLATION ON 1./D.F.
C
      IF (NDF .LE. MAXDF) GO TO 20
      CALL MVNPRD(AA,BB,BPD,ERB2,N,INF,IERC,HNC,F(1),BND,IFLT)
      NF    =  MAXDF / 2
      CALL SIMPSN(NF,A,B,BPD,ERB2,N,INF,D,IERC,HNC,F(3),BND,IFLT)
      NF    =  NF * 2
      CALL SIMPSN(NF,A,B,BPD,ERB2,N,INF,D,IERC,HNC,F(2),BND,IFLT)
      XX    =  DBLE(NF) / DBLE(NDF)
      AX    =  F(3) - F(2)*TWO + F(1)
      BX    =  F(2)*FOUR - F(3) - F(1)*THREE
      PROB  =  F(1) + XX * (AX * XX + BX) * HALF
      RETURN
   20 CALL SIMPSN (NDF,A,B,BPD,ERB2,N,INF,D,IERC,HNC,PROB,BND,IFLT)
      RETURN
      END
      SUBROUTINE SIMPSN (NDF,A,B,BPD,ERRB,N,INF,D,IERC,HNC,PROB,
     *   BND,IFLT)
      implicit none
C
C        STUDENTIZES A MULTIVARIATE INTEGRAL USING SIMPSON'S RULE.
C
      double precision :: A,B,BPD,INF,D,
     *   FV,F1T,F2T,F3T,
     *   LDIR,PSUM,ESTT,ERRR,GV,G1T,G2T,
     *   G3T,GSUM
      double precision :: PROB, BOUNDA, BOUNDG,sTART, DAX, ERB2,ERRB,
     $     EPS1, F0,G0, HNC, ERROR,DA,Z3,WT, CONTRG,DX,Z2,Z4,ESTL,ESTR,
     $     ESTGL,ESTGR,CONTRB,TSUM,SUMG,DLG,ERRL,EXCESS,BND
      double precision :: ZERO,HALF,ONE,ONEP5,TWO,FOUR,SIX,DXMIN
      INTEGER :: IFLAG, IER, NDF,N, IERC,LVL,IFLT
      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))
      TSUM     =  ESTL  +  ESTR
      SUMG    =  ESTGL +  ESTGR
      DLG     =  ABS(CONTRB - TSUM)
      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
      TSUM     =  TSUM + PSUM(LVL)
      SUMG    =  SUMG + GSUM(LVL)
      LVL     =  LVL - 1
      IF (LVL .GT. 0) GO TO 60
      CONTRB  =  TSUM
      CONTRG  =  SUMG
      LVL     =  1
      DLG     =  ERROR
      GO TO 80
C
C        RESTORE SAVED INFORMATION TO PROCESS RIGHT HALF.
C
   70 PSUM(LVL)  =  TSUM
      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
      DOUBLE PRECISION FUNCTION SDIST(Y,N)
      implicit none
C
C        COMPUTE Y**(N/2 - 1) EXP(-Y) / GAMMA(N/2)
C
C                (Revised: 1994-01-19)
C
      DOUBLE PRECISION :: Y,XN,TEST, ZERO, HALF, ONE, X23,SQRTPI
      INTEGER :: N,JJ,JK,JKP,J
      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         =  DBLE(N) * HALF
      TEST       =  LOG(Y) - Y / DBLE(JKP)
      IF ( TEST .LT. X23 ) THEN
         SDIST = ZERO
         RETURN
      ENDIF
      SDIST = LOG ( SDIST )
      DO 10 J = 1, JKP
         XN    =   XN - ONE
         SDIST   =  SDIST + TEST - LOG(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)
      implicit none
      double precision :: ZERO, ONE, TWO, SMALL, Z, arg, term , f0, g0
     $     ,df,ERB2,HNC,BND,PROB
      INTEGER NN,NDF,N, I, IER,IERC,IFLT
      PARAMETER (NN=100)
      DOUBLE precision :: A,B,H,HL,BPD,D, SDIST
      integer :: INF
      DIMENSION A(NN),B(NN),H(*),HL(*),BPD(*),INF(*),D(*)
      DATA  ZERO, ONE, TWO, SMALL / 0.0, 1.0, 2.0, 1.0E-08 /
      external SDIST
      F0    =  ZERO
      G0    =  ZERO
      IF (Z .LE. -ONE .OR. Z .GE. ONE) RETURN
      DF    =  DBLE(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 * * * * * * * * * * * * * * * * * * * * * * * * * * * *


