C     ALGORITHM 599, COLLECTED ALGORITHMS FROM ACM.
C     ALGORITHM APPEARED IN ACM-TRANS. MATH. SOFTWARE, VOL.9, NO. 2,
C     JUN., 1983, P. 255-257.
C**********************************************************************
C**********************************************************************
C**********************************************************************
C
C
C
C     F O R T R A N  SOFTWARE PACKAGE FOR RANDOM NUMBER GENERATION
C
C
C
C**********************************************************************
C**********************************************************************
C**********************************************************************
C
C
C
C     CONTENTS:
C
C     1) SUNIF  -  0,1 -UNIFORM DISTRIBUTION
C
C     2) SEXPO  - (STANDARD-) EXPONENTIAL DISTRIBUTION
C
C     3) SNORM  - (STANDARD-) NORMAL DISTRIBUTION
C
C     4) SGAMMA - (STANDARD-) GAMMA DISTRIBUTION
C
C     5) KPOISS - POISSON DISTRIBUTION
C
C
C     THIS PACKAGE CONSTITUTES A FORTRAN-77 DOCUMENTATION OF A SET OF
C     ASSEMBLER FUNCTIONS FOR SAMPLING FROM THE ABOVE DISTRIBUTIONS.
C     ALL ROUTINES MAKE AMPLE USE OF BINARY REPRESENTATIONS OF NUMBERS,
C     THEY ARE AMONG THE MOST ACCURATE AND FAST SAMPLING FUNCTIONS
C     KNOWN. THE FORTRAN PROGRAMS BELOW YIELD THE SAME RANDOM NUMBER
C     SEQUENCES AS THE ONES FROM OUR ASSEMBLER PACKAGE, BUT THEY ARE
C     OF COURSE MUCH SLOWER (BY FACTORS 5-8 ON OUR SIEMENS 7760
C     COMPUTER.)
C     THE SET OF ROUTINES WILL ALSO BE ACCEPTABLE TO FORTRAN IV
C     COMPILERS WHICH ALLOW DATA STATEMENTS FOR ARRAYS WITHOUT
C     IMPLICIT DO-LOOPS.
C
C
C     REMARKS:
C
C     -  NO CARE IS TAKEN TO ENSURE THAT THE PARAMETER VALUES LIE
C        IN THE ALLOWED RANGE (E.G. A/MU > 0.0 FOR SGAMMA/KPOISS).
C
C     -  THE PARAMETER 'IR' MUST BE SET TO SOME  4*K+1 > 0  BEFORE
C        THE FIRST CALL OF ANY OF THE GENERATORS. THEREAFTER IR
C        MUST NOT BE ALTERED UNTIL A NEW INITIALIZATION IS DESIRED.
C
C     -  THE PACKAGE PROVIDES RANDOM DEVIATES OF 6-7 DIGITS ACCURACY.
C        ON MORE ACCURATE COMPUTERS THE CONSTANTS IN SEXPO, SNORM,
C        SGAMMA AND KPOISS OUGHT TO BE ADJUSTED ACCORDING TO LOCAL
C        COMMENTS OR WITH THE AID OF THE TABLES IN THE LITERATURE
C        QUOTED AT THE BEGINNING OF EACH FUNCTION.
C
C
C**********************************************************************C
C**********************************************************************C
C                                                                      C
C                                                                      C
C       0 , 1   - U N I F O R M  DISTRIBUTION                          C
C                                                                      C
C                                                                      C
C**********************************************************************C
C**********************************************************************C
C                                                                      C
C     FOR DETAILS SEE:                                                 C
C                                                                      C
C               AHRENS, J.H., DIETER, U. AND GRUBE, A.                 C
C               PSEUDO-RANDOM NUMBERS:  A NEW PROPOSAL                 C
C                     FOR THE CHOICE OF MULTIPLICATORS                 C
C               COMPUTING, 6 (1970), 121 - 138                         C
C                                                                      C
C**********************************************************************C
C
      REAL FUNCTION SUNIF(IR)

