      SUBROUTINE FCDF(X,NU1,NU2,CDF)
C
C     PURPOSE--THIS SUBROUTINE COMPUTES THE CUMULATIVE DISTRIBUTION
C              FUNCTION VALUE FOR THE F DISTRIBUTION
C              WITH INTEGER DEGREES OF FREEDOM
C              PARAMETERS = NU1 AND NU2.
C              THIS DISTRIBUTION IS DEFINED FOR ALL NON-NEGATIVE X.
C              THE PROBABILITY DENSITY FUNCTION IS GIVEN
C              IN THE REFERENCES BELOW. 
C     INPUT  ARGUMENTS--X      = THE SINGLE PRECISION VALUE AT
C                                WHICH THE CUMULATIVE DISTRIBUTION
C                                FUNCTION IS TO BE EVALUATED.
C                                X SHOULD BE NON-NEGATIVE.
C                     --NU1    = THE INTEGER DEGREES OF FREEDOM
C                                FOR THE NUMERATOR OF THE F RATIO.
C                                NU1 SHOULD BE POSITIVE.
C                     --NU2    = THE INTEGER DEGREES OF FREEDOM
C                                FOR THE DENOMINATOR OF THE F RATIO.
C                                NU2 SHOULD BE POSITIVE.
C     OUTPUT ARGUMENTS--CDF    = THE SINGLE PRECISION CUMULATIVE
C                                DISTRIBUTION FUNCTION VALUE.
C     OUTPUT--THE SINGLE PRECISION CUMULATIVE DISTRIBUTION
C             FUNCTION VALUE CDF FOR THE F DISTRIBUTION
C             WITH DEGREES OF FREEDOM
C             PARAMETERS = NU1 AND NU2. 
C     PRINTING--NONE UNLESS AN INPUT ARGUMENT ERROR CONDITION EXISTS. 
C     RESTRICTIONS--X SHOULD BE NON-NEGATIVE.
C                 --NU1 SHOULD BE A POSITIVE INTEGER VARIABLE.
C                 --NU2 SHOULD BE A POSITIVE INTEGER VARIABLE.
C     OTHER DATAPAC   SUBROUTINES NEEDED--NORCDF,CHSCDF.
C     FORTRAN LIBRARY SUBROUTINES NEEDED--DSQRT, DATAN.
C     MODE OF INTERNAL OPERATIONS--DOUBLE PRECISION.
C     LANGUAGE--ANSI FORTRAN. 
C     REFERENCES--NATIONAL BUREAU OF STANDARDS APPLIED MATHEMATICS
C                 SERIES 55, 1964, PAGES 946-947, 
C                 FORMULAE 26.6.4, 26.6.5, 26.6.8, AND 26.6.15.
C               --JOHNSON AND KOTZ, CONTINUOUS UNIVARIATE
C                 DISTRIBUTIONS--2, 1970, PAGE 83, FORMULA 20,
C                 AND PAGE 84, THIRD FORMULA.
C               --PAULSON, AN APPROXIMATE NORMAILIZATION
C                 OF THE ANALYSIS OF VARIANCE DISTRIBUTION, 
C                 ANNALS OF MATHEMATICAL STATISTICS, 1942,
C                 NUMBER 13, PAGES 233-135.
C               --SCHEFFE AND TUKEY, A FORMULA FOR SAMPLE SIZES
C                 FOR POPULATION TOLERANCE LIMITS, 1944,
C                 NUMBER 15, PAGE 217.
C     WRITTEN BY--JAMES J. FILLIBEN
C                 STATISTICAL ENGINEERING LABORATORY (205.03)
C                 NATIONAL BUREAU OF STANDARDS
C                 WASHINGTON, D. C. 20234
C                 PHONE:  301-921-2315
C     ORIGINAL VERSION--AUGUST    1972. 
C     UPDATED         --SEPTEMBER 1975. 
C     UPDATED         --NOVEMBER  1975. 
C     UPDATED         --OCTOBER   1976. 
C
C---------------------------------------------------------------------
C
      DOUBLE PRECISION DX,PI,ANU1,ANU2,Z,SUM,TERM,AI,COEF1,COEF2,ARG
      DOUBLE PRECISION COEF
      DOUBLE PRECISION THETA,SINTH,COSTH,A,B
      DOUBLE PRECISION DSQRT,DATAN
      DOUBLE PRECISION DFACT1,DFACT2,DNUM,DDEN
      DOUBLE PRECISION DPOW1,DPOW2
      DOUBLE PRECISION DNU1,DNU2
      DOUBLE PRECISION TERM1,TERM2,TERM3
      DATA PI/3.14159265358979D0/
      DATA DPOW1,DPOW2/0.33333333333333D0,0.66666666666667D0/
      DATA NUCUT1,NUCUT2/100,1000/
