* CDFNML
c
c     Latest revision  -  03/15/90  (JRD)
c
c     This routine computes the normal cumulative distribution
c     function from the error function as described in chapter 13
c     of distributions in statistics - continuous univariate
c     distributions - 1, by Johnson and Kotz.
c
c     Written by - Janet R. Donaldson
c                  Statistical Engineering Division
c                  National Bureau of Standards, Boulder, Colorado
c
c     Creation date - December 7, 1981
c
c-----------------------------------------------------------------------
      REAL FUNCTION CDFNML(X)
c  variable declarations
      IMPLICIT NONE
c  scalar arguments
      REAL X

c  external functions
      REAL ERF
      EXTERNAL ERF

c     variable definitions (alphabetically)
c
c     REAL X 
c        the percent point from the chi squared distribution.
c
c     commence body of routine
c
      CDFNML = 0.5E0 * (1.0E0 +  ERF(X/SQRT(2.0E0)))
C
      RETURN
      END


* CDFT
c
c     Latest revision  -  03/15/90  (JRD)
c
c     Purpose--This subroutine computes the cumulative distribution
c              function value for Student's t distribution
c              with integer degrees of freedom parameter = idf.
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 value at
c                                which the cumulative distribution
c                                function is to be evaluated.
c                                X should be non-negative.
c                     --IDF     = the integer number of degrees
c                                of freedom.
c                                IDF 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 = IDF.
c     Printing--None unless an input argument error condition exists.
c     Restrictions--IDF should be a positive integer variable.
c     Other datapac subroutines needed--NORCDFC [? is now CDFNML (AJA, 2016)]
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     Original version--June      1972.
c     Updated         --May       1974.
c     Updated         --September 1975.
c     Updated         --November  1975.
c     Updated         --October   1976.
c
c---------------------------------------------------------------------
      REAL FUNCTION CDFT (X, IDF)

c  Variable declarations
      IMPLICIT NONE

c  Scalar arguments
      REAL X
      INTEGER IDF

c  Local scalars
      REAL
     $   B11,B21,B22,B23,B24,B25,B31,B32,B33,B34,B35,B36,B37,C,CSQ,D1,
     $   D11,D3,D5,D7,D9,DCONST,DF,FPSPM,PI,SD,SUM,TERM,TERM1,TERM2,
     $   TERM3,Z
      INTEGER I,IDFCUT,IEVODD,IMAX,IMIN

c  External functions
      REAL CDFNML, RMACHINE

c     Variable definitions (alphabetically)
C
C     REAL B11, B21, B22, B23, B24, B25
c        constants used in the computations.
C     REAL B31, B32, B33, B34, B35, B36, B37
c        constants used in the computations.
C     REAL DF        the degrees of freedom.
C     REAL D1, D11, D3, D5, D7, D9
c        constants used in the computations.
C     REAL FPSPM       the floating point smallest positive magnitude.
C     INTEGER I        an index.
C     INTEGER IDF        the degrees of freedom.
C     REAL X        the t statistic.

      DATA IDFCUT /1000/
      DATA DCONST /0.3989422804E0/
      DATA B11 /0.25E0/
      DATA B21 /96.0E0/
C     DATA B21 /0.01041666666667E0/
      DATA B22, B23, B24, B25 /3.0E0,-7.0E0,-5.0E0,-3.0E0/
      DATA B31 /0.00260416666667E0/
      DATA B32, B33, B34, B35, B36, B37
     $    /1.0E0,-11.0E0,14.0E0,6.0E0,-3.0E0,-15.0E0/

      PI = ACOS(-1.)
      FPSPM = rmachine(2)

c     Check the input arguments for errors.
      IF (IDF.LE.0) GO TO 10
      GO TO 20
   10 PRINT 1000
      PRINT 1010, IDF
      CDFT = 0.0E0
      RETURN
   20 CONTINUE

c-----start point-----------------------------------------------------

      DF = IDF

