SUBROUTINE POIPPF(P,ALAMBA,PPF)
C
C     PURPOSE--THIS SUBROUTINE COMPUTES THE PERCENT POINT
C              FUNCTION VALUE AT THE SINGLE PRECISION VALUE P
C              FOR THE POISSON DISTRIBUTION
C              WITH SINGLE PRECISION
C              TAIL LENGTH PARAMETER = ALAMBA.
C              THE POISSON DISTRIBUTION USED
C              HEREIN HAS MEAN = ALAMBA 
C              AND STANDARD DEVIATION = SQRT(ALAMBA).
C              THIS DISTRIBUTION IS DEFINED FOR
C              ALL DISCRETE NON-NEGATIVE INTEGER  X--X = 0, 1, 2, ... .
C              THIS DISTRIBUTION HAS THE PROBABILITY FUNCTION
C              F(X) = EXP(-ALAMBA) * ALAMBA**X / X!.
C              THE POISSON DISTRIBUTION IS THE
C              DISTRIBUTION OF THE NUMBER OF EVENTS
C              IN THE INTERVAL (0,ALAMBA) WHEN
C              THE WAITING TIME BETWEEN EVENTS
C              IS EXPONENTIALLY DISTRIBUTED
C              WITH MEAN = 1 AND STANDARD DEVIATION = 1.
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                     --ALAMBA = THE SINGLE PRECISION VALUE 
C                                OF THE TAIL LENGTH PARAMETER.
C                                ALAMBA SHOULD BE POSITIVE. 
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 POISSON DISTRIBUTION
C             WITH TAIL LENGTH PARAMETER = ALAMBA.
C     PRINTING--NONE UNLESS AN INPUT ARGUMENT ERROR CONDITION EXISTS. 
C     RESTRICTIONS--ALAMBA SHOULD BE POSITIVE.
C                 --P SHOULD BE BETWEEN 0.0 (INCLUSIVELY)
C                   AND 1.0 (EXCLUSIVELY).
C     OTHER DATAPAC   SUBROUTINES NEEDED--NORPPF, POICDF.
C     FORTRAN LIBRARY SUBROUTINES NEEDED--SQRT, DEXP.
C     MODE OF INTERNAL OPERATIONS--SINGLE PRECISION AND DOUBLE PRECISION.
C     LANGUAGE--ANSI FORTRAN. 
C     COMMENT--THE SINGLE PRECISION TAIL LENGTH
C              PARAMETER ALAMBA IS     NOT     RESTRICTED
C              TO ONLY INTEGER VALUES.
C              ALAMBA CAN BE SET TO ANY POSITIVE REAL
C              VALUE--INTEGER OR NON-INTEGER.
C            --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 87-121,
C                 ESPECIALLY PAGE 102, FORMULA 36.1.
C               --HASTINGS AND PEACOCK, STATISTICAL
C                 DISTRIBUTIONS--A HANDBOOK FOR
C                 STUDENTS AND PRACTITIONERS, 1975,
C                 PAGES 108-113.
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 146-154. 
C               --COX AND MILLER, THE THEORY OF STOCHASTIC
C                 PROCESSES, 1965, PAGE 7.
C               --GENERAL ELECTRIC COMPANY, TABLES OF THE
C                 INDIVIDUAL AND CUMULATIVE TERMS OF POISSON
C                 DISTRIBUTION, 1962.
C               --OWEN, HANDBOOK OF STATISTICAL
C                 TABLES, 1962, PAGES 259-261.
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 DLAMBA 
C
      IPR=6
C
C     CHECK THE INPUT ARGUMENTS FOR ERRORS
C
      IF(P.LT.0.0.OR.P.GE.1.0)GOTO50
      IF(ALAMBA.LE.0.0)GOTO55 
      GOTO90
   50 WRITE(IPR,1)
      WRITE(IPR,46)P
      PPF=0.0
      RETURN
   55 WRITE(IPR,15) 
      WRITE(IPR,46)ALAMBA
      PPF=0.0
      RETURN
   90 CONTINUE
    1 FORMAT(1H ,115H***** FATAL ERROR--THE FIRST  INPUT ARGUMENT TO THE
     1 POIPPF SUBROUTINE IS OUTSIDE THE ALLOWABLE (0,1) INTERVAL *****)
   15 FORMAT(1H , 91H***** FATAL ERROR--THE SECOND INPUT ARGUMENT TO THE
     1 POIPPF SUBROUTINE IS NON-POSITIVE *****)
   46 FORMAT(1H , 35H***** THE VALUE OF THE ARGUMENT IS ,E15.8,6H *****)
C
C-----START POINT-----------------------------------------------------
C
      DLAMBA=ALAMBA 
      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) PPF = 0
C
      IF(P.EQ.0.0)GOTO110
      PF0=DEXP(-DLAMBA)
      IF(P.LE.PF0)GOTO110
      GOTO190
  110 PPF=0.0
      RETURN
  190 CONTINUE
C
C     DETERMINE AN INITIAL APPROXIMATION TO THE POISSON
C     PERCENT POINT BY USE OF THE NORMAL APPROXIMATION
C     TO THE POISSON.
C     (SEE JOHNSON AND KOTZ, DISCRETE DISTRIBUTIONS,
C     PAGE 102, FORMULA 36.1).
C
      AMEAN=ALAMBA
      SD=SQRT(ALAMBA)
      CALL NORPPF(P,ZPPF)
      X2=AMEAN-1.0+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 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 POICDF(X2,ALAMBA,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 POICDF(X2,ALAMBA,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 POICDF(X2,ALAMBA,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
CCCCC 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 POISSON PROBABILITIES FOR THE
C     DERIVED LOWER AND UPPER BOUNDS.
C
      X0=IX0
      X1=IX1
      CALL POICDF(X0,ALAMBA,P0)
      CALL POICDF(X1,ALAMBA,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 POICDF(X2,ALAMBA,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)ALAMBA
      RETURN
C
  222 FORMAT(1H ,43HNO UPPER BOUND FOUND AFTER 10**7 ITERATIONS)
  240 FORMAT(1H ,9HIX0    = ,I8,10X,5HP0 = ,F14.7)
  241 FORMAT(1H ,9HIX1    = ,I8,10X,5HP1 = ,F14.7)
  242 FORMAT(1H ,9HIX2    = ,I8,10X,5HP2 = ,F14.7)
  244 FORMAT(1H ,9HP      = ,F14.7)
  245 FORMAT(1H ,9HALAMBA = ,F14.7)
  249 FORMAT(1H ,47H***** INTERNAL ERROR IN POIPPF 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