      SUBROUTINE PROBF(F,V1,V2,NWDS,P,THETA,Z,B1) 
C 
C     * F. ZWIERS  -- JAN 18/85 
C 
C     * EXACT CALCULATION OF PROB(F>EF) WHERE A HAS V1 AND V2 DEGREES OF
C     * FREEDOM.
C     * 
C     * ARGUEMENT LIST: 
C     *  INPUT -
C     *    F    - ARRAY OF F VALUES 
C     * 
C     *           ***************************************************** 
C     *           *                                                   * 
C     *           *  THE CONTENTS OF ARRAY F ARE MODIFIED BY THIS     * 
C     *           *  SUBROUTINE                                       * 
C     *           *                                                   * 
C     *           ***************************************************** 
C     * 
C     *    V1   - DEGREES OF FREEDOM FOR THE NUMERATOR OF THE F-RATIOS
C     *           (SCALAR)
C     *    V2   - DEGREES OF FREEDOM FOR THE DENOMINATOR OF THE F-
C     *           RATIOS (SCALAR) 
C     *    NWDS - LENGTH OF ARRAY F 
C     * 
C     *  OUTPUT - 
C     *    P    - ARRAY OF PROBABILITIES OF OBSERVING F-RATIOS WITH V1
C     *           AND V2 DEGREES OF FREEDOM WHICH ARE LARGER THAN THE 
C     *           CORRESPONDING F-VALUES CONTAINED IN ARRAY F.
C     * 
C     *  WORK SPACE - 
C     *    THETA,Z,B1 - WORK ARRAYS OF LENGTH NWDS
C     * 
C     * 
C     * REF - LACKRITZ, J.R., 1984: EXACT P VALUES FOR F AND T TESTS, 
C     *       THE AMERICAN STATISTICIAN, 38, 312-314. 
C 
      IMPLICIT REAL (A-H,O-Z),
     +INTEGER (I-N)
      REAL F(1),P(1),THETA(1),Z(1),B1(1)
      INTEGER V1,V2 
      LOGICAL N1EVEN,N2EVEN 
      DATA PI/3.1415927E0/,EPS/1.E-6/ 
C 
C-----------------------------------------------------------------------------
C 
      IF(V1.LT.1 .OR. V2.LT.1)THEN
         WRITE(6,6010) V1,V2
         CALL XIT('PROBF',-1) 
      ENDIF 
C 
C 
      DO 100 I=1,NWDS 
         Z(I)=MAX(EPS,F(I)) 
  100 CONTINUE
      DO 105 I=1,NWDS 
         F(I)=Z(I)
  105 CONTINUE
C 
C 
      N1=V1 
      N2=V2 
      IF(V2.LT.V1)THEN
         DO 110 I=1,NWDS
           F(I)=1.0E0/F(I)
  110    CONTINUE 
         N1=V2
         N2=V1
      ENDIF 
C 
C 
      FN1=FLOAT(N1) 
      FN2=FLOAT(N2) 
      M1=N1/2 
      M2=N2/2 
      N1EVEN=.TRUE.
      IF(M1*2.LT.N1) N1EVEN=.FALSE. 
      N2EVEN=.TRUE.
      IF(M2*2.LT.N2) N2EVEN=.FALSE. 
      FN1BN2=FN1/FN2
      DO 120 I=1,NWDS 
        THETA(I)=F(I)*FN1BN2
  120 CONTINUE
C 
C     * V1 EVEN 
C 
      IF(N1EVEN)THEN
         DO 130 I=1,NWDS
            F(I) = THETA(I)/(1.0E0 + THETA(I))
            Z(I)  = 1.0E0 
            P(I)  = 1.0E0 
  130    CONTINUE 
         IF(M1.GT.1)THEN
            DO 200 I=1,M1-1 
               FTWOI=FLOAT(I+I) 
               FACTI=(FN2 + FTWOI - 2.0E0)/FTWOI
               DO 140 J=1,NWDS
                  Z(J)=Z(J)*F(J)*FACTI
                  P(J)=P(J) + Z(J)
  140          CONTINUE 
  200       CONTINUE
         ENDIF
         POWER=FN2/2.0E0
         DO 250 I=1,NWDS
            P(I)=P(I)/((1.0E0 + THETA(I))**POWER) 
  250    CONTINUE 