C
      IPR=6
C
C     CHECK THE INPUT ARGUMENTS FOR ERRORS
C
      IF(NU1.LE.0)GOTO50
      IF(NU2.LE.0)GOTO55
      IF(X.LT.0.0)GOTO60
      GOTO90
   50 WRITE(IPR,15) 
      WRITE(IPR,47)NU1
      CDF=0.0
      RETURN
   55 WRITE(IPR,23) 
      WRITE(IPR,47)NU2
      CDF=0.0
      RETURN
   60 WRITE(IPR,4)
      WRITE(IPR,46)X
      CDF=0.0
      RETURN
   90 CONTINUE
    4 FORMAT(1H , 96H***** NON-FATAL DIAGNOSTIC--THE FIRST  INPUT ARGUME
     1NT TO THE FCDF   SUBROUTINE IS NEGATIVE *****)
   15 FORMAT(1H , 91H***** FATAL ERROR--THE SECOND INPUT ARGUMENT TO THE
     1 FCDF   SUBROUTINE IS NON-POSITIVE *****)
   23 FORMAT(1H , 91H***** FATAL ERROR--THE THIRD  INPUT ARGUMENT TO THE
     1 FCDF   SUBROUTINE IS NON-POSITIVE *****)
   46 FORMAT(1H , 35H***** THE VALUE OF THE ARGUMENT IS ,E15.8,6H *****)
   47 FORMAT(1H , 35H***** THE VALUE OF THE ARGUMENT IS ,I8   ,6H *****)
C
C-----START POINT-----------------------------------------------------
C
      DX=X
      M=NU1
      N=NU2
      ANU1=NU1
      ANU2=NU2
      DNU1=NU1
      DNU2=NU2
C
C     IF X IS NON-POSITIVE, SET CDF = 0.0 AND RETURN.
C     IF NU2 IS 5 THROUGH 9 AND X IS MORE THAN 3000
C     STANDARD DEVIATIONS BELOW THE MEAN,
C     SET CDF = 0.0 AND RETURN.
C     IF NU2 IS 10 OR LARGER AND X IS MORE THAN 150
C     STANDARD DEVIATIONS BELOW THE MEAN,
C     SET CDF = 0.0 AND RETURN.
C     IF NU2 IS 5 THROUGH 9 AND X IS MORE THAN 3000
C     STANDARD DEVIATIONS ABOVE THE MEAN,
C     SET CDF = 1.0 AND RETURN.
C     IF NU2 IS 10 OR LARGER AND X IS MORE THAN 150
C     STANDARD DEVIATIONS ABOVE THE MEAN,
C     SET CDF = 1.0 AND RETURN.
C
      IF(X.LE.0.0)GOTO105
      IF(NU2.LE.4)GOTO109
      T1=2.0/ANU1
      T2=ANU2/(ANU2-2.0)
      T3=(ANU1+ANU2-2.0)/(ANU2-4.0)
      AMEAN=T2
      SD=SQRT(T1*T2*T2*T3)
      ZRATIO=(X-AMEAN)/SD
      IF(NU2.LT.10.AND.ZRATIO.LT.-3000.0)GOTO105
      IF(NU2.GE.10.AND.ZRATIO.LT.-150.0)GOTO105
      IF(NU2.LT.10.AND.ZRATIO.GT.3000.0)GOTO107
      IF(NU2.GE.10.AND.ZRATIO.GT.150.0)GOTO107
      GOTO109
  105 CDF=0.0
      RETURN
  107 CDF=1.0
      RETURN
  109 CONTINUE