c*********************************************************************72
c
cc SUNIF samples the standard uniform distribution on [0,1].
c
c  Modified:
c
c    22 April 2013
c
c  Author:
c
c    Joachim Ahrens, Klaus-Dieter Kohrt, Ulrich Dieter
c
c  Reference:
c
c    Joachim Ahrens, Klaus-Dieter Kohrt, Ulrich Dieter,
c    Algorithm 599:
c    Sampling from Gamma and Poisson Distributions,
c    ACM Transactions on Mathematical Software,
c    Volume 9, Number 2, June 1983, pages 255-257.
c
c  Parameters:
c
      DOUBLE PRECISION R,FACTOR,TWO28
      SAVE R
C
C     FACTOR - INTEGER OF THE FORM 8*K+5 AS CLOSE AS POSSIBLE
C              TO  2**26 * (SQRT(5)-1)/2     (GOLDEN SECTION)
C     TWO28  = 2**28  (I.E. 28 SIGNIFICANT BITS FOR DEVIATES)
C
      DATA FACTOR /41475557.0D0/, TWO28 /268435456.0D0/
C
C     RETURNS SAMPLE U FROM THE  0,1 -UNIFORM DISTRIBUTION
C     BY A MULTIPLICATIVE CONGRUENTIAL GENERATOR OF THE FORM
C        R := R * FACTOR (MOD 1) .
C     IN THE FIRST CALL R IS INITIALIZED TO
C        R := IR / 2**28 ,
C     WHERE IR MUST BE OF THE FORM  IR = 4*K+1.
C     THEN R ASSUMES ALL VALUES  0 < (4*K+1)/2**28 < 1 DURING
C     A FULL PERIOD 2**26 OF SUNIF.
C     THE PARAMETER IR IS USED ONLY IN THE FIRST CALL FOR
C     INITIALIZATION OF SUNIF. THEREAFTER (WHEN NEGATIVE)
C     IR BECOMES A DUMMY VARIABLE.
C
      IF (IR .GE. 0) GO TO 1
C
C     STANDARD CASE:  SAMPLING
C
      R=DMOD(R*FACTOR,1.0D0)
      SUNIF=SNGL(R)
      RETURN
C
C     FIRST CALL: INITIALIZATION
C
1     R=DBLE(FLOAT(IR))/TWO28
      R=DMOD(R*FACTOR,1.0D0)
      SUNIF=SNGL(R)
      IR=-1
      RETURN
      END
C
C**********************************************************************C
C**********************************************************************C
C                                                                      C
C                                                                      C
C     (STANDARD-)  E X P O N E N T I A L   DISTRIBUTION                C
C                                                                      C
C                                                                      C
C**********************************************************************C
C**********************************************************************C
C                                                                      C
C     FOR DETAILS SEE:                                                 C
C                                                                      C
C               AHRENS, J.H. AND DIETER, U.                            C
C               COMPUTER METHODS FOR SAMPLING FROM THE                 C
C               EXPONENTIAL AND NORMAL DISTRIBUTIONS.                  C
C               COMM. ACM, 15,10 (OCT. 1972), 873 - 882.               C
C                                                                      C
C     ALL STATEMENT NUMBERS CORRESPOND TO THE STEPS OF ALGORITHM       C
C     'SA' IN THE ABOVE PAPER (SLIGHTLY MODIFIED IMPLEMENTATION)       C
C                                                                      C
C**********************************************************************C
C
      REAL FUNCTION SEXPO(IR)

c*********************************************************************72
c
cc SEXPO samples the standard exponential distribution.
c
c  Modified:
c
c    22 April 2013
c
c  Author:
c
c    Joachim Ahrens, Klaus-Dieter Kohrt, Ulrich Dieter
c
c  Reference:
c
c    Joachim Ahrens, Klaus-Dieter Kohrt, Ulrich Dieter,
c    Algorithm 599:
c    Sampling from Gamma and Poisson Distributions,
c    ACM Transactions on Mathematical Software,
c    Volume 9, Number 2, June 1983, pages 255-257.
c
c  Parameters:
c
      DIMENSION Q(8)
      EQUIVALENCE (Q(1),Q1)
C
C     Q(N) = SUM(ALOG(2.0)**K/K])    K=1,..,N ,      THE HIGHEST N
C     (HERE 8) IS DETERMINED BY Q(N)=1.0 WITHIN STANDARD PRECISION
C
      DATA Q/.6931472,.9333737,.9888778,.9984959,
     ,.9998293,.9999833,.9999986,.9999999/
