1      REAL FUNCTION sgamma(a)
2C**********************************************************************C
3C                                                                      C
4C                                                                      C
5C     (STANDARD-)  G A M M A  DISTRIBUTION                             C
6C                                                                      C
7C                                                                      C
8C**********************************************************************C
9C**********************************************************************C
10C                                                                      C
11C               PARAMETER  A >= 1.0  !                                 C
12C                                                                      C
13C**********************************************************************C
14C                                                                      C
15C     FOR DETAILS SEE:                                                 C
16C                                                                      C
17C               AHRENS, J.H. AND DIETER, U.                            C
18C               GENERATING GAMMA VARIATES BY A                         C
19C               MODIFIED REJECTION TECHNIQUE.                          C
20C               COMM. ACM, 25,1 (JAN. 1982), 47 - 54.                  C
21C                                                                      C
22C     STEP NUMBERS CORRESPOND TO ALGORITHM 'GD' IN THE ABOVE PAPER     C
23C                                 (STRAIGHTFORWARD IMPLEMENTATION)     C
24C                                                                      C
25C     Modified by Barry W. Brown, Feb 3, 1988 to use RANF instead of   C
26C     SUNIF.  The argument IR thus goes away.                          C
27C                                                                      C
28C**********************************************************************C
29C                                                                      C
30C               PARAMETER  0.0 < A < 1.0  !                            C
31C                                                                      C
32C**********************************************************************C
33C                                                                      C
34C     FOR DETAILS SEE:                                                 C
35C                                                                      C
36C               AHRENS, J.H. AND DIETER, U.                            C
37C               COMPUTER METHODS FOR SAMPLING FROM GAMMA,              C
38C               BETA, POISSON AND BINOMIAL DISTRIBUTIONS.              C
39C               COMPUTING, 12 (1974), 223 - 246.                       C
40C                                                                      C
41C     (ADAPTED IMPLEMENTATION OF ALGORITHM 'GS' IN THE ABOVE PAPER)    C
42C                                                                      C
43C**********************************************************************C
44C
45C
46C     INPUT: A =PARAMETER (MEAN) OF THE STANDARD GAMMA DISTRIBUTION
47C     OUTPUT: SGAMMA = SAMPLE FROM THE GAMMA-(A)-DISTRIBUTION
48C
49C     COEFFICIENTS Q(K) - FOR Q0 = SUM(Q(K)*A**(-K))
50C     COEFFICIENTS A(K) - FOR Q = Q0+(T*T/2)*SUM(A(K)*V**K)
51C     COEFFICIENTS E(K) - FOR EXP(Q)-1 = SUM(E(K)*Q**K)
52C
53C     .. Scalar Arguments ..
54      REAL a
55C     ..
56C     .. Local Scalars .. (JJV added B0 to fix rare and subtle bug)
57      REAL a1,a2,a3,a4,a5,a6,a7,aa,aaa,b,b0,c,d,e,e1,e2,e3,e4,e5,p,q,q0,
58     +     q1,q2,q3,q4,q5,q6,q7,r,s,s2,si,sqrt32,t,u,v,w,x
59C     ..
60C     .. External Functions ..
61      REAL ranf,sexpo,snorm
62      EXTERNAL ranf,sexpo,snorm
63C     ..
64C     .. Intrinsic Functions ..
65      INTRINSIC abs,alog,exp,sign,sqrt
66C     ..
67C     .. Save statement ..
68C     JJV added Save statement for vars in Data satatements
69      SAVE aa,aaa,s2,s,d,q0,b,si,c,q1,q2,q3,q4,q5,q6,q7,a1,a2,a3,a4,a5,
70     +     a6,a7,e1,e2,e3,e4,e5,sqrt32
71C     ..
72C     .. Data statements ..
73C
74C     PREVIOUS A PRE-SET TO ZERO - AA IS A', AAA IS A"
75C     SQRT32 IS THE SQUAREROOT OF 32 = 5.656854249492380
76C
77      DATA q1,q2,q3,q4,q5,q6,q7/.04166669,.02083148,.00801191,.00144121,
78     +     -.00007388,.00024511,.00024240/
79      DATA a1,a2,a3,a4,a5,a6,a7/.3333333,-.2500030,.2000062,-.1662921,
80     +     .1423657,-.1367177,.1233795/
81      DATA e1,e2,e3,e4,e5/1.,.4999897,.1668290,.0407753,.0102930/
82      DATA aa/0.0/,aaa/0.0/,sqrt32/5.656854/
83C     ..
84C     .. Executable Statements ..
85C
86      IF (a.EQ.aa) GO TO 10
87      IF (a.LT.1.0) GO TO 130
88C
89C     STEP  1:  RECALCULATIONS OF S2,S,D IF A HAS CHANGED
90C
91      aa = a
92      s2 = a - 0.5
93      s = sqrt(s2)
94      d = sqrt32 - 12.0*s
95C
96C     STEP  2:  T=STANDARD NORMAL DEVIATE,
97C               X=(S,1/2)-NORMAL DEVIATE.
98C               IMMEDIATE ACCEPTANCE (I)
99C
100   10 t = snorm()
101      x = s + 0.5*t
102      sgamma = x*x
103      IF (t.GE.0.0) RETURN
104C
105C     STEP  3:  U= 0,1 -UNIFORM SAMPLE. SQUEEZE ACCEPTANCE (S)
106C
107      u = ranf()
108      IF (d*u.LE.t*t*t) RETURN
109C
110C     STEP  4:  RECALCULATIONS OF Q0,B,SI,C IF NECESSARY
111C
112      IF (a.EQ.aaa) GO TO 40
113      aaa = a
114      r = 1.0/a
115      q0 = ((((((q7*r+q6)*r+q5)*r+q4)*r+q3)*r+q2)*r+q1)*r
116C
117C               APPROXIMATION DEPENDING ON SIZE OF PARAMETER A
118C               THE CONSTANTS IN THE EXPRESSIONS FOR B, SI AND
119C               C WERE ESTABLISHED BY NUMERICAL EXPERIMENTS
120C
121      IF (a.LE.3.686) GO TO 30
122      IF (a.LE.13.022) GO TO 20
123C
124C               CASE 3:  A .GT. 13.022
125C
126      b = 1.77
127      si = .75
128      c = .1515/s
129      GO TO 40
130C
131C               CASE 2:  3.686 .LT. A .LE. 13.022
132C
133   20 b = 1.654 + .0076*s2
134      si = 1.68/s + .275
135      c = .062/s + .024
136      GO TO 40
137C
138C               CASE 1:  A .LE. 3.686
139C
140   30 b = .463 + s + .178*s2
141      si = 1.235
142      c = .195/s - .079 + .16*s
143C
144C     STEP  5:  NO QUOTIENT TEST IF X NOT POSITIVE
145C
146   40 IF (x.LE.0.0) GO TO 70
147C
148C     STEP  6:  CALCULATION OF V AND QUOTIENT Q
149C
150      v = t/ (s+s)
151      IF (abs(v).LE.0.25) GO TO 50
152      q = q0 - s*t + 0.25*t*t + (s2+s2)*alog(1.0+v)
153      GO TO 60
154
155   50 q = q0 + 0.5*t*t* ((((((a7*v+a6)*v+a5)*v+a4)*v+a3)*v+a2)*v+a1)*v
156C
157C     STEP  7:  QUOTIENT ACCEPTANCE (Q)
158C
159   60 IF (alog(1.0-u).LE.q) RETURN
160C
161C     STEP  8:  E=STANDARD EXPONENTIAL DEVIATE
162C               U= 0,1 -UNIFORM DEVIATE
163C               T=(B,SI)-DOUBLE EXPONENTIAL (LAPLACE) SAMPLE
164C
165   70 e = sexpo()
166      u = ranf()
167      u = u + u - 1.0
168      t = b + sign(si*e,u)
169C
170C     STEP  9:  REJECTION IF T .LT. TAU(1) = -.71874483771719
171C
172   80 IF (t.LT. (-.7187449)) GO TO 70
173C
174C     STEP 10:  CALCULATION OF V AND QUOTIENT Q
175C
176      v = t/ (s+s)
177      IF (abs(v).LE.0.25) GO TO 90
178      q = q0 - s*t + 0.25*t*t + (s2+s2)*alog(1.0+v)
179      GO TO 100
180
181   90 q = q0 + 0.5*t*t* ((((((a7*v+a6)*v+a5)*v+a4)*v+a3)*v+a2)*v+a1)*v
182C
183C     STEP 11:  HAT ACCEPTANCE (H) (IF Q NOT POSITIVE GO TO STEP 8)
184C
185  100 IF (q.LE.0.0) GO TO 70
186      IF (q.LE.0.5) GO TO 110
187C
188C     JJV modified the code through line 125 to handle large Q case
189C
190      IF (q.LT.15.0) GO TO 105
191C
192C     JJV Here Q is large enough that Q = log(exp(Q) - 1.0) (for real Q)
193C     JJV so reformulate test at 120 in terms of one EXP, if not too big
194C     JJV 87.49823 is close to the largest real which can be
195C     JJV exponentiated (87.49823 = log(1.0E38))
196C
197      IF ((q+e-0.5*t*t).GT.87.49823) GO TO 125
198      IF (c*abs(u).GT.exp(q+e-0.5*t*t)) GO TO 70
199      GO TO 125
200
201 105  w = exp(q) - 1.0
202      GO TO 120
203
204  110 w = ((((e5*q+e4)*q+e3)*q+e2)*q+e1)*q
205C
206C               IF T IS REJECTED, SAMPLE AGAIN AT STEP 8
207C
208  120 IF (c*abs(u).GT.w*exp(e-0.5*t*t)) GO TO 70
209 125  x = s + 0.5*t
210      sgamma = x*x
211      RETURN
212C
213C     ALTERNATE METHOD FOR PARAMETERS A BELOW 1  (.3678794=EXP(-1.))
214C
215C     JJV changed B to B0 (which was added to declarations for this)
216C     JJV in 130 to END to fix rare and subtle bug.
217C     JJV Line: '130 aa = 0.0' was removed (unnecessary, wasteful).
218C     JJV Reasons: the state of AA only serves to tell the A .GE. 1.0
219C     JJV case if certain A-dependant constants need to be recalculated.
220C     JJV The A .LT. 1.0 case (here) no longer changes any of these, and
221C     JJV the recalculation of B (which used to change with an
222C     JJV A .LT. 1.0 call) is governed by the state of AAA anyway.
223C
224 130  b0 = 1.0 + .3678794*a
225  140 p = b0*ranf()
226      IF (p.GE.1.0) GO TO 150
227      sgamma = exp(alog(p)/a)
228      IF (sexpo().LT.sgamma) GO TO 140
229      RETURN
230
231  150 sgamma = -alog((b0-p)/a)
232      IF (sexpo().LT. (1.0-a)*alog(sgamma)) GO TO 140
233      RETURN
234
235      END
236