c     If IDF is 3 through 9 and X is more than 3000
c     standard deviations below the mean,
c     set CDFT = 0.0e0 and return.
c     If IDF is 10 or larger and x is more than 150
c     standard deviations below the mean,
c     set CDFT = 0.0e0 and return.
c     If IDF is 3 through 9 and x is more than 3000
c     standard deviations above the mean,
c     set CDFT = 1.0e0 and return.
c     If IDF is 10 or larger and x is more than 150
c     standard deviations above the mean,
c     set CDFT = 1.0e0 and return.

      IF (IDF .LE. 2) GO TO 50
      SD = SQRT(DF/(DF-2.0E0))
      Z = X/SD
      IF (IDF.LT.10 .AND. Z.LT.(-3000.0E0)) GO TO 30
      IF (IDF.GE.10 .AND. Z.LT.(-150.0E0)) GO TO 30
      IF (IDF.LT.10 .AND. Z.GT.3000.0E0) GO TO 40
      IF (IDF.GE.10 .AND. Z.GT.150.0E0) GO TO 40
      GO TO 50
   30 CDFT = 0.0E0
      RETURN
   40 CDFT = 1.0E0
      RETURN
   50 CONTINUE

c     Distinguish between the small and moderate
c     degrees of freedom case versus the
c     large degrees of freedom case

      IF (IDF.LT.IDFCUT) GO TO 60
      GO TO 120

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).

   60 CONTINUE
      CSQ = DF/(X*X+DF)
      C = SQRT(CSQ)
      IMAX = IDF - 2
      IEVODD = IDF - 2*(IDF/2)
      IF (IEVODD.NE.0) THEN
         IF (IDF.EQ.1) THEN
            SUM = 0.0E0
         ELSE
            SUM = C
         END IF
         TERM = C
         IMIN = 3
      ELSE
         SUM = 1.0E0
         TERM = 1.0E0
         IMIN = 2
      END IF

      DO 90 I=IMIN,IMAX,2
         IF (TERM.NE.0.0E0) THEN
            IF (LOG(TERM)+LOG((I-1.0E0)/I)+LOG(CSQ).GE.LOG(FPSPM)) THEN
               TERM = TERM*((I-1.0E0)/I)*CSQ
               SUM = SUM + TERM
            ELSE
               TERM = 0.0E0
            END IF
         END IF
   90 CONTINUE

      IF (SUM.EQ.0.0E0 .OR. X.EQ.0.0E0) THEN
         SUM = 0.0E0
      ELSE
         IF (LOG(SUM)+LOG(ABS(X))-0.5*LOG(X*X+DF) .LT. LOG(FPSPM)) THEN
            SUM = 0.0E0
         ELSE
            SUM = SUM*X/SQRT(X*X+DF)
         END IF
      END IF
      IF (IEVODD.EQ.0) GO TO 110
      SUM = (2.0E0/PI)*(ATAN(X/SQRT(DF))+SUM)
  110 CDFT = 0.5E0 + SUM/2.0E0
      RETURN

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).

  120 CONTINUE
      D1 = X
      D3 = X**3
      D5 = X**5
      D7 = X**7
      D9 = X**9
      D11 = X**11
      TERM1 = B11*(D3+D1)/DF
C     TERM2 = B21*(B22*D7+B23*D5+B24*D3+B25*D1)/(DF**2)
      TERM2 = (B22*D7+B23*D5+B24*D3+B25*D1)/(DF**2) / B21
      TERM3 = B31*(B32*D11+B33*D9+B34*D7+B35*D5+B36*D3+B37*D1)/(DF**3)
      CDFT = TERM1 + TERM2 + TERM3
      CDFT = CDFNML(X) - (DCONST*(EXP(-X*X/2.0E0)))*CDFT
      RETURN

C     Format statements

 1000 FORMAT (' ', '***** FATAL ERROR--THE SECOND INPUT ARGUMENT TO ',
     $   'THE CDFT SUBROUTINE IS NON-POSITIVE *****')
 1010 FORMAT (' ', '***** THE VALUE OF THE ARGUMENT IS ', I8, ' *****')
      END