C
C     DISTINGUISH BETWEEN 6 SEPARATE REGIONS
C     OF THE (NU1,NU2) SPACE. 
C     BRANCH TO THE PROPER COMPUTATIONAL METHOD
C     DEPENDING ON THE REGION.
C     NUCUT1 HAS THE VALUE 100.
C     NUCUT2 HAS THE VALUE 1000.
C
      IF(NU1.LT.NUCUT2.AND.NU2.LT.NUCUT2)GOTO1000 
      IF(NU1.GE.NUCUT2.AND.NU2.GE.NUCUT2)GOTO2000 
      IF(NU1.LT.NUCUT1.AND.NU2.GE.NUCUT2)GOTO3000 
      IF(NU1.GE.NUCUT1.AND.NU2.GE.NUCUT2)GOTO2000 
      IF(NU1.GE.NUCUT2.AND.NU2.LT.NUCUT1)GOTO5000 
      IF(NU1.GE.NUCUT2.AND.NU2.GE.NUCUT1)GOTO2000 
      IBRAN=5
      WRITE(IPR,99)IBRAN
   99 FORMAT(1H ,42H*****INTERNAL ERROR IN   FCDF SUBROUTINE--,
     146HIMPOSSIBLE BRANCH CONDITION AT BRANCH POINT = ,I8) 
      RETURN
C
C     TREAT THE CASE WHEN NU1 AND NU2
C     ARE BOTH SMALL OR MODERATE
C     (THAT IS, BOTH ARE SMALLER THAN 1000).
C     METHOD UTILIZED--EXACT FINITE SUM 
C     (SEE AMS 55, PAGE 946, FORMULAE 26.6.4, 26.6.5,
C     AND 26.6.8).
C
 1000 CONTINUE
      Z=ANU2/(ANU2+ANU1*DX)
      IFLAG1=NU1-2*(NU1/2)
      IFLAG2=NU2-2*(NU2/2)
      IF(IFLAG1.EQ.0)GOTO120
      IF(IFLAG2.EQ.0)GOTO150
      GOTO250
C
C     DO THE NU1 EVEN AND NU2 EVEN OR ODD CASE
C
  120 SUM=0.0D0
      TERM=1.0D0
      IMAX=(M-2)/2
      IF(IMAX.LE.0)GOTO110
      DO100I=1,IMAX 
      AI=I
      COEF1=2.0D0*(AI-1.0D0)
      COEF2=2.0D0*AI
      TERM=TERM*((ANU2+COEF1)/COEF2)*(1.0D0-Z)
      SUM=SUM+TERM
  100 CONTINUE
C
  110 SUM=SUM+1.0D0 
      SUM=(Z**(ANU2/2.0D0))*SUM
      CDF=1.0D0-SUM 
      RETURN
C
C     DO THE NU1 ODD AND NU2 EVEN CASE
C
  150 SUM=0.0D0
      TERM=1.0D0
      IMAX=(N-2)/2
      IF(IMAX.LE.0)GOTO210
      DO200I=1,IMAX 
      AI=I
      COEF1=2.0D0*(AI-1.0D0)
      COEF2=2.0D0*AI
      TERM=TERM*((ANU1+COEF1)/COEF2)*Z
      SUM=SUM+TERM
  200 CONTINUE
C
  210 SUM=SUM+1.0D0 
      CDF=((1.0D0-Z)**(ANU1/2.0D0))*SUM 
      RETURN
C
C     DO THE NU1 ODD AND NU2 ODD CASE
C
  250 SUM=0.0D0
      TERM=1.0D0
      ARG=DSQRT((ANU1/ANU2)*DX)
      THETA=DATAN(ARG)
      SINTH=ARG/DSQRT(1.0D0+ARG*ARG)
      COSTH=1.0D0/DSQRT(1.0D0+ARG*ARG)
      IF(N.EQ.1)GOTO320
      IF(N.EQ.3)GOTO310
      IMAX=N-2
      DO300I=3,IMAX,2
      AI=I
      COEF1=AI-1.0D0
      COEF2=AI
      TERM=TERM*(COEF1/COEF2)*(COSTH*COSTH)
      SUM=SUM+TERM
  300 CONTINUE
C
  310 SUM=SUM+1.0D0 
      SUM=SUM*SINTH*COSTH