C
      A=0.0
      U=SUNIF(IR)
      GO TO 2
   3  A=A+Q1
   2  U=U+U
      IF (U.LE.1.0) GO TO 3
      U=U-1.0
      IF (U.GT.Q1) GO TO 6
      SEXPO=A+U
      RETURN
   6  I=1
      USTAR=SUNIF(IR)
      UMIN=USTAR
   7  USTAR=SUNIF(IR)
      IF (USTAR.LT.UMIN) UMIN=USTAR
      I=I+1
      IF (U.GT.Q(I)) GO TO 7
      SEXPO=A+UMIN*Q1
      RETURN
      END
C
C**********************************************************************C
C**********************************************************************C
C                                                                      C
C                                                                      C
C     (STANDARD-)  N O R M A L  DISTRIBUTION                           C
C                                                                      C
C                                                                      C
C**********************************************************************C
C**********************************************************************C
C                                                                      C
C     FOR DETAILS SEE:                                                 C
C                                                                      C
C               AHRENS, J.H. AND DIETER, U.                            C
C               EXTENSIONS OF FORSYTHE'S METHOD FOR RANDOM             C
C               SAMPLING FROM THE NORMAL DISTRIBUTION.                 C
C               MATH. COMPUT., 27,124 (OCT. 1973), 927 - 937.          C
C                                                                      C
C     ALL STATEMENT NUMBERS CORRESPOND TO THE STEPS OF ALGORITHM 'FL'  C
C     (M=5) IN THE ABOVE PAPER     (SLIGHTLY MODIFIED IMPLEMENTATION)  C
C                                                                      C
C**********************************************************************C
C
      REAL FUNCTION SNORM(IR)

c*********************************************************************72
c
cc SNORM samples the standard normal distribution.
c
c  Modified:
c
c    22 April 2013
c
c  Author:
c
c    Joachim Ahrens, Klaus-Dieter Kohrt, Ulrich Dieter
c
c  Reference:
c
c    Joachim Ahrens, Klaus-Dieter Kohrt, Ulrich Dieter,
c    Algorithm 599:
c    Sampling from Gamma and Poisson Distributions,
c    ACM Transactions on Mathematical Software,
c    Volume 9, Number 2, June 1983, pages 255-257.
c
c  Parameters:
c
      DIMENSION A(32),D(31),T(31),H(31)
C
C     THE DEFINITIONS OF THE CONSTANTS A(K), D(K), T(K) AND
C     H(K) ARE ACCORDING TO THE ABOVEMENTIONED ARTICLE
C
      DATA A/0.0,.3917609E-1,.7841241E-1,.1177699,.1573107,
     ,.1970991,.2372021,.2776904,.3186394,.3601299,.4022501,
     ,.4450965,.4887764,.5334097,.5791322,.6260990,.6744898,
     ,.7245144,.7764218,.8305109,.8871466,.9467818,1.009990,
     ,1.077516,1.150349,1.229859,1.318011,1.417797,1.534121,
     ,1.675940,1.862732,2.153875/
      DATA D/5*0.0,.2636843,.2425085,.2255674,.2116342,.1999243,
     ,.1899108,.1812252,.1736014,.1668419,.1607967,.1553497,
     ,.1504094,.1459026,.1417700,.1379632,.1344418,.1311722,
     ,.1281260,.1252791,.1226109,.1201036,.1177417,.1155119,
     ,.1134023,.1114027,.1095039/
      DATA T/.7673828E-3,.2306870E-2,.3860618E-2,.5438454E-2,
     ,.7050699E-2,.8708396E-2,.1042357E-1,.1220953E-1,.1408125E-1,
     ,.1605579E-1,.1815290E-1,.2039573E-1,.2281177E-1,.2543407E-1,
     ,.2830296E-1,.3146822E-1,.3499233E-1,.3895483E-1,.4345878E-1,
     ,.4864035E-1,.5468334E-1,.6184222E-1,.7047983E-1,.8113195E-1,
     ,.9462444E-1,.1123001,.1364980,.1716886,.2276241,.3304980,
     ,.5847031/
      DATA H/.3920617E-1,.3932705E-1,.3950999E-1,.3975703E-1,
     ,.4007093E-1,.4045533E-1,.4091481E-1,.4145507E-1,.4208311E-1,
     ,.4280748E-1,.4363863E-1,.4458932E-1,.4567523E-1,.4691571E-1,
     ,.4833487E-1,.4996298E-1,.5183859E-1,.5401138E-1,.5654656E-1,
     ,.5953130E-1,.6308489E-1,.6737503E-1,.7264544E-1,.7926471E-1,
     ,.8781922E-1,.9930398E-1,.1155599,.1404344,.1836142,.2790016,
     ,.7010474/
