SUBROUTINE NBPPF(P,PPAR,N,PPF)
C
C     PURPOSE--THIS SUBROUTINE COMPUTES THE PERCENT POINT
C              FUNCTION VALUE AT THE SINGLE PRECISION VALUE P
C              FOR THE NEGATIVE BINOMIAL DISTRIBUTION
C              WITH SINGLE PRECISION 'BERNOULLI PROBABILITY'
C              PARAMETER = PPAR,
C              AND INTEGER 'NUMBER OF SUCCESSES IN BERNOULLI TRIALS'
C              PARAMETER = N. 
C              THE NEGATIVE BINOMIAL DISTRIBUTION USED
C              HEREIN HAS MEAN = N*(1-PPAR)/PPAR
C              AND STANDARD DEVIATION = SQRT(N*(1-PPAR)/(PPAR*PPAR))).
C              THIS DISTRIBUTION IS DEFINED FOR
C              ALL NON-NEGATIVE INTEGER X--X = 0, 1, 2, ... .
C              THIS DISTRIBUTION HAS THE PROBABILITY FUNCTION
C              F(X) = C(N+X-1,N) * PPAR**N * (1-PPAR)**X.
C              WHERE C(N+X-1,N) IS THE COMBINATORIAL FUNCTION
C              EQUALING THE NUMBER OF COMBINATIONS OF N+X-1 ITEMS
C              TAKEN N AT A TIME.
C              THE NEGATIVE BINOMIAL DISTRIBUTION IS THE
C              DISTRIBUTION OF THE NUMBER OF FAILURES
C              BEFORE OBTAINING N SUCCESSES IN AN 
C              INDEFINITE SEQUENCE OF BERNOULLI (0,1)
C              TRIALS WHERE THE PROBABILITY OF SUCCESS
C              IN A SINGLE TRIAL = PPAR.
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                     --PPAR   = THE SINGLE PRECISION VALUE 
C                                OF THE 'BERNOULLI PROBABILITY'
C                                PARAMETER FOR THE NEGATIVE BINOMIAL
C                                DISTRIBUTION.
C                                PPAR SHOULD BE BETWEEN
C                                0.0 (EXCLUSIVELY) AND
C                                1.0 (EXCLUSIVELY).
C                     --N      = THE INTEGER VALUE
C                                OF THE 'NUMBER OF SUCCESSES
C                                IN BERNOULLI TRIALS' PARAMETER.
C                                N SHOULD BE A POSITIVE INTEGER.
C     OUTPUT ARGUMENTS--PPF    = THE SINGLE PRECISION PERCENT
C                                POINT FUNCTION VALUE.
C     OUTPUT--THE SINGLE PRECISION PERCENT POINT  .
C             FUNCTION VALUE PPF
C             FOR THE NEGATIVE BINOMIAL DISTRIBUTION
C             WITH 'BERNOULLI PROBABILITY' PARAMETER = PPAR 
C             AND 'NUMBER OF SUCCESSES IN BERNOULLI TRIALS' 
C             PARAMETER = N.
C     PRINTING--NONE UNLESS AN INPUT ARGUMENT ERROR CONDITION EXISTS. 
C     RESTRICTIONS--PPAR SHOULD BE BETWEEN 0.0 (EXCLUSIVELY)
C                   AND 1.0 (EXCLUSIVELY).
C                 --N SHOULD BE A POSITIVE INTEGER.
C                 --P SHOULD BE BETWEEN 0.0 (INCLUSIVELY)
C                   AND 1.0 (EXCLUSIVELY).
C     OTHER DATAPAC   SUBROUTINES NEEDED--NORPPF, NBCDF.
C     FORTRAN LIBRARY SUBROUTINES NEEDED--SQRT, EXP, ALOG.
C     MODE OF INTERNAL OPERATIONS--SINGLE PRECISION AND DOUBLE PRECISION.
C     LANGUAGE--ANSI FORTRAN. 
C     COMMENT--NOTE THAT EVEN THOUGH THE OUTPUT
C              FROM THIS DISCRETE DISTRIBUTION
C              PERCENT POINT FUNCTION
C              SUBROUTINE MUST NECESSARILY BE A
C              DISCRETE INTEGER VALUE,
C              THE OUTPUT VARIABLE PPF IS SINGLE
C              PRECISION IN MODE.
C              PPF HAS BEEN SPECIFIED AS SINGLE
C              PRECISION SO AS TO CONFORM WITH THE DATAPAC
C              CONVENTION THAT ALL OUTPUT VARIABLES FROM ALL
C              DATAPAC SUBROUTINES ARE SINGLE PRECISION.
C              THIS CONVENTION IS BASED ON THE BELIEF THAT
C              1) A MIXTURE OF MODES (FLOATING POINT
C              VERSUS INTEGER) IS INCONSISTENT AND
C              AN UNNECESSARY COMPLICATION
C              IN A DATA ANALYSIS; AND
C              2) FLOATING POINT MACHINE ARITHMETIC
C              (AS OPPOSED TO INTEGER ARITHMETIC) 
C              IS THE MORE NATURAL MODE FOR DOING 
C              DATA ANALYSIS. 
C     REFERENCES--JOHNSON AND KOTZ, DISCRETE
C                 DISTRIBUTIONS, 1969, PAGES 122-142,
C                 ESPECIALLY PAGE 127, FORMULA 22.
C               --HASTINGS AND PEACOCK, STATISTICAL
C                 DISTRIBUTIONS--A HANDBOOK FOR
C                 STUDENTS AND PRACTITIONERS, 1975,
C                 PAGES 92-95.
C               --NATIONAL BUREAU OF STANDARDS APPLIED MATHEMATICS
C                 SERIES 55, 1964, PAGE 929.
C               --FELLER, AN INTRODUCTION TO PROBABILITY
C                 THEORY AND ITS APPLICATIONS, VOLUME 1,
C                 EDITION 2, 1957, PAGES 155-157, 210.
C               --KENDALL AND STUART, THE ADVANCED THEORY OF
C                 STATISTICS, VOLUME 1, EDITION 2, 1963, PAGES 130-131.
C               --WILLIAMSON AND BRETHERTON, TABLES OF
C                 THE NEGATIVE BINOMIAL PROBABILITY
C                 DISTRIBUTION, 1963.
C               --OWEN, HANDBOOK OF STATISTICAL
C                 TABLES, 1962, PAGE 304.
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--NOVEMBER  1975. 
C
C---------------------------------------------------------------------
C
      DOUBLE PRECISION DPPAR
