1 /* This is third-party software that is distributed with UTILIB.
2 * For licensing information concerning this file, see the UTILIB home page:
3 * http://software.sandia.gov/trac/utilib
4 */
5
6 /* sgamma.c
7 *
8 */
9
10 #ifdef HAVE_FSIGN
11
12 #include <utilib_config.h>
13 #include <math.h>
14 #include <utilib/Random.h>
15
16
sgamma(double a)17 double sgamma(double a)
18 /*
19 **********************************************************************
20
21
22 (STANDARD-) G A M M A DISTRIBUTION
23
24
25 **********************************************************************
26 **********************************************************************
27
28 PARAMETER A >= 1.0 !
29
30 **********************************************************************
31
32 FOR DETAILS SEE:
33
34 AHRENS, J.H. AND DIETER, U.
35 GENERATING GAMMA VARIATES BY A
36 MODIFIED REJECTION TECHNIQUE.
37 COMM. ACM, 25,1 (JAN. 1982), 47 - 54.
38
39 STEP NUMBERS CORRESPOND TO ALGORITHM 'GD' IN THE ABOVE PAPER
40 (STRAIGHTFORWARD IMPLEMENTATION)
41
42 Modified by Barry W. Brown, Feb 3, 1988 to use RANF instead of
43 SUNIF. The argument IR thus goes away.
44
45 **********************************************************************
46
47 PARAMETER 0.0 < A < 1.0 !
48
49 **********************************************************************
50
51 FOR DETAILS SEE:
52
53 AHRENS, J.H. AND DIETER, U.
54 COMPUTER METHODS FOR SAMPLING FROM GAMMA,
55 BETA, POISSON AND BINOMIAL DISTRIBUTIONS.
56 COMPUTING, 12 (1974), 223 - 246.
57
58 (ADAPTED IMPLEMENTATION OF ALGORITHM 'GS' IN THE ABOVE PAPER)
59
60 **********************************************************************
61 INPUT: A =PARAMETER (MEAN) OF THE STANDARD GAMMA DISTRIBUTION
62 OUTPUT: SGAMMA = SAMPLE FROM THE GAMMA-(A)-DISTRIBUTION
63 COEFFICIENTS Q(K) - FOR Q0 = SUM(Q(K)*A**(-K))
64 COEFFICIENTS A(K) - FOR Q = Q0+(T*T/2)*SUM(A(K)*V**K)
65 COEFFICIENTS E(K) - FOR EXP(Q)-1 = SUM(E(K)*Q**K)
66 PREVIOUS A PRE-SET TO ZERO - AA IS A', AAA IS A"
67 SQRT32 IS THE SQUAREROOT OF 32 = 5.656854249492380
68 */
69 {
70 extern double fsign( double num, double sign );
71 static double q1 = 4.166669E-2;
72 static double q2 = 2.083148E-2;
73 static double q3 = 8.01191E-3;
74 static double q4 = 1.44121E-3;
75 static double q5 = -7.388E-5;
76 static double q6 = 2.4511E-4;
77 static double q7 = 2.424E-4;
78 static double a1 = 0.3333333;
79 static double a2 = -0.250003;
80 static double a3 = 0.2000062;
81 static double a4 = -0.1662921;
82 static double a5 = 0.1423657;
83 static double a6 = -0.1367177;
84 static double a7 = 0.1233795;
85 static double e1 = 1.0;
86 static double e2 = 0.4999897;
87 static double e3 = 0.166829;
88 static double e4 = 4.07753E-2;
89 static double e5 = 1.0293E-2;
90 static double aa = 0.0;
91 static double aaa = 0.0;
92 static double sqrt32 = 5.656854;
93 static double sgamma,s2,s,d,t,x,u,r,q0,b,si,c,v,q,e,w,p;
94 if(a == aa) goto S10;
95 if(a < 1.0) goto S120;
96 /*
97 STEP 1: RECALCULATIONS OF S2,S,D IF A HAS CHANGED
98 */
99 aa = a;
100 s2 = a-0.5;
101 s = sqrt(s2);
102 d = sqrt32-12.0*s;
103 S10:
104 /*
105 STEP 2: T=STANDARD NORMAL DEVIATE,
106 X=(S,1/2)-NORMAL DEVIATE.
107 IMMEDIATE ACCEPTANCE (I)
108 */
109 t = snorm();
110 x = s+0.5*t;
111 sgamma = x*x;
112 if(t >= 0.0) return sgamma;
113 /*
114 STEP 3: U= 0,1 -UNIFORM SAMPLE. SQUEEZE ACCEPTANCE (S)
115 */
116 u = ranf();
117 if(d*u <= t*t*t) return sgamma;
118 /*
119 STEP 4: RECALCULATIONS OF Q0,B,SI,C IF NECESSARY
120 */
121 if(a == aaa) goto S40;
122 aaa = a;
123 r = 1.0/ a;
124 q0 = ((((((q7*r+q6)*r+q5)*r+q4)*r+q3)*r+q2)*r+q1)*r;
125 /*
126 APPROXIMATION DEPENDING ON SIZE OF PARAMETER A
127 THE CONSTANTS IN THE EXPRESSIONS FOR B, SI AND
128 C WERE ESTABLISHED BY NUMERICAL EXPERIMENTS
129 */
130 if(a <= 3.686) goto S30;
131 if(a <= 13.022) goto S20;
132 /*
133 CASE 3: A .GT. 13.022
134 */
135 b = 1.77;
136 si = 0.75;
137 c = 0.1515/s;
138 goto S40;
139 S20:
140 /*
141 CASE 2: 3.686 .LT. A .LE. 13.022
142 */
143 b = 1.654+7.6E-3*s2;
144 si = 1.68/s+0.275;
145 c = 6.2E-2/s+2.4E-2;
146 goto S40;
147 S30:
148 /*
149 CASE 1: A .LE. 3.686
150 */
151 b = 0.463+s-0.178*s2;
152 si = 1.235;
153 c = 0.195/s-7.9E-2+1.6E-2*s;
154 S40:
155 /*
156 STEP 5: NO QUOTIENT TEST IF X NOT POSITIVE
157 */
158 if(x <= 0.0) goto S70;
159 /*
160 STEP 6: CALCULATION OF V AND QUOTIENT Q
161 */
162 v = t/(s+s);
163 if(fabs(v) <= 0.25) goto S50;
164 q = q0-s*t+0.25*t*t+(s2+s2)*log(1.0+v);
165 goto S60;
166 S50:
167 q = q0+0.5*t*t*((((((a7*v+a6)*v+a5)*v+a4)*v+a3)*v+a2)*v+a1)*v;
168 S60:
169 /*
170 STEP 7: QUOTIENT ACCEPTANCE (Q)
171 */
172 if(log(1.0-u) <= q) return sgamma;
173 S70:
174 /*
175 STEP 8: E=STANDARD EXPONENTIAL DEVIATE
176 U= 0,1 -UNIFORM DEVIATE
177 T=(B,SI)-DOUBLE EXPONENTIAL (LAPLACE) SAMPLE
178 */
179 e = sexpo();
180 u = ranf();
181 u += (u-1.0);
182 t = b+fsign(si*e,u);
183 /*
184 STEP 9: REJECTION IF T .LT. TAU(1) = -.71874483771719
185 */
186 if(t < -0.7187449) goto S70;
187 /*
188 STEP 10: CALCULATION OF V AND QUOTIENT Q
189 */
190 v = t/(s+s);
191 if(fabs(v) <= 0.25) goto S80;
192 q = q0-s*t+0.25*t*t+(s2+s2)*log(1.0+v);
193 goto S90;
194 S80:
195 q = q0+0.5*t*t*((((((a7*v+a6)*v+a5)*v+a4)*v+a3)*v+a2)*v+a1)*v;
196 S90:
197 /*
198 STEP 11: HAT ACCEPTANCE (H) (IF Q NOT POSITIVE GO TO STEP 8)
199 */
200 if(q <= 0.0) goto S70;
201 if(q <= 0.5) goto S100;
202 w = exp(q)-1.0;
203 goto S110;
204 S100:
205 w = ((((e5*q+e4)*q+e3)*q+e2)*q+e1)*q;
206 S110:
207 /*
208 IF T IS REJECTED, SAMPLE AGAIN AT STEP 8
209 */
210 if(c*fabs(u) > w*exp(e-0.5*t*t)) goto S70;
211 x = s+0.5*t;
212 sgamma = x*x;
213 return sgamma;
214 S120:
215 /*
216 ALTERNATE METHOD FOR PARAMETERS A BELOW 1 (.3678794=EXP(-1.))
217 */
218 aa = 0.0;
219 b = 1.0+0.3678794*a;
220 S130:
221 p = b*ranf();
222 if(p >= 1.0) goto S140;
223 sgamma = exp(log(p)/ a);
224 if(sexpo() < sgamma) goto S130;
225 return sgamma;
226 S140:
227 sgamma = -log((b-p)/ a);
228 if(sexpo() < (1.0-a)*log(sgamma)) goto S130;
229 return sgamma;
230 }
231 #endif
232
233