C
      U=SUNIF(IR)
      S=0.0
      IF (U.GE.0.5) S=1.0
      U=U+U-S
      U=32.0*U
      I=INT(U)
      IF (I.EQ.0) GO TO 9
C
C  START CENTER
C
      USTAR=U-FLOAT(I)
      AA=A(I)
   4  IF (USTAR.LE.T(I)) GO TO 5
      W=(USTAR-T(I))*H(I)
C
C  EXIT   (BOTH CASES)
C
  17  Y=AA+W
      SNORM=Y
      IF (S.EQ.1.0) SNORM=-Y
      RETURN
C
C  CENTER CONTINUED
C
   5  U=SUNIF(IR)
      W=U*(A(I+1)-AA)
      TT=(0.5*W+AA)*W
      GO TO 6
   8  TT=U
      USTAR=SUNIF(IR)
   6  IF (USTAR.GT.TT) GO TO 17
      U=SUNIF(IR)
      IF (USTAR.GE.U) GO TO 8
      USTAR=SUNIF(IR)
      GO TO 4
C
C  START TAIL
C
   9  I=6
      AA=A(32)
      GO TO 10
  11  AA=AA+D(I)
      I=I+1
  10  U=U+U
      IF (U.LT.1.0) GO TO 11
      U=U-1.0
  13  W=U*D(I)
      TT=(0.5*W+AA)*W
      GO TO 14
  16  TT=U
  14  USTAR=SUNIF(IR)
      IF (USTAR.GT.TT) GO TO 17
      U=SUNIF(IR)
      IF (USTAR.GE.U) GO TO 16
      U=SUNIF(IR)
      GO TO 13
      END
C
C**********************************************************************C
C**********************************************************************C
C                                                                      C
C                                                                      C
C     (STANDARD-)  G A M M A  DISTRIBUTION                             C
C                                                                      C
C                                                                      C
C**********************************************************************C
C**********************************************************************C
C                                                                      C
C               PARAMETER  A >= 1.0  ]                                 C
C                                                                      C
C**********************************************************************C
C                                                                      C
C     FOR DETAILS SEE:                                                 C
C                                                                      C
C               AHRENS, J.H. AND DIETER, U.                            C
C               GENERATING GAMMA VARIATES BY A                         C
C               MODIFIED REJECTION TECHNIQUE.                          C
C               COMM. ACM, 25,1 (JAN. 1982), 47 - 54.                  C
C                                                                      C
C     STEP NUMBERS CORRESPOND TO ALGORITHM 'GD' IN THE ABOVE PAPER     C
C                                 (STRAIGHTFORWARD IMPLEMENTATION)     C
C                                                                      C
C**********************************************************************C
C                                                                      C
C               PARAMETER  0.0 < A < 1.0  ]                            C
C                                                                      C
C**********************************************************************C
C                                                                      C
C     FOR DETAILS SEE:                                                 C
C                                                                      C
C               AHRENS, J.H. AND DIETER, U.                            C
C               COMPUTER METHODS FOR SAMPLING FROM GAMMA,              C
C               BETA, POISSON AND BINOMIAL DISTRIBUTIONS.              C
C               COMPUTING, 12 (1974), 223 - 246.                       C
C                                                                      C
C     (ADAPTED IMPLEMENTATION OF ALGORITHM 'GS' IN THE ABOVE PAPER)    C
C                                                                      C
C**********************************************************************C
C
      REAL FUNCTION SGAMMA(IR,A)

