SUBROUTINE BINPPF(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 BINOMIAL DISTRIBUTION
C              WITH SINGLE PRECISION 'BERNOULLI PROBABILITY'
C              PARAMETER = PPAR,
C              AND INTEGER 'NUMBER OF BERNOULLI TRIALS'
C              PARAMETER = N. 
C              THE BINOMIAL DISTRIBUTION USED
C              HEREIN HAS MEAN = N*PPAR 
C              AND STANDARD DEVIATION = SQRT(N*PPAR*(1-PPAR)).
C              THIS DISTRIBUTION IS DEFINED FOR ALL
C              DISCRETE INTEGER X BETWEEN 0 (INCLUSIVELY)
C              AND N (INCLUSIVELY).
C              THIS DISTRIBUTION HAS THE PROBABILITY FUNCTION
C              F(X) = C(N,X) * PPAR**X * (1-PPAR)**(N-X).
C              WHERE C(N,X) IS THE COMBINATORIAL FUNCTION
C              EQUALING THE NUMBER OF COMBINATIONS OF N ITEMS
C              TAKEN X AT A TIME.
C              THE BINOMIAL DISTRIBUTION IS THE
C              DISTRIBUTION OF THE NUMBER OF
C              SUCCESSES IN N 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 (INCLUSIVELY))
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 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 BERNOULLI TRIALS'
C                                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 BINOMIAL DISTRIBUTION
C             WITH 'BERNOULLI PROBABILITY' PARAMETER = PPAR 
C             AND 'NUMBER OF BERNOULLI TRIALS' 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 (INCLUSIVELY).
C     OTHER DATAPAC   SUBROUTINES NEEDED--NORPPF, BINCDF.
C     FORTRAN LIBRARY SUBROUTINES NEEDED--SQRT.
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 50-86,
C                 ESPECIALLY PAGE 64, FORMULA 36. 
C               --HASTINGS AND PEACOCK, STATISTICAL
C                 DISTRIBUTIONS--A HANDBOOK FOR
C                 STUDENTS AND PRACTITIONERS, 1975,
C                 PAGES 36-41.
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 135-142. 
C               --KENDALL AND STUART, THE ADVANCED THEORY OF
C                 STATISTICS, VOLUME 1, EDITION 2, 1963, PAGES 120-125.
C               --MOOD AND GRABLE, INTRODUCTION TO THE THEORY
C                 OF STATISTICS, EDITION 2, 1963, PAGES 64-69.
C               --OWEN, HANDBOOK OF STATISTICAL
C                 TABLES, 1962, PAGES 264-272.
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.GT.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 BINPPF SUBROUTINE IS OUTSIDE THE ALLOWABLE (0,1) INTERVAL *****)
   11 FORMAT(1H ,115H***** FATAL ERROR--THE SECOND INPUT ARGUMENT TO THE
     1 BINPPF SUBROUTINE IS OUTSIDE THE ALLOWABLE (0,1) INTERVAL *****)
   25 FORMAT(1H , 91H***** FATAL ERROR--THE THIRD  INPUT ARGUMENT TO THE
     1 BINPPF 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 OR 1.0
C     2) P = 0.5 AND PPAR = 0.5
C     3) PPF = 0 OR N
C
      IF(P.EQ.0.0)GOTO110
      IF(P.EQ.1.0)GOTO120
      IF(P.EQ.0.5.AND.PPAR.EQ.0.5)GOTO130
      PF0=(1.0D0-DPPAR)**N
      QFN=1.0D0-(DPPAR**N)
      IF(P.LE.PF0)GOTO110
      IF(P.GT.QFN)GOTO120
      GOTO190
  110 PPF=0.0
      RETURN
  120 PPF=N
      RETURN
  130 PPF=N/2
      RETURN
  190 CONTINUE
C
C     DETERMINE AN INITIAL APPROXIMATION TO THE BINOMIAL
C     PERCENT POINT BY USE OF THE NORMAL APPROXIMATION
C     TO THE BINOMIAL.
C     (SEE JOHNSON AND KOTZ, DISCRETE DISTRIBUTIONS,
C     PAGE 64, FORMULA 36).
C
      AMEAN=AN*PPAR 
      SD=SQRT(AN*PPAR*(1.0-PPAR))
      CALL NORPPF(P,ZPPF)
      X2=AMEAN-0.5+ZPPF*SD
      IX2=X2
C
C     CHECK AND MODIFY (IF NECESSARY) THIS INITIAL
C     ESTIMATE OF THE PERCENT POINT
C     TO ASSURE THAT IT BE IN THE CLOSED INTERVAL 0 TO N.
C
      IF(IX2.LT.0)IX2=0
      IF(IX2.GT.N)IX2=N
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=N
      ISD=SD+1.0
      X2=IX2
      CALL BINCDF(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 BINCDF(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 BINCDF(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 BINOMIAL PROBABILITIES FOR THE
C     DERIVED LOWER AND UPPER BOUNDS.
C
      X0=IX0
      X1=IX1
      CALL BINCDF(X0,PPAR,N,P0)
      CALL BINCDF(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 BINCDF(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 BINPPF 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