C
      IPR=6
C
C     CHECK THE INPUT ARGUMENTS FOR ERRORS
C
      IF(P.LT.0.0.OR.P.GE.1.0)GOTO50
      IF(PPAR.LE.0.0.OR.PPAR.GE.1.0)GOTO55
      IF(N.LT.1)GOTO60
      GOTO90
   50 WRITE(IPR,1)
      WRITE(IPR,46)P
      PPF=0.0
      RETURN
   55 WRITE(IPR,11) 
      WRITE(IPR,46)PPAR
      PPF=0.0
      RETURN
   60 WRITE(IPR,25) 
      WRITE(IPR,47)N
      PPF=0.0
      RETURN
   90 CONTINUE
    1 FORMAT(1H ,115H***** FATAL ERROR--THE FIRST  INPUT ARGUMENT TO THE
     1 NBPPF  SUBROUTINE IS OUTSIDE THE ALLOWABLE (0,1) INTERVAL *****)
   11 FORMAT(1H ,115H***** FATAL ERROR--THE SECOND INPUT ARGUMENT TO THE
     1 NBPPF  SUBROUTINE IS OUTSIDE THE ALLOWABLE (0,1) INTERVAL *****)
   25 FORMAT(1H , 91H***** FATAL ERROR--THE THIRD  INPUT ARGUMENT TO THE
     1 NBPPF  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
      AN=N
      DPPAR=PPAR
      PPF=0.0
      IX0=0
      IX1=0
      IX2=0
      P0=0.0
      P1=0.0
      P2=0.0
C
C     TREAT CERTAIN SPECIAL CASES IMMEDIATELY--
C     1) P = 0.0
C     2) P = 0.5 AND PPAR = 0.5
C     3) PPF = 0
C
      IF(P.EQ.0.0)GOTO110
      IF(P.EQ.0.5.AND.PPAR.EQ.0.5)GOTO130
      PF0=DPPAR**N
      IF(P.LE.PF0)GOTO110
      GOTO190
  110 PPF=0.0
      RETURN
  130 PPF=N-1
      RETURN
  190 CONTINUE
