SUBROUTINE TPPF(P,NU,PPF)
C
C     PURPOSE--THIS SUBROUTINE COMPUTES THE PERCENT POINT
C              FUNCTION VALUE FOR THE STUDENT'S T DISTRIBUTION
C              WITH INTEGER DEGREES OF FREEDOM PARAMETER = NU.
C              THE STUDENT'S T DISTRIBUTION USED
C              HEREIN IS DEFINED FOR ALL X,
C              AND ITS PROBABILITY DENSITY FUNCTION IS GIVEN
C              IN THE REFERENCES BELOW. 
C              NOTE THAT THE PERCENT POINT FUNCTION OF A DISTRIBUTION 
C              IS IDENTICALLY THE SAME AS THE INVERSE CUMULATIVE
C              DISTRIBUTION FUNCTION OF THE DISTRIBUTION.
C     INPUT  ARGUMENTS--P      = THE SINGLE PRECISION VALUE 
C                                (BETWEEN 0.0 (EXCLUSIVELY) 
C                                AND 1.0 (EXCLUSIVELY))
C                                AT WHICH THE PERCENT POINT 
C                                FUNCTION IS TO BE EVALUATED.
C                     --NU     = THE INTEGER NUMBER OF DEGREES
C                                OF FREEDOM.
C                                NU SHOULD BE POSITIVE.
C     OUTPUT ARGUMENTS--PPF    = THE SINGLE PRECISION PERCENT
C                                POINT FUNCTION VALUE.
C     OUTPUT--THE SINGLE PRECISION PERCENT POINT FUNCTION . 
C             VALUE PPF 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                 --P SHOULD BE BETWEEN 0.0 (EXCLUSIVELY)
C                   AND 1.0 (EXCLUSIVELY).
C     OTHER DATAPAC   SUBROUTINES NEEDED--NORPPF. 
C     FORTRAN LIBRARY SUBROUTINES NEEDED--DSIN, DCOS, DSQRT, DATAN.
C     MODE OF INTERNAL OPERATIONS--DOUBLE PRECISION.
C     LANGUAGE--ANSI FORTRAN. 
C     COMMENT--FOR NU = 1 AND NU = 2, THE PERCENT POINT FUNCTION
C              FOR THE T DISTRIBUTION EXISTS IN SIMPLE CLOSED FORM
C              AND SO THE COMPUTED PERCENT POINTS ARE EXACT.
C            --FOR OTHER SMALL VALUES OF NU (NU BETWEEN 3 AND 6,
C              INCLUSIVELY), THE APPROXIMATION
C              OF THE T PERCENT POINT BY THE FORMULA
C              GIVEN IN THE REFERENCE BELOW IS AUGMENTED
C              BY 3 ITERATIONS OF NEWTON'S METHOD FOR
C              ROOT DETERMINATION.
C              THIS IMPROVES THE ACCURACY--ESPECIALLY FOR
C              VALUES OF P NEAR 0 OR 1. 
C     REFERENCES--NATIONAL BUREAU OF STANDARDS APPLIED MATHMATICS
C                 SERIES 55, 1964, PAGE 949, FORMULA 26.7.5.
C               --JOHNSON AND KOTZ, CONTINUOUS UNIVARIATE
C                 DISTRIBUTIONS--2, 1970, PAGE 102,
C                 FORMULA 11. 
C               --FEDERIGHI, 'EXTENDED TABLES OF THE
C                 PERCENTAGE POINTS OF STUDENT'S T
C                 DISTRIBUTION, JOURNAL OF THE
C                 AMERICAN STATISTICAL ASSOCIATION,
C                 1969, PAGES 683-688.
C               --HASTINGS AND PEACOCK, STATISTICAL
C                 DISTRIBUTIONS--A HANDBOOK FOR
C                 STUDENTS AND PRACTITIONERS, 1975,
C                 PAGES 120-123.
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--OCTOBER   1975. 
C     UPDATED         --NOVEMBER  1975. 
C
C---------------------------------------------------------------------
C
      DOUBLE PRECISION PI
      DOUBLE PRECISION SQRT2
      DOUBLE PRECISION DP
      DOUBLE PRECISION DNU
      DOUBLE PRECISION TERM1,TERM2,TERM3,TERM4,TERM5
      DOUBLE PRECISION DPPFN
      DOUBLE PRECISION DPPF,DCON,DARG,Z,S,C
      DOUBLE PRECISION B21
      DOUBLE PRECISION B31,B32,B33,B34
      DOUBLE PRECISION B41,B42,B43,B44,B45
      DOUBLE PRECISION B51,B52,B53,B54,B55,B56
      DOUBLE PRECISION D1,D3,D5,D7,D9
      DATA PI/3.14159265358979D0/
      DATA SQRT2/1.414213562D0/
      DATA B21/0.25D0/
      DATA B31,B32,B33,B34/0.01041666666667D0,5.0D0,16.0D0,3.0D0/
      DATA B41,B42,B43,B44,B45/0.00260416666667D0,3.0D0,19.0D0,17.0D0,
     1                         -15.0D0/ 
      DATA B51,B52,B53,B54,B55,B56/0.00001085069444D0,79.0D0,776.0D0, 
     1                             1482.0D0,-1920.0D0,-945.0D0/