C 
C     * V1 ODD, V2 EVEN 
C 
      ELSEIF(N2EVEN)THEN
         DO 260 I=1,NWDS
            THETA(I)=1.0E0/THETA(I) 
            Z(I)=1.0E0
            P(I)=1.0E0
            F(I)=THETA(I)/(1.0E0 + THETA(I))
  260    CONTINUE 
C 
         IF(M2.GT.1)THEN
            DO 300 I=1,M2-1 
               FTWOI=FLOAT(I+I) 
               FACTI=(FN1 + FTWOI - 2.0E0)/FTWOI
               DO 270 J=1,NWDS
                  Z(J)=Z(J)*F(J)*FACTI
                  P(J)=P(J) + Z(J)
  270          CONTINUE 
  300       CONTINUE
         ENDIF
C 
         POWER=FN1/2.0E0
         DO 350 I=1,NWDS
            P(I)=1.0E0 - P(I)/((1.0E0 + THETA(I))**POWER) 
  350    CONTINUE 
C 
C     * V1 ODD, V2 ODD
C 
      ELSE
C 
         IF(N1.EQ.1 .AND. N2.EQ.1)THEN
            DO 360 I=1,NWDS 
               P(I)=2.0E0*ATAN(1.0E0/SQRT(THETA(I)))/PI 
  360       CONTINUE
         ELSE 
            DO 370 I=1,NWDS 
               F(I)=2.0E0/(1.0E0 + THETA(I))
               Z(I)=F(I)
               P(I)=F(I)
  370       CONTINUE
            IF(M2.GT.1)THEN 
               DO 400 I=2,M2
                  FTWOI=FLOAT(I+I)
                  FACTI=FLOAT(I-1)/(FTWOI-1.0E0)
                  DO 380 J=1,NWDS 
                     Z(J)=Z(J)*F(J)*FACTI 
                     P(J)=P(J) + Z(J) 
  380             CONTINUE
  400          CONTINUE 
            ENDIF 
            DO 410 I=1,NWDS 
               B1(I)=SQRT(THETA(I))*P(I)/PI 
  410       CONTINUE
            IF(N1.EQ.1)THEN 
               DO 420 I=1,NWDS
                  P(I)=2.0E0*ATAN(1.0E0/SQRT(THETA(I)))/PI - B1(I)
  420          CONTINUE 
            ELSE
               FACT=1.0E0 
               DO 430 I=1,M2-1
                  FACT=FACT*FLOAT(I)/(FLOAT(I)+0.5E0) 
  430          CONTINUE 
               FACT=2.0E0*FACT
               DO 440 I=1,NWDS
                  F(I)=2.0E0*THETA(I)/(1.0E0 + THETA(I))
                  Z(I)=1.0E0
                  P(I)=0.0E0
  440          CONTINUE 
               DO 500 I=1,M1
                  FACTI=FLOAT(M2 + I - 1)/FLOAT(I + I - 1)
                  DO 450 J=1,NWDS 
                     Z(J)=Z(J)*F(J)*FACTI 
                     P(J)=P(J) + Z(J) 
  450             CONTINUE
  500          CONTINUE 
               DO 510 I=1,NWDS
                  P(I)=P(I)*FACT/ 
     1                   (PI*SQRT(THETA(I))*(1.0E0 + THETA(I))**M2) 
                         P(I)=P(I) - B1(I) + 2.0E0*ATAN(1.0E0/
     &                     SQRT(THETA(I)))/PI 
  510          CONTINUE 
            ENDIF 
         ENDIF
      ENDIF 
C 
C 
C 
      IF(V2.LT.V1) THEN 
         DO 520 I=1,NWDS
            P(I)=1.0E0 - P(I) 
  520    CONTINUE 
      ENDIF 
      RETURN
C 
C-----------------------------------------------------------------------------
C 
 6010 FORMAT('0PROBF CALLED WITH DEGREES OF FREEDOM ',2I5)
      END 
