SUBROUTINE CHSPPF(P,NU,PPF)
C
C     PURPOSE--THIS SUBROUTINE COMPUTES THE PERCENT POINT
C              FUNCTION VALUE FOR THE CHI-SQUARED DISTRIBUTION
C              WITH INTEGER DEGREES OF FREEDOM PARAMETER = NU.
C              THE CHI-SQUARED DISTRIBUTION USED
C              HEREIN IS DEFINED FOR ALL NON-NEGATIVE X,
C              AND ITS PROBABILITY DENSITY FUNCTION IS GIVEN
C              IN REFERENCES 2, 3, AND 4 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 (INCLUSIVELY) 
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 CHI-SQUARED 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 (INCLUSIVELY)
C                   AND 1.0 (EXCLUSIVELY).
C     OTHER DATAPAC   SUBROUTINES NEEDED--NONE.
C     FORTRAN LIBRARY SUBROUTINES NEEDED--DEXP, DLOG.
C     MODE OF INTERNAL OPERATIONS--DOUBLE PRECISION.
C     LANGUAGE--ANSI FORTRAN. 
C     ACCURACY--(ON THE UNIVAC 1108, EXEC 8 SYSTEM AT NBS)
C               COMPARED TO THE KNOWN NU = 2 (EXPONENTIAL)
C               RESULTS, AGREEMENT WAS HAD OUT TO 6 SIGNIFICANT
C               DIGITS FOR ALL TESTED P IN THE RANGE P = .001 TO
C               P = .999.  FOR P = .95 AND SMALLER, THE AGREEMENT
C               WAS EVEN BETTER--7 SIGNIFICANT DIGITS.
C               (NOTE THAT THE TABULATED VALUES GIVEN IN THE WILK,
C               GNANADESIKAN, AND HUYETT REFERENCE BELOW, PAGE 20,
C               ARE IN ERROR FOR AT LEAST THE GAMMA = 1 CASE--
C               THE WORST DETECTED ERROR WAS AGREEMENT TO ONLY 3
C               SIGNIFICANT DIGITS (IN THEIR 8 SIGNIFICANT DIGIT TABLE)
C               FOR P = .999.)
C     REFERENCES--WILK, GNANADESIKAN, AND HUYETT, 'PROBABILITY
C                 PLOTS FOR THE GAMMA DISTRIBUTION',
C                 TECHNOMETRICS, 1962, PAGES 1-15,
C                 ESPECIALLY PAGES 3-5. 
C               --NATIONAL BUREAU OF STANDARDS APPLIED MATHEMATICS
C                 SERIES 55, 1964, PAGE 257, FORMULA 6.1.41,
C                 AND PAGES 940-943.
C               --JOHNSON AND KOTZ, CONTINUOUS UNIVARIATE
C                 DISTRIBUTIONS--1, 1970, PAGES 166-206.
C               --HASTINGS AND PEACOCK, STATISTICAL
C                 DISTRIBUTIONS--A HANDBOOK FOR
C                 STUDENTS AND PRACTITIONERS, 1975,
C                 PAGES 46-51.
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--SEPTEMBER 1975. 
C     UPDATED         --NOVEMBER  1975. 
C
C---------------------------------------------------------------------
C
      DOUBLE PRECISION DP,DGAMMA
      DOUBLE PRECISION Z,Z2,Z3,Z4,Z5,DEN,A,B,C,D,G
      DOUBLE PRECISION XMIN0,XMIN,AI,XMAX,DX,PCALC,XMID
      DOUBLE PRECISION XLOWER,XUPPER,XDEL
      DOUBLE PRECISION SUM,TERM,CUT1,CUT2,AJ,CUTOFF,T
      DOUBLE PRECISION DEXP,DLOG
      DIMENSION D(10)
      DATA C/ .918938533204672741D0/
      DATA D(1),D(2),D(3),D(4),D(5)
     1                 /+.833333333333333333D-1,-.277777777777777778D-2,
     1+.793650793650793651D-3,-.595238095238095238D-3,+.8417508417508417
     151D-3/
      DATA D(6),D(7),D(8),D(9),D(10)
     1     /-.191752691752691753D-2,+.641025641025641025D-2,-.2955065359
     147712418D-1,+.179644372368830573D0,-.139243221690590111D1/
C
      IPR=6