C
      IPR=6
C
C     CHECK THE INPUT ARGUMENTS FOR ERRORS
C
      IF(P.LE.0.0.OR.P.GE.1.0)GOTO50
      GOTO90
   50 WRITE(IPR,1)
      WRITE(IPR,46)P
      RETURN
   90 CONTINUE
    1 FORMAT(1H ,115H***** FATAL ERROR--THE FIRST  INPUT ARGUMENT TO THE
     1 TPPF   SUBROUTINE IS OUTSIDE THE ALLOWABLE (0,1) INTERVAL *****)
   46 FORMAT(1H , 35H***** THE VALUE OF THE ARGUMENT IS ,E15.8,6H *****)
C
C-----START POINT-----------------------------------------------------
C
      DNU=NU
      DP=P
      MAXIT=5
C
      IF(NU.GE.3)GOTO250
      IF(NU.EQ.1)GOTO100
      IF(NU.EQ.2)GOTO200
      WRITE(IPR,105)
  105 FORMAT(1H ,33HINTERNAL ERROR IN TPPF SUBROUTINE)
      PPF=0.0
      RETURN
C
C     TREAT THE NU = 1 (CAUCHY) CASE
C
  100 DARG=PI*DP
      PPF=-DCOS(DARG)/DSIN(DARG)
      RETURN
C
C     TREAT THE NU = 2 CASE
C
  200 TERM1=SQRT2/2.0D0
      TERM2=2.0D0*DP-1.0D0
      TERM3=DSQRT(DP*(1.0D0-DP))
      PPF=TERM1*TERM2/TERM3
      RETURN
C
C     TREAT THE NU GREATER THAN OR EQUAL TO 3 CASE
C
  250 CALL NORPPF(P,PPFN)
      DPPFN=PPFN
      D1=DPPFN
      D3=DPPFN**3
      D5=DPPFN**5
      D7=DPPFN**7
      D9=DPPFN**9
      TERM1=D1
      TERM2=B21*(D3+D1)/DNU
      TERM3=B31*(B32*D5+B33*D3+B34*D1)/(DNU**2)
      TERM4=B41*(B42*D7+B43*D5+B44*D3+B45*D1)/(DNU**3)
      TERM5=B51*(B52*D9+B53*D7+B54*D5+B55*D3+B56*D1)/(DNU**4)
      DPPF=TERM1+TERM2+TERM3+TERM4+TERM5
      PPF=DPPF
      IF(NU.GE.7)RETURN
      IF(NU.EQ.3)GOTO300
      IF(NU.EQ.4)GOTO400
      IF(NU.EQ.5)GOTO500
      IF(NU.EQ.6)GOTO600
      RETURN
C
C     AUGMENT THE RESULTS FOR THE NU = 3 CASE
C
  300 DCON=PI*(DP-0.5D0)
      DARG=DPPF/DSQRT(DNU)
      Z=DATAN(DARG) 
      DO350IPASS=1,MAXIT
      S=DSIN(Z)
      C=DCOS(Z)
      Z=Z-(Z+S*C-DCON)/(2.0D0*C*C)
  350 CONTINUE
      PPF=DSQRT(DNU)*S/C
      RETURN
C
C     AUGMENT THE RESULTS FOR THE NU = 4 CASE
C
  400 DCON=2.0D0*(DP-0.5D0)
      DARG=DPPF/DSQRT(DNU)
      Z=DATAN(DARG) 
      DO450IPASS=1,MAXIT
      S=DSIN(Z)
      C=DCOS(Z)
      Z=Z-((1.0D0+0.5D0*C*C)*S-DCON)/(1.5D0*C*C*C)
  450 CONTINUE
      PPF=DSQRT(DNU)*S/C
      RETURN
C
C     AUGMENT THE RESULTS FOR THE NU = 5 CASE
C
  500 DCON=PI*(DP-0.5D0)
      DARG=DPPF/DSQRT(DNU)
      Z=DATAN(DARG) 
      DO550IPASS=1,MAXIT
      S=DSIN(Z)
      C=DCOS(Z)
      Z=Z-(Z+(C+(2.0D0/3.0D0)*C*C*C)*S-DCON)/((8.0D0/3.0D0)*C**4)
  550 CONTINUE
      PPF=DSQRT(DNU)*S/C
      RETURN
C
C     AUGMENT THE RESULTS FOR THE NU = 6 CASE
C
  600 DCON=2.0D0*(DP-0.5D0)
      DARG=DPPF/DSQRT(DNU)
      Z=DATAN(DARG) 
      DO650IPASS=1,MAXIT
      S=DSIN(Z)
      C=DCOS(Z)
      Z=Z-((1.0D0+0.5D0*C*C+0.375D0*C**4)*S-DCON)/((15.0D0/8.0D0)*C**5)
  650 CONTINUE
      PPF=DSQRT(DNU)*S/C
      RETURN
C
      END