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