c*********************************************************************72
c
cc SGAMMA samples the standard gamma distribution.
c
c  Modified:
c
c    22 April 2013
c
c  Author:
c
c    Joachim Ahrens, Klaus-Dieter Kohrt, Ulrich Dieter
c
c  Reference:
c
c    Joachim Ahrens, Klaus-Dieter Kohrt, Ulrich Dieter,
c    Algorithm 599:
c    Sampling from Gamma and Poisson Distributions,
c    ACM Transactions on Mathematical Software,
c    Volume 9, Number 2, June 1983, pages 255-257.
c
c  Parameters:
c
C     INPUT:  IR=CURRENT STATE OF BASIC RANDOM NUMBER GENERATOR
C             A =PARAMETER (MEAN) OF THE STANDARD GAMMA DISTRIBUTION
C     OUTPUT: SGAMMA = SAMPLE FROM THE GAMMA-(A)-DISTRIBUTION
C
C     COEFFICIENTS Q(K) - FOR Q0 = SUM(Q(K)*A**(-K))
C     COEFFICIENTS A(K) - FOR Q = Q0+(T*T/2)*SUM(A(K)*V**K)
C     COEFFICIENTS E(K) - FOR EXP(Q)-1 = SUM(E(K)*Q**K)
C
      DATA Q1,Q2,Q3,Q4,Q5,Q6,Q7 /.04166669,.02083148,
     ,.00801191,.00144121,-.00007388,.00024511,.00024240/
      DATA A1,A2,A3,A4,A5,A6,A7 /.3333333,-.2500030,
     ,.2000062,-.1662921,.1423657,-.1367177,.1233795/
      DATA E1,E2,E3,E4,E5 /1.,.4999897,.1668290,.0407753,.0102930/
C
C     PREVIOUS A PRE-SET TO ZERO - AA IS A', AAA IS A"
C     SQRT32 IS THE SQUAREROOT OF 32 = 5.656854249492380
C
      DATA AA /0.0/, AAA /0.0/, SQRT32 /5.656854/
C
      IF (A .EQ. AA) GO TO 1
      IF (A .LT. 1.0) GO TO 12
C
C     STEP  1:  RECALCULATIONS OF S2,S,D IF A HAS CHANGED
C
      AA=A
      S2=A-0.5
      S=SQRT(S2)
      D=SQRT32-12.0*S
C
C     STEP  2:  T=STANDARD NORMAL DEVIATE,
C               X=(S,1/2)-NORMAL DEVIATE.
C               IMMEDIATE ACCEPTANCE (I)
C
   1  T=SNORM(IR)
      X=S+0.5*T
      SGAMMA=X*X
      IF (T .GE. 0.0) RETURN
C
C     STEP  3:  U= 0,1 -UNIFORM SAMPLE. SQUEEZE ACCEPTANCE (S)
C
      U=SUNIF(IR)
      IF (D*U .LE. T*T*T) RETURN
C
C     STEP  4:  RECALCULATIONS OF Q0,B,SI,C IF NECESSARY
C
      IF (A .EQ. AAA) GO TO 4
      AAA=A
      R=1.0/A
      Q0=((((((Q7*R+Q6)*R+Q5)*R+Q4)*R+Q3)*R+Q2)*R+Q1)*R
C
C               APPROXIMATION DEPENDING ON SIZE OF PARAMETER A
C               THE CONSTANTS IN THE EXPRESSIONS FOR B, SI AND
C               C WERE ESTABLISHED BY NUMERICAL EXPERIMENTS
C
      IF (A .LE. 3.686) GO TO 3
      IF (A .LE. 13.022) GO TO 2
C
C               CASE 3:  A .GT. 13.022
C
      B=1.77
      SI=.75
      C=.1515/S
      GO TO 4
C
C               CASE 2:  3.686 .LT. A .LE. 13.022
C
   2  B=1.654+.0076*S2
      SI=1.68/S+.275
      C=.062/S+.024
      GO TO 4
C
C               CASE 1:  A .LE. 3.686
C
   3  B=.463+S-.178*S2
      SI=1.235
      C=.195/S-.079+.016*S
