SUBROUTINE TCDF(X,NU,CDF)
C
C     PURPOSE--THIS SUBROUTINE COMPUTES THE CUMULATIVE DISTRIBUTION
C              FUNCTION VALUE FOR STUDENT'S T DISTRIBUTION
C              WITH INTEGER DEGREES OF FREEDOM PARAMETER = NU.
C              THIS DISTRIBUTION IS DEFINED FOR ALL 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                     --NU     = THE INTEGER NUMBER OF DEGREES
C                                OF FREEDOM.
C                                NU 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 STUDENT'S T DISTRIBUTION
C             WITH DEGREES OF FREEDOM PARAMETER = NU.
C     PRINTING--NONE UNLESS AN INPUT ARGUMENT ERROR CONDITION EXISTS. 
C     RESTRICTIONS--NU SHOULD BE A POSITIVE INTEGER VARIABLE.
C     OTHER DATAPAC   SUBROUTINES NEEDED--NORCDF. 
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 MATHMATICS
C                 SERIES 55, 1964, PAGE 948, FORMULAE 26.7.3 AND 26.7.4.
C               --JOHNSON AND KOTZ, CONTINUOUS UNIVARIATE
C                 DISTRIBUTIONS--2, 1970, PAGES 94-129.
C               --FEDERIGHI, EXTENDED TABLES OF THE
C                 PERCENTAGE POINTS OF STUDENT'S
C                 T-DISTRIBUTION, JOURNAL OF THE
C                 AMERICAN STATISTICAL ASSOCIATION,
C                 1959, PAGES 683-688.
C               --OWEN, HANDBOOK OF STATISTICAL TABLES,
C                 1962, PAGES 27-30.
C               --PEARSON AND HARTLEY, BIOMETRIKA TABLES
C                 FOR STATISTICIANS, VOLUME 1, 1954,
C                 PAGES 132-134.
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--JUNE      1972. 
C     UPDATED         --MAY       1974. 
C     UPDATED         --SEPTEMBER 1975. 
C     UPDATED         --NOVEMBER  1975. 
C     UPDATED         --OCTOBER   1976. 
C
C---------------------------------------------------------------------
C
      DOUBLE PRECISION DX,DNU,PI,C,CSQ,S,SUM,TERM,AI
      DOUBLE PRECISION DSQRT,DATAN
      DOUBLE PRECISION DCONST 
      DOUBLE PRECISION TERM1,TERM2,TERM3
      DOUBLE PRECISION DCDFN
      DOUBLE PRECISION DCDF
      DOUBLE PRECISION B11
      DOUBLE PRECISION B21,B22,B23,B24,B25
      DOUBLE PRECISION B31,B32,B33,B34,B35,B36,B37
      DOUBLE PRECISION D1,D3,D5,D7,D9,D11
      DATA NUCUT/1000/
      DATA PI/3.14159265358979D0/
      DATA DCONST/0.3989422804D0/
      DATA B11/0.25D0/
      DATA B21/0.01041666666667D0/
      DATA B22,B23,B24,B25/3.0D0,-7.0D0,-5.0D0,-3.0D0/
      DATA B31/0.00260416666667D0/
      DATA B32,B33,B34,B35,B36,B37/1.0D0,-11.0D0,14.0D0,6.0D0,
     1                            -3.0D0,-15.0D0/ 
C
      IPR=6
C
C     CHECK THE INPUT ARGUMENTS FOR ERRORS
C
      IF(NU.LE.0)GOTO50
      GOTO90
   50 WRITE(IPR,15) 
      WRITE(IPR,47)NU
      CDF=0.0
      RETURN
   90 CONTINUE
   15 FORMAT(1H , 91H***** FATAL ERROR--THE SECOND INPUT ARGUMENT TO THE
     1 TCDF   SUBROUTINE IS NON-POSITIVE *****)
   47 FORMAT(1H , 35H***** THE VALUE OF THE ARGUMENT IS ,I8   ,6H *****)
