1 #include <math.h>
2 #include "grand.h"
3 #include "core_math.h"
4
5
fsign1(double x,double y)6 double fsign1 (double x, double y)
7 {
8 if (y >= 0.0)
9 {
10 return Abs(x);
11 }
12 else
13 {
14 return -Abs(x);
15 }
16 }
17
C2F(sgamma)18 double C2F(sgamma)(double *a)
19 /*
20 **********************************************************************
21
22
23 (STANDARD-) G A M M A DISTRIBUTION
24
25
26 **********************************************************************
27 **********************************************************************
28
29 PARAMETER A >= 1.0 !
30
31 **********************************************************************
32 This source code was taken in the project "freemat"(BSD license)
33 This source code was modified by Gaüzère Sabine according to the
34 modifications done by JJV
35
36 FOR DETAILS SEE:
37
38 AHRENS, J.H. AND DIETER, U.
39 GENERATING GAMMA VARIATES BY A
40 MODIFIED REJECTION TECHNIQUE.
41 COMM. ACM, 25,1 (JAN. 1982), 47 - 54.
42
43 STEP NUMBERS CORRESPOND TO ALGORITHM 'GD' IN THE ABOVE PAPER
44 (STRAIGHTFORWARD IMPLEMENTATION)
45
46 Modified by Barry W. Brown, Feb 3, 1988 to use RANF instead of
47 SUNIF. The argument IR thus goes away.
48
49 **********************************************************************
50
51 PARAMETER 0.0 < A < 1.0 !
52
53 **********************************************************************
54
55 FOR DETAILS SEE:
56
57 AHRENS, J.H. AND DIETER, U.
58 COMPUTER METHODS FOR SAMPLING FROM GAMMA,
59 BETA, POISSON AND BINOMIAL DISTRIBUTIONS.
60 COMPUTING, 12 (1974), 223 - 246.
61
62 (ADAPTED IMPLEMENTATION OF ALGORITHM 'GS' IN THE ABOVE PAPER)
63
64 **********************************************************************
65 INPUT: A =PARAMETER (MEAN) OF THE STANDARD GAMMA DISTRIBUTION
66 OUTPUT: SGAMMA = SAMPLE FROM THE GAMMA-(A)-DISTRIBUTION
67 COEFFICIENTS Q(K) - FOR Q0 = SUM(Q(K)*A**(-K))
68 COEFFICIENTS A(K) - FOR Q = Q0+(T*T/2)*SUM(A(K)*V**K)
69 COEFFICIENTS E(K) - FOR EXP(Q)-1 = SUM(E(K)*Q**K)
70 PREVIOUS A PRE-SET TO ZERO - AA IS A', AAA IS A"
71 SQRT32 IS THE SQUAREROOT OF 32 = 5.656854249492380
72 */
73 {
74 //extern float sign( float num, float sign );
75 static double q1 = 4.166669026017189026E-2;
76 static double q2 = 2.083148062229156494E-2;
77 static double q3 = 8.01191013306379318E-3;
78 static double q4 = 1.44121004268527031E-3;
79 static double q5 = -7.388000085484236E-5;
80 static double q6 = 2.4510998628102243E-4;
81 static double q7 = 2.4239999765995890E-4;
82 static double a1 = 0.33333331346511840820;
83 static double a2 = -0.25000301003456115723;
84 static double a3 = 0.20000620186328887939;
85 static double a4 = -0.16629210114479064941;
86 static double a5 = 0.14236569404602050781;
87 static double a6 = -0.13671770691871643066;
88 static double a7 = 0.12337949872016906738;
89 static double e1 = 1.0;
90 static double e2 = 0.49998968839645385742;
91 static double e3 = 0.16682900488376617432;
92 static double e4 = 4.077529907226562500E-2;
93 static double e5 = 1.029300037771463394E-2;
94 static double aa = 0.0;
95 static double aaa = 0.0;
96 static double sqrt32 = 5.65685415267944335938;
97 static double sgamma, s2, s, d, t, x, u, r, q0, b, si, c, v, q, e, w, p, b0;
98 if (*a == aa)
99 {
100 goto S10;
101 }
102 if (*a < 1.0)
103 {
104 goto S130;
105 }
106 /*
107 STEP 1: RECALCULATIONS OF S2,S,D IF A HAS CHANGED
108 */
109 aa = *a;
110 s2 = *a - 0.5;
111 s = sqrt(s2);
112 d = sqrt32 - 12.0 * s;
113 S10:
114 /*
115 STEP 2: T=STANDARD NORMAL DEVIATE,
116 X=(S,1/2)-NORMAL DEVIATE.
117 IMMEDIATE ACCEPTANCE (I)
118 */
119 t = C2F(snorm)();
120 x = s + 0.5 * t;
121 sgamma = x * x;
122 if (t >= 0.0)
123 {
124 return sgamma;
125 }
126 /*
127 STEP 3: U= 0,1 -UNIFORM SAMPLE. SQUEEZE ACCEPTANCE (S)
128 */
129 u = C2F(ranf)();
130 if (d * u <= t * t * t)
131 {
132 return sgamma;
133 }
134 /*
135 STEP 4: RECALCULATIONS OF Q0,B,SI,C IF NECESSARY
136 */
137 if (*a == aaa)
138 {
139 goto S40;
140 }
141 aaa = *a;
142 r = 1.0 / *a;
143 q0 = ((((((q7 * r + q6) * r + q5) * r + q4) * r + q3) * r + q2) * r + q1) * r;
144 /*
145 APPROXIMATION DEPENDING ON SIZE OF PARAMETER A
146 THE CONSTANTS IN THE EXPRESSIONS FOR B, SI AND
147 C WERE ESTABLISHED BY NUMERICAL EXPERIMENTS
148 */
149 if (*a <= 3.686)
150 {
151 goto S30; //3.68600010871887207031
152 }
153 if (*a <= 13.022)
154 {
155 goto S20; //13.02200031280517578125
156 }
157 /*
158 CASE 3: A .GT. 13.022
159 */
160 b = 1.76999998092651367188;
161 si = 0.75;
162 c = 0.15150000154972076416 / s;
163 goto S40;
164 S20:
165 /*
166 CASE 2: 3.686 .LT. A .LE. 13.022
167 */
168 b = 1.65400004386901855469 + 7.60000012814998627E-3 * s2;
169 si = 1.67999994754791259766 / s + 0.27500000596046447754;
170 c = 6.199999898672103882E-2 / s + 2.400000020861625671E-2;
171 goto S40;
172 S30:
173 /*
174 CASE 1: A .LE. 3.686
175 */
176 b = 0.46299999952316284180 + s + 0.17800000309944152832 * s2;
177 si = 1.23500001430511474609;
178 c = 0.19499999284744262695 / s - 7.900000363588333130E-2 + 1.5999999642372131348E-1 * s;
179 S40:
180 /*
181 STEP 5: NO QUOTIENT TEST IF X NOT POSITIVE
182 */
183 if (x <= 0.0)
184 {
185 goto S70;
186 }
187 /*
188 STEP 6: CALCULATION OF V AND QUOTIENT Q
189 */
190 v = t / (s + s);
191 if (Abs(v) <= 0.25)
192 {
193 goto S50;
194 }
195 q = q0 - s * t + 0.25 * t * t + (s2 + s2) * log(1.0 + v);
196 goto S60;
197 S50:
198 q = q0 + 0.5 * t * t * ((((((a7 * v + a6) * v + a5) * v + a4) * v + a3) * v + a2) * v + a1) * v;
199 S60:
200 /*
201 STEP 7: QUOTIENT ACCEPTANCE (Q)
202 */
203 if (log(1.0 - u) <= q)
204 {
205 return sgamma;
206 }
207 S70:
208 /*
209 STEP 8: E=STANDARD EXPONENTIAL DEVIATE
210 U= 0,1 -UNIFORM DEVIATE
211 T=(B,SI)-DOUBLE EXPONENTIAL (LAPLACE) SAMPLE
212 */
213 e = C2F(sexpo)();
214 u = C2F(ranf)();
215 u += (u - 1.0);
216 t = b + fsign1(si * e, u);
217 /*
218 STEP 9: REJECTION IF T .LT. TAU(1) = -.71874483771719
219 */
220 if (t < -0.7187449)
221 {
222 goto S70; //-0.71874487400054931641
223 }
224 /*
225 STEP 10: CALCULATION OF V AND QUOTIENT Q
226 */
227 v = t / (s + s);
228 if (Abs(v) <= 0.25)
229 {
230 goto S90;
231 }
232 q = q0 - s * t + 0.25 * t * t + (s2 + s2) * log(1.0 + v);
233 goto S100;
234 S90:
235 q = q0 + 0.5 * t * t * ((((((a7 * v + a6) * v + a5) * v + a4) * v + a3) * v + a2) * v + a1) * v;
236 S100:
237 /*
238 STEP 11: HAT ACCEPTANCE (H) (IF Q NOT POSITIVE GO TO STEP 8)
239 */
240 if (q <= 0.0)
241 {
242 goto S70;
243 }
244 if (q <= 0.5)
245 {
246 goto S110;
247 }
248 // JJV modified the code through line 125 to handle large Q case
249 if (q < 15.0)
250 {
251 goto S105;
252 }
253 // JJV Here Q is large enough that Q = log(exp(Q) - 1.0) (for DOUBLE PRECISION Q)
254 // JJV so reformulate test at 120 in terms of one EXP, if not too big
255 // JJV 87.49823 is close to the largest DOUBLE PRECISION which can be
256 // JJV exponentiated (87.49823 = log(1.0E38))
257 if ((q + e - 0.5 * t * t) > 87.49823)
258 {
259 goto S125; //87.498222998046875
260 }
261 if (c * Abs(u) > exp(q + e - 0.5 * t * t))
262 {
263 goto S70;
264 }
265 goto S125;
266 S105:
267 w = exp(q) - 1.0;
268 goto S120;
269 S110:
270 w = ((((e5 * q + e4) * q + e3) * q + e2) * q + e1) * q;
271 S120:
272 /*
273 IF T IS REJECTED, SAMPLE AGAIN AT STEP 8
274 */
275 if (c * Abs(u) > w * exp(e - 0.5 * t * t))
276 {
277 goto S70;
278 }
279 S125:
280 x = s + 0.5 * t;
281 sgamma = x * x;
282 return sgamma;
283 /*
284 ALTERNATE METHOD FOR PARAMETERS A BELOW 1 (.3678794=EXP(-1.))
285 */
286 // JJV changed B to B0 (which was added to declarations for this)
287 // JJV in 130 to END to fix rare and subtle bug.
288 // JJV Line: '130 aa = 0.0' was removed (unnecessary, wasteful).
289 // JJV Reasons: the state of AA only serves to tell the A .GE. 1.0
290 // JJV case if certain A-dependant constants need to be recalculated.
291 // JJV The A .LT. 1.0 case (here) no longer changes any of these, and
292 // JJV the recalculation of B (which used to change with an
293 // JJV A .LT. 1.0 call) is governed by the state of AAA anyway.
294 S130:
295 b0 = 1.0 + 0.3678794 * (*a); //0.36787939071655273438
296 S140:
297 p = b0 * C2F(ranf)();
298 if (p >= 1.0)
299 {
300 goto S150;
301 }
302 sgamma = exp(log(p) / *a);
303 if (C2F(sexpo)() < sgamma)
304 {
305 goto S140;
306 }
307 return sgamma;
308 S150:
309 sgamma = -log((b0 - p) / *a);
310 if (C2F(sexpo)() < (1.0 - *a)*log(sgamma))
311 {
312 goto S140;
313 }
314 return sgamma;
315 }
316