C
C     CHECK THE INPUT ARGUMENTS FOR ERRORS
C
      IF(P.LT.0.0.OR.P.GE.1.0)GOTO50
      IF(NU.LT.1)GOTO55
      GOTO90
   50 WRITE(IPR,1)
      WRITE(IPR,46)P
      PPF=0.0
      RETURN
   55 WRITE(IPR,15) 
      WRITE(IPR,47)NU
      PPF=0.0
      RETURN
   90 CONTINUE
    1 FORMAT(1H ,115H***** FATAL ERROR--THE FIRST  INPUT ARGUMENT TO THE
     1 CHSPPF SUBROUTINE IS OUTSIDE THE ALLOWABLE (0,1) INTERVAL *****)
   15 FORMAT(1H , 91H***** FATAL ERROR--THE SECOND INPUT ARGUMENT TO THE
     1 CHSPPF 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
C     EXPRESS THE CHI-SQUARED DISTRIBUTION PERCENT POINT
C     FUNCTION IN TERMS OF THE EQUIVALENT GAMMA
C     DISTRIBUTION PERCENT POINT FUNCTION,
C     AND THEN EVALUATE THE LATTER.
C
      ANU=NU
      GAMMA=ANU/2.0 
      DP=P
      DNU=NU
      DGAMMA=DNU/2.0D0
      MAXIT=10000
C
C     COMPUTE THE GAMMA FUNCTION USING THE ALGORITHM IN THE 
C     NBS APPLIED MATHEMATICS SERIES REFERENCE.
C     THIS GAMMA FUNCTION NEED BE CALCULATED ONLY ONCE.
C     IT IS USED IN THE CALCULATION OF THE CDF BASED ON
C     THE TENTATIVE VALUE OF THE PPF IN THE ITERATION.
C
      Z=DGAMMA
      DEN=1.0D0
  150 IF(Z.GE.10.0D0)GOTO160
      DEN=DEN*Z
      Z=Z+1.0D0
      GOTO150
  160 Z2=Z*Z
      Z3=Z*Z2
      Z4=Z2*Z2
      Z5=Z2*Z3
      A=(Z-0.5D0)*DLOG(Z)-Z+C 
      B=D(1)/Z+D(2)/Z3+D(3)/Z5+D(4)/(Z2*Z5)+D(5)/(Z4*Z5)+
     1D(6)/(Z*Z5*Z5)+D(7)/(Z3*Z5*Z5)+D(8)/(Z5*Z5*Z5)+D(9)/(Z2*Z5*Z5*Z5)
      G=DEXP(A+B)/DEN
C
C     DETERMINE LOWER AND UPPER LIMITS ON THE DESIRED 100P
C     PERCENT POINT.
C
      ILOOP=1
      XMIN0=(DP*DGAMMA*G)**(1.0D0/DGAMMA)
      XMIN=XMIN0
      ICOUNT=1
  350 AI=ICOUNT
      XMAX=AI*XMIN0 
      DX=XMAX
      GOTO1000
  360 IF(PCALC.GE.DP)GOTO370
      XMIN=XMAX
      ICOUNT=ICOUNT+1
      IF(ICOUNT.LE.30000)GOTO350
  370 XMID=(XMIN+XMAX)/2.0D0
C
C     NOW ITERATE BY BISECTION UNTIL THE DESIRED ACCURACY IS ACHIEVED.
C
      ILOOP=2
      XLOWER=XMIN
      XUPPER=XMAX
      ICOUNT=0
  550 DX=XMID
      GOTO1000
  560 IF(PCALC.EQ.DP)GOTO570
      IF(PCALC.GT.DP)GOTO580
      XLOWER=XMID
      XMID=(XMID+XUPPER)/2.0D0
      GOTO590
  580 XUPPER=XMID
      XMID=(XMID+XLOWER)/2.0D0
  590 XDEL=XMID-XLOWER
      IF(XDEL.LT.0.0D0)XDEL=-XDEL
      ICOUNT=ICOUNT+1
      IF(XDEL.LT.0.0000000001D0.OR.ICOUNT.GT.100)GOTO570
      GOTO550
  570 PPF=2.0D0*XMID
      RETURN
C
C******************************************************************** 
C     THIS SECTION BELOW IS LOGICALLY SEPARATE FROM THE ABOVE.
C     THIS SECTION COMPUTES A CDF VALUE FOR ANY GIVEN TENTATIVE
C     PERCENT POINT X VALUE AS DEFINED IN EITHER OF THE 2
C     ITERATION LOOPS IN THE ABOVE CODE.
C
C     COMPUTE T-SUB-Q AS DEFINED ON PAGE 4 OF THE WILK, GNANADESIKAN, 
C     AND HUYETT REFERENCE
C
 1000 SUM=1.0D0/DGAMMA
      TERM=1.0D0/DGAMMA
      CUT1=DX-DGAMMA
      CUT2=DX*10000000000.0D0 
      DO700J=1,MAXIT
      AJ=J
      TERM=DX*TERM/(DGAMMA+AJ)
      SUM=SUM+TERM
      CUTOFF=CUT1+(CUT2*TERM/SUM)
      IF(AJ.GT.CUTOFF)GOTO750 
  700 CONTINUE
      WRITE(IPR,705)MAXIT
      WRITE(IPR,706)P
      WRITE(IPR,707)NU
      WRITE(IPR,708)
      PPF=0.0
      RETURN
C
  750 T=SUM
      PCALC=(DX**DGAMMA)*(DEXP(-DX))*T/G
      IF(ILOOP.EQ.1)GOTO360
      GOTO560
C
  705 FORMAT(1H ,48H*****ERROR IN INTERNAL OPERATIONS IN THE CHSPPF , 
     1 45HSUBROUTINE--THE NUMBER OF ITERATIONS EXCEEDS ,I7) 
  706 FORMAT(1H ,33H     THE INPUT VALUE OF P     IS ,E15.8)
  707 FORMAT(1H ,33H     THE INPUT VALUE OF NU    IS ,I8)
  708 FORMAT(1H ,48H     THE OUTPUT VALUE OF PPF HAS BEEN SET TO 0.0) 
C
      END