C
C     DETERMINE AN INITIAL APPROXIMATION TO THE NEGATIVE BINOMIAL
C     PERCENT POINT BY USE OF THE HYPERBOLIC ARCSIN
C     TRANSFORMATION OF THE NEGATIVE BINOMIAL
C     TO APPROXIMATE NORMALITY.
C     (SEE JOHNSON AND KOTZ, DISCRETE DISTRIBUTIONS,
C     PAGE 127, FORMULA 22).
C
      AMEAN=AN*(1.0-PPAR)/PPAR
      SD=SQRT(AN*(1.0-PPAR)/(PPAR*PPAR))
      ARG=SQRT((AMEAN+0.375)/(AN-0.75)) 
      ARCSH=ALOG(ARG+SQRT(ARG*ARG+1.0)) 
      YMEAN=(SQRT(AN-0.5))*ARCSH
      YSD=0.5
      CALL NORPPF(P,ZPPF)
      YPPF=YMEAN+ZPPF*YSD
      ARG=YPPF/SQRT(AN-0.5)
      E=EXP(ARG)
      SINH=(E-1.0/E)/2.0
      X2=-0.375+(AN-0.75)*SINH*SINH
      X2=X2+0.5
      IX2=X2
C
C     CHECK AND MODIFY (IF NECESSARY) THIS INITIAL
C     ESTIMATE OF THE PERCENT POINT
C     TO ASSURE THAT IT BE NON-NEGATIVE.
C
      IF(IX2.LT.0)IX2=0
C
C     DETERMINE UPPER AND LOWER BOUNDS ON THE DESIRED
C     PERCENT POINT BY ITERATING OUT (BOTH BELOW AND ABOVE) 
C     FROM THE ORIGINAL APPROXIMATION AT STEPS
C     OF 1 STANDARD DEVIATION.
C     THE RESULTING BOUNDS WILL BE AT MOST
C     1 STANDARD DEVIATION APART.
C
      IX0=0
      IX1=10**10
      ISD=SD+1.0
      X2=IX2
      CALL NBCDF(X2,PPAR,N,P2)
C
      IF(P2.LT.P)GOTO210
      GOTO250
C
  210 IX0=IX2
      DO220I=1,100000
      IX2=IX0+ISD
      IF(IX2.GE.IX1)GOTO275
      X2=IX2
      CALL NBCDF(X2,PPAR,N,P2)
      IF(P2.GE.P)GOTO230
      IX0=IX2
  220 CONTINUE
      WRITE(IPR,249)
      WRITE(IPR,222)
      GOTO950
  230 IX1=IX2
      GOTO275
C
  250 IX1=IX2
      DO260I=1,100000
      IX2=IX1-ISD
      IF(IX2.LE.IX0)GOTO275
      X2=IX2
      CALL NBCDF(X2,PPAR,N,P2)
      IF(P2.LT.P)GOTO270
      IX1=IX2
  260 CONTINUE
      WRITE(IPR,249)
      WRITE(IPR,262)
      GOTO950
  270 IX0=IX2
C
  275 IF(IX0.EQ.IX1)GOTO280
      GOTO295
  280 IF(IX0.EQ.0)GOTO285
      IF(IX0.EQ.N)GOTO290
      WRITE(IPR,249)
      WRITE(IPR,282)
      GOTO950
  285 IX1=IX1+1
      GOTO295
  290 IX0=IX0-1
  295 CONTINUE
C
C     COMPUTE NEGATIVE BINOMIAL PROBABILITIES FOR THE
C     DERIVED LOWER AND UPPER BOUNDS.
C
      X0=IX0
      X1=IX1
      CALL NBCDF(X0,PPAR,N,P0)
      CALL NBCDF(X1,PPAR,N,P1)
C
C     CHECK THE PROBABILITIES FOR PROPER ORDERING 
C
      IF(P0.LT.P.AND.P.LE.P1)GOTO490
      IF(P0.EQ.P)GOTO410
      IF(P1.EQ.P)GOTO420
      IF(P0.GT.P1)GOTO430
      IF(P0.GT.P)GOTO440
      IF(P1.LT.P)GOTO450
      WRITE(IPR,249)
      WRITE(IPR,401)
      GOTO950
  410 PPF=IX0
      RETURN
  420 PPF=IX1
      RETURN
  430 WRITE(IPR,249)
      WRITE(IPR,431)
      GOTO950
  440 WRITE(IPR,249)
      WRITE(IPR,441)
      GOTO950
  450 WRITE(IPR,249)
      WRITE(IPR,451)
      GOTO950
  490 CONTINUE