C
C     STEP  5:  NO QUOTIENT TEST IF X NOT POSITIVE
C
   4  IF (X .LE. 0.0) GO TO 7
C
C     STEP  6:  CALCULATION OF V AND QUOTIENT Q
C
      V=T/(S+S)
      IF (ABS(V) .LE. 0.25) GO TO 5
      Q=Q0-S*T+0.25*T*T+(S2+S2)*ALOG(1.0+V)
      GO TO 6
   5  Q=Q0+0.5*T*T*((((((A7*V+A6)*V+A5)*V+A4)*V+A3)*V+A2)*V+A1)*V
C
C     STEP  7:  QUOTIENT ACCEPTANCE (Q)
C
   6  IF (ALOG(1.0-U) .LE. Q) RETURN
C
C     STEP  8:  E=STANDARD EXPONENTIAL DEVIATE
C               U= 0,1 -UNIFORM DEVIATE
C               T=(B,SI)-DOUBLE EXPONENTIAL (LAPLACE) SAMPLE
C
   7  E=SEXPO(IR)
      U=SUNIF(IR)
      U=U+U-1.0
      T=B+SIGN(SI*E,U)
C
C     STEP  9:  REJECTION IF T .LT. TAU(1) = -.71874483771719
C
      IF (T .LT. (-.7187449)) GO TO 7
C
C     STEP 10:  CALCULATION OF V AND QUOTIENT Q
C
      V=T/(S+S)
      IF (ABS(V) .LE. 0.25) GO TO 8
      Q=Q0-S*T+0.25*T*T+(S2+S2)*ALOG(1.0+V)
      GO TO 9
   8  Q=Q0+0.5*T*T*((((((A7*V+A6)*V+A5)*V+A4)*V+A3)*V+A2)*V+A1)*V
C
C     STEP 11:  HAT ACCEPTANCE (H) (IF Q NOT POSITIVE GO TO STEP 8)
C
   9  IF (Q .LE. 0.0) GO TO 7
      IF (Q .LE. 0.5) GO TO 10
      W=EXP(Q)-1.0
      GO TO 11
  10  W=((((E5*Q+E4)*Q+E3)*Q+E2)*Q+E1)*Q
C
C               IF T IS REJECTED, SAMPLE AGAIN AT STEP 8
C
  11  IF (C*ABS(U) .GT. W*EXP(E-0.5*T*T)) GO TO 7
      X=S+0.5*T
      SGAMMA=X*X
      RETURN
C
C     ALTERNATE METHOD FOR PARAMETERS A BELOW 1  (.3678794=EXP(-1.))
C
  12  AA=0.0
      B=1.0+.3678794*A
  13  P=B*SUNIF(IR)
      IF (P .GE. 1.0) GO TO 14
      SGAMMA=EXP(ALOG(P)/A)
      IF (SEXPO(IR) .LT. SGAMMA) GO TO 13
      RETURN
  14  SGAMMA=-ALOG((B-P)/A)
      IF (SEXPO(IR) .LT. (1.0-A)*ALOG(SGAMMA)) GO TO 13
      RETURN
      END
C
C**********************************************************************C
C**********************************************************************C
C                                                                      C
C                                                                      C
C     P O I S S O N  DISTRIBUTION                                      C
C                                                                      C
C                                                                      C
C**********************************************************************C
C**********************************************************************C
C                                                                      C
C     FOR DETAILS SEE:                                                 C
C                                                                      C
C               AHRENS, J.H. AND DIETER, U.                            C
C               COMPUTER GENERATION OF POISSON DEVIATES                C
C               FROM MODIFIED NORMAL DISTRIBUTIONS.                    C
C               ACM TRANS. MATH. SOFTWARE, 8,2 (JUNE 1982), 163 - 179. C
C                                                                      C
C     (SLIGHTLY MODIFIED VERSION OF THE PROGRAM IN THE ABOVE ARTICLE)  C
C                                                                      C
C**********************************************************************C
C
      INTEGER FUNCTION KPOISS(IR,MU)