C
C-----START POINT-----------------------------------------------------
C
      DX=X
      ANU=NU
      DNU=NU
C
C     IF NU IS 3 THROUGH 9 AND X IS MORE THAN 3000
C     STANDARD DEVIATIONS BELOW THE MEAN,
C     SET CDF = 0.0 AND RETURN.
C     IF NU 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 NU IS 3 THROUGH 9 AND X IS MORE THAN 3000
C     STANDARD DEVIATIONS ABOVE THE MEAN,
C     SET CDF = 1.0 AND RETURN.
C     IF NU 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(NU.LE.2)GOTO109
      SD=SQRT(ANU/(ANU-2.0))
      Z=X/SD
      IF(NU.LT.10.AND.Z.LT.-3000.0)GOTO107
      IF(NU.GE.10.AND.Z.LT.-150.0)GOTO107
      IF(NU.LT.10.AND.Z.GT.3000.0)GOTO108
      IF(NU.GE.10.AND.Z.GT.150.0)GOTO108
      GOTO109
  107 CDF=0.0
      RETURN
  108 CDF=1.0
      RETURN
  109 CONTINUE
C
C     DISTINGUISH BETWEEN THE SMALL AND MODERATE
C     DEGREES OF FREEDOM CASE VERSUS THE
C     LARGE DEGREES OF FREEDOM CASE
C
      IF(NU.LT.NUCUT)GOTO110
      GOTO250
C
C     TREAT THE SMALL AND MODERATE DEGREES OF FREEDOM CASE
C     METHOD UTILIZED--EXACT FINITE SUM 
C     (SEE AMS 55, PAGE 948, FORMULAE 26.7.3 AND 26.7.4).
C
  110 CONTINUE
      C=DSQRT(DNU/(DX*DX+DNU))
      CSQ=DNU/(DX*DX+DNU)
      S=DX/DSQRT(DX*DX+DNU)
      IMAX=NU-2
      IEVODD=NU-2*(NU/2)
      IF(IEVODD.EQ.0)GOTO120
C
      SUM=C
      IF(NU.EQ.1)SUM=0.0D0
      TERM=C
      IMIN=3
      GOTO130
C
  120 SUM=1.0D0
      TERM=1.0D0
      IMIN=2
C
  130 IF(IMIN.GT.IMAX)GOTO160 
      DO100I=IMIN,IMAX,2
      AI=I
      TERM=TERM*((AI-1.0D0)/AI)*CSQ
      SUM=SUM+TERM
  100 CONTINUE
C
  160 SUM=SUM*S
      IF(IEVODD.EQ.0)GOTO170
      SUM=(2.0D0/PI)*(DATAN(DX/DSQRT(DNU))+SUM)
  170 CDF=0.5D0+SUM/2.0D0
      RETURN
C
C     TREAT THE LARGE DEGREES OF FREEDOM CASE.
C     METHOD UTILIZED--TRUNCATED ASYMPTOTIC EXPANSION
C     (SEE JOHNSON AND KOTZ, VOLUME 2, PAGE 102, FORMULA 10;
C     SEE FEDERIGHI, PAGE 687).
C
  250 CONTINUE
      CALL NORCDF(X,CDFN)
      DCDFN=CDFN
      D1=DX
      D3=DX**3
      D5=DX**5
      D7=DX**7
      D9=DX**9
      D11=DX**11
      TERM1=B11*(D3+D1)/DNU
      TERM2=B21*(B22*D7+B23*D5+B24*D3+B25*D1)/(DNU**2)
      TERM3=B31*(B32*D11+B33*D9+B34*D7+B35*D5+B36*D3+B37*D1)/(DNU**3) 
      DCDF=TERM1+TERM2+TERM3
      DCDF=DCDFN-(DCONST*(DEXP(-DX*DX/2.0D0)))*DCDF
      CDF=DCDF
      RETURN
C
      END