C
C     THE STOPPING CRITERION IS THAT THE LOWER BOUND
C     AND UPPER BOUND ARE EXACTLY 1 UNIT APART.
C     CHECK TO SEE IF IX1 = IX0 + 1;
C     IF SO, THE ITERATIONS ARE COMPLETE;
C     IF NOT, THEN BISECT, COMPUTE PROBABILIIES,
C     CHECK PROBABILITIES, AND CONTINUE ITERATING 
C     UNTIL IX1 = IX0 + 1.
C
  300 IX0P1=IX0+1
      IF(IX1.EQ.IX0P1)GOTO690 
      IX2=(IX0+IX1)/2
      IF(IX2.EQ.IX0)GOTO610
      IF(IX2.EQ.IX1)GOTO620
      X2=IX2
      CALL NBCDF(X2,PPAR,N,P2)
      IF(P0.LT.P2.AND.P2.LT.P1)GOTO630
      IF(P2.LE.P0)GOTO640
      IF(P2.GE.P1)GOTO650
  610 WRITE(IPR,249)
      WRITE(IPR,611)
      GOTO950
  620 WRITE(IPR,249)
      WRITE(IPR,611)
      GOTO950
  630 IF(P2.LE.P)GOTO635
      IX1=IX2
      P1=P2
      GOTO300
  635 IX0=IX2
      P0=P2
      GOTO300
  640 WRITE(IPR,249)
      WRITE(IPR,641)
      GOTO950
  650 WRITE(IPR,249)
      WRITE(IPR,651)
      GOTO950
  690 PPF=IX1
      IF(P0.EQ.P)PPF=IX0
      RETURN
C
  950 WRITE(IPR,240)IX0,P0
      WRITE(IPR,241)IX1,P1
      WRITE(IPR,242)IX2,P2
      WRITE(IPR,244)P
      WRITE(IPR,245)PPAR,N
      RETURN
C
  222 FORMAT(1H ,43HNO UPPER BOUND FOUND AFTER 10**7 ITERATIONS)
  240 FORMAT(1H ,7HIX0  = ,I8,10X,5HP0 = ,F14.7)
  241 FORMAT(1H ,7HIX1  = ,I8,10X,5HP1 = ,F14.7)
  242 FORMAT(1H ,7HIX2  = ,I8,10X,5HP2 = ,F14.7)
  244 FORMAT(1H ,7HP    = ,F14.7)
  245 FORMAT(1H ,7HPPAR = ,F14.7,10X,5HN  = ,I8)
  249 FORMAT(1H ,47H***** INTERNAL ERROR IN NBPPF  SUBROUTINE *****)
  262 FORMAT(1H ,43HNO LOWER BOUND FOUND AFTER 10**7 ITERATIONS)
  282 FORMAT(1H ,31HLOWER AND UPPER BOUND IDENTICAL)
  401 FORMAT(1H ,39HIMPOSSIBLE BRANCH CONDITION ENCOUNTERED)
  431 FORMAT(1H ,42HLOWER BOUND PROBABILITY (P0) GREATER THAN ,
     1 28HUPPER BOUND PROBABILITY (P1)) 
  441 FORMAT(1H ,42HLOWER BOUND PROBABILITY (P0) GREATER THAN ,
     1 21HINPUT PROBABILITY (P))
  451 FORMAT(1H ,42HUPPER BOUND PROBABILITY (P1) LESS    THAN ,
     1 21HINPUT PROBABILITY (P))
  611 FORMAT(1H ,39HBISECTION VALUE (X2) = LOWER BOUND (X0))
  621 FORMAT(1H ,39HBISECTION VALUE (X2) = UPPER BOUND (X1))
  641 FORMAT(1H ,33HBISECTION VALUE PROBABILITY (P2) ,
     1 38HLESS THAN LOWER BOUND PROBABILITY (P0)) 
  651 FORMAT(1H ,33HBISECTION VALUE PROBABILITY (P2) ,
     1 41HGREATER THAN UPPER BOUND PROBABILITY (P1))
C
      END