c*********************************************************************72
c
cc KPOISSON samples the poisson distribution.
c
c  Modified:
c
c    22 April 2013
c
c  Author:
c
c    Joachim Ahrens, Klaus-Dieter Kohrt, Ulrich Dieter
c
c  Reference:
c
c    Joachim Ahrens, Klaus-Dieter Kohrt, Ulrich Dieter,
c    Algorithm 599:
c    Sampling from Gamma and Poisson Distributions,
c    ACM Transactions on Mathematical Software,
c    Volume 9, Number 2, June 1983, pages 255-257.
c
c  Parameters:
c
C     INPUT:  IR=CURRENT STATE OF BASIC RANDOM NUMBER GENERATOR
C             MU=MEAN MU OF THE POISSON DISTRIBUTION
C     OUTPUT: KPOISS=SAMPLE FROM THE POISSON-(MU)-DISTRIBUTION
C
      REAL MU, MUPREV, MUOLD
C
C     MUPREV=PREVIOUS MU, MUOLD=MU AT LAST EXECUTION OF STEP P OR B.
C     TABLES: COEFFICIENTS A0-A7 FOR STEP F. FACTORIALS FACT
C     COEFFICIENTS A(K) - FOR PX = FK*V*V*SUM(A(K)*V**K)-DEL
C
      DIMENSION FACT(10), PP(35)
      DATA MUPREV,MUOLD /0.,0./
      DATA A0,A1,A2,A3,A4,A5,A6,A7 /-.5,.3333333,-.2500068,
     ,.2000118,-.1661269,.1421878,-.1384794,.1250060/
      DATA FACT /1.,1.,2.,6.,24.,120.,720.,5040.,40320.,362880./
C
C     SEPARATION OF CASES A AND B
C
      IF (MU .EQ. MUPREV) GO TO 1
      IF (MU .LT. 10.0) GO TO 12
C
C     C A S E  A. (RECALCULATION OF S,D,L IF MU HAS CHANGED)
C
      MUPREV=MU
      S=SQRT(MU)
      D=6.0*MU*MU
C
C             THE POISSON PROBABILITIES PK EXCEED THE DISCRETE NORMAL
C             PROBABILITIES FK WHENEVER K >= M(MU). L=IFIX(MU-1.1484)
C             IS AN UPPER BOUND TO M(MU) FOR ALL MU >= 10 .
C
      L=IFIX(MU-1.1484)
C
C     STEP N. NORMAL SAMPLE - SNORM(IR) FOR STANDARD NORMAL DEVIATE
C
   1  G=MU+S*SNORM(IR)
      IF (G .LT. 0.0) GO TO 2
      KPOISS=IFIX(G)
C
C     STEP I. IMMEDIATE ACCEPTANCE IF KPOISS IS LARGE ENOUGH
C
      IF (KPOISS .GE. L) RETURN
C
C     STEP S. SQUEEZE ACCEPTANCE - SUNIF(IR) FOR (0,1)-SAMPLE U
C
      FK=FLOAT(KPOISS)
      DIFMUK=MU-FK
      U=SUNIF(IR)
      IF (D*U .GE. DIFMUK*DIFMUK*DIFMUK) RETURN
C
C     STEP P. PREPARATIONS FOR STEPS Q AND H.
C             (RECALCULATIONS OF PARAMETERS IF NECESSARY)
C             .3989423=(2*PI)**(-.5)  .416667E-1=1./24.  .1428571=1./7.
C             THE QUANTITIES B1, B2, C3, C2, C1, C0 ARE FOR THE HERMITE
C             APPROXIMATIONS TO THE DISCRETE NORMAL PROBABILITIES FK.
C             C=.1069/MU GUARANTEES MAJORIZATION BY THE 'HAT'-FUNCTION.
C
   2  IF (MU .EQ. MUOLD) GO TO 3
      MUOLD=MU
      OMEGA=.3989423/S
      B1=.4166667E-1/MU
      B2=.3*B1*B1
      C3=.1428571*B1*B2
      C2=B2-15.*C3
      C1=B1-6.*B2+45.*C3
      C0=1.-B1+3.*B2-15.*C3
      C=.1069/MU
   3  IF (G .LT. 0.0) GO TO 5
C
C             'SUBROUTINE' F IS CALLED (KFLAG=0 FOR CORRECT RETURN)
C
      KFLAG=0
      GO TO 7
