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