C
  320 A=(2.0D0/PI)*(THETA+SUM)
  350 SUM=0.0D0
      TERM=1.0D0
      IF(M.EQ.1)B=0.0D0
      IF(M.EQ.1)GOTO450
      IF(M.EQ.3)GOTO410
      IMAX=M-3
      DO400I=1,IMAX,2
      AI=I
      COEF1=AI
      COEF2=AI+2.0D0
      TERM=TERM*((ANU2+COEF1)/COEF2)*(SINTH*SINTH)
      SUM=SUM+TERM
  400 CONTINUE
C
  410 SUM=SUM+1.0D0 
      SUM=SUM*SINTH*(COSTH**N)
      COEF=1.0D0
      IEVODD=N-2*(N/2)
      IMIN=3
      IF(IEVODD.EQ.0)IMIN=2
      IF(IMIN.GT.N)GOTO420
      DO430I=IMIN,N,2
      AI=I
      COEF=((AI-1.0D0)/AI)*COEF
  430 CONTINUE
C
  420 COEF=COEF*ANU2
      IF(IEVODD.EQ.0)GOTO440
      COEF=COEF*(2.0D0/PI)
C
  440 B=COEF*SUM
C
  450 CDF=A-B
      RETURN
C
C     TREAT THE CASE WHEN NU1 AND NU2
C     ARE BOTH LARGE
C     (THAT IS, BOTH ARE EQUAL TO OR LARGER THAN 1000);
C     OR WHEN NU1 IS MODERATE AND NU2 IS LARGE
C     (THAT IS, WHEN NU1 IS EQUAL TO OR GREATER THAN 100
C     BUT SMALLER THAN 1000,
C     AND NU2 IS EQUAL TO OR LARGER THAN 1000);
C     OR WHEN NU2 IS MODERATE AND NU1 IS LARGE
C     (THAT IS WHEN NU2 IS EQUAL TO OR GREATER THAN 100
C     BUT SMALLER THAN 1000,
C     AND NU1 IS EQUAL TO OR LARGER THAN 1000).
C     METHOD UTILIZED--PAULSON APPROXIMATION
C     (SEE AMS 55, PAGE 947, FORMULA 26.6.15).
C
 2000 CONTINUE
      DFACT1=1.0D0/(4.5D0*DNU1)
      DFACT2=1.0D0/(4.5D0*DNU2)
      DNUM=((1.0D0-DFACT2)*(DX**DPOW1))-(1.0D0-DFACT1)
      DDEN=DSQRT((DFACT2*(DX**DPOW2))+DFACT1)
      U=DNUM/DDEN
      CALL NORCDF(U,GCDF)
      CDF=GCDF
      RETURN
C
C     TREAT THE CASE WHEN NU1 IS SMALL
C     AND NU2 IS LARGE
C     (THAT IS, WHEN NU1 IS SMALLER THAN 100,
C     AND NU2 IS EQUAL TO OR LARGER THAN 1000).
C     METHOD UTILIZED--SHEFFE-TUKEY APPROXIMATION 
C     (SEE JOHNSON AND KOTZ, VOLUME 2, PAGE 84, THIRD FORMULA).
C
 3000 CONTINUE
      TERM1=DNU1
      TERM2=(DNU1/DNU2)*(0.5D0*DNU1-1.0D0)
      TERM3=-(DNU1/DNU2)*0.5D0
      U=(TERM1+TERM2)/((1.0D0/DX)-TERM3)
      CALL CHSCDF(U,NU1,CCDF) 
      CDF=CCDF
      RETURN
C
C     TREAT THE CASE WHEN NU2 IS SMALL
C     AND NU1 IS LARGE
C     (THAT IS, WHEN NU2 IS SMALLER THAN 100,
C     AND NU1 IS EQUAL TO OR LARGER THAN 1000).
C     METHOD UTILIZED--SHEFFE-TUKEY APPROXIMATION 
C     (SEE JOHNSON AND KOTZ, VOLUME 2, PAGE 84, THIRD FORMULA).
C
 5000 CONTINUE
      TERM1=DNU2
      TERM2=(DNU2/DNU1)*(0.5D0*DNU2-1.0D0)
      TERM3=-(DNU2/DNU1)*0.5D0
      U=(TERM1+TERM2)/(DX-TERM3)
      CALL CHSCDF(U,NU2,CCDF) 
      CDF=1.0-CCDF
      RETURN
C
      END 