C
C     STEP Q. QUOTIENT ACCEPTANCE (RARE CASE)
C
   4  IF (FY-U*FY .LE. PY*EXP(PX-FX)) RETURN
C
C     STEP E. EXPONENTIAL SAMPLE - SEXPO(IR) FOR STANDARD EXPONENTIAL
C             DEVIATE E AND SAMPLE T FROM THE LAPLACE 'HAT'
C             (IF T <= -.6744 THEN PK < FK FOR ALL MU >= 10.)
C
   5  E=SEXPO(IR)
      U=SUNIF(IR)
      U=U+U-1.0
      T=1.8+SIGN(E,U)
      IF (T .LE. (-.6744)) GO TO 5
      KPOISS=IFIX(MU+S*T)
      FK=FLOAT(KPOISS)
      DIFMUK=MU-FK
C
C             'SUBROUTINE' F IS CALLED (KFLAG=1 FOR CORRECT RETURN)
C
      KFLAG=1
      GO TO 7
C
C     STEP H. HAT ACCEPTANCE (E IS REPEATED ON REJECTION)
C
   6  IF (C*ABS(U) .GT. PY*EXP(PX+E)-FY*EXP(FX+E)) GO TO 5
      RETURN
C
C     STEP F. 'SUBROUTINE' F. CALCULATION OF PX,PY,FX,FY.
C             CASE KPOISS .LT. 10 USES FACTORIALS FROM TABLE FACT
C
   7  IF (KPOISS .GE. 10) GO TO 8
      PX=-MU
      PY=MU**KPOISS/FACT(KPOISS+1)
      GO TO 11
C
C             CASE KPOISS .GE. 10 USES POLYNOMIAL APPROXIMATION
C             A0-A7 FOR ACCURACY WHEN ADVISABLE
C             .8333333E-1=1./12.  .3989423=(2*PI)**(-.5)
C
   8  DEL=.8333333E-1/FK
      DEL=DEL-4.8*DEL*DEL*DEL
      V=DIFMUK/FK
      IF (ABS(V) .LE. 0.25) GO TO 9
      PX=FK*ALOG(1.0+V)-DIFMUK-DEL
      GO TO 10
   9  PX=FK*V*V*(((((((A7*V+A6)*V+A5)*V+A4)*V+A3)*V+A2)*V+A1)*V+A0)-DEL
  10  PY=.3989423/SQRT(FK)
  11  X=(0.5-DIFMUK)/S
      XX=X*X
      FX=-0.5*XX
      FY=OMEGA*(((C3*XX+C2)*XX+C1)*XX+C0)
      IF (KFLAG) 4,4,6
C
C     C A S E  B. (START NEW TABLE AND CALCULATE P0 IF NECESSARY)
C
  12  MUPREV=0.0
      IF (MU .EQ. MUOLD) GO TO 13
      MUOLD=MU
      M=MAX0(1,IFIX(MU))
      L=0
      P=EXP(-MU)
      Q=P
      P0=P
C
C     STEP U. UNIFORM SAMPLE FOR INVERSION METHOD
C
  13  U=SUNIF(IR)
      KPOISS=0
      IF (U .LE. P0) RETURN
C
C     STEP T. TABLE COMPARISON UNTIL THE END PP(L) OF THE
C             PP-TABLE OF CUMULATIVE POISSON PROBABILITIES
C             (0.458=PP(9) FOR MU=10)
C
      IF (L .EQ. 0) GO TO 15
      J=1
      IF (U .GT. 0.458) J=MIN0(L,M)
      DO 14 K=J,L
      IF (U .LE. PP(K)) GO TO 18
  14  CONTINUE
      IF (L .EQ. 35) GO TO 13
C
C     STEP C. CREATION OF NEW POISSON PROBABILITIES P
C             AND THEIR CUMULATIVES Q=PP(K)
C
  15  L=L+1
      DO 16 K=L,35
      P=P*MU/FLOAT(K)
      Q=Q+P
      PP(K)=Q
      IF (U .LE. Q) GO TO 17
  16  CONTINUE
      L=35
      GO TO 13
  17  L=K
  18  KPOISS=K
      RETURN
      END