1 #include <config.h>
2 #include <rng/RmathRNG.h>
3
4 #include <cmath>
5 #include <cfloat>
6 #include <stdexcept>
7
8 #define repeat for(;;)
9
10 using std::string;
11 using std::log;
12 using std::exp;
13 using std::sin;
14 using std::cos;
15 using std::sqrt;
16 using std::fabs;
17 using std::logic_error;
18
19 #define PI 3.141592653589793238462643383280
20
21 namespace jags {
22
fmin2(double x,double y)23 static inline double fmin2(double x, double y) {
24 return (x < y) ? x : y;
25 }
26
fmax2(double x,double y)27 static inline double fmax2(double x, double y)
28 {
29 return (x < y) ? y : x;
30 }
31
RmathRNG(string const & name,NormKind N01_kind)32 RmathRNG::RmathRNG(string const &name, NormKind N01_kind)
33 : RNG(name), _N01_kind(N01_kind), _BM_norm_keep(0)
34 {}
35
exponential()36 double RmathRNG::exponential()
37 {
38 /* q[k-1] = sum(log(2)^k / k!) k=1,..,n, */
39 /* The highest n (here 8) is determined by q[n-1] = 1.0 */
40 /* within standard precision */
41 const static double q[] =
42 {
43 0.6931471805599453,
44 0.9333736875190459,
45 0.9888777961838675,
46 0.9984959252914960,
47 0.9998292811061389,
48 0.9999833164100727,
49 0.9999985691438767,
50 0.9999998906925558,
51 0.9999999924734159,
52 0.9999999995283275,
53 0.9999999999728814,
54 0.9999999999985598,
55 0.9999999999999289,
56 0.9999999999999968,
57 0.9999999999999999,
58 1.0000000000000000
59 };
60 double a, u, ustar, umin;
61 int i;
62
63 a = 0.;
64 /* precaution if u = 0 is ever returned */
65 u = uniform();
66 while(u <= 0.0 || u >= 1.0) u = uniform();
67 for (;;) {
68 u += u;
69 if (u > 1.0)
70 break;
71 a += q[0];
72 }
73 u -= 1.;
74
75 if (u <= q[0])
76 return a + u;
77
78 i = 0;
79 ustar = uniform();
80 umin = ustar;
81 do {
82 ustar = uniform();
83 if (ustar < umin)
84 umin = ustar;
85 i++;
86 } while (u > q[i]);
87 return a + umin * q[0];
88 }
89
90 /*
91 * REFERENCE
92 *
93 * Ahrens, J.H. and Dieter, U.
94 * Extensions of Forsythe's method for random sampling from
95 * the normal distribution.
96 * Math. Comput. 27, 927-937.
97 *
98 * The definitions of the constants a[k], d[k], t[k] and
99 * h[k] are according to the abovementioned article
100 */
normal()101 double RmathRNG::normal()
102 {
103
104 const static double a[32] =
105 {
106 0.0000000, 0.03917609, 0.07841241, 0.1177699,
107 0.1573107, 0.19709910, 0.23720210, 0.2776904,
108 0.3186394, 0.36012990, 0.40225010, 0.4450965,
109 0.4887764, 0.53340970, 0.57913220, 0.6260990,
110 0.6744898, 0.72451440, 0.77642180, 0.8305109,
111 0.8871466, 0.94678180, 1.00999000, 1.0775160,
112 1.1503490, 1.22985900, 1.31801100, 1.4177970,
113 1.5341210, 1.67594000, 1.86273200, 2.1538750
114 };
115
116 const static double d[31] =
117 {
118 0.0000000, 0.0000000, 0.0000000, 0.0000000,
119 0.0000000, 0.2636843, 0.2425085, 0.2255674,
120 0.2116342, 0.1999243, 0.1899108, 0.1812252,
121 0.1736014, 0.1668419, 0.1607967, 0.1553497,
122 0.1504094, 0.1459026, 0.1417700, 0.1379632,
123 0.1344418, 0.1311722, 0.1281260, 0.1252791,
124 0.1226109, 0.1201036, 0.1177417, 0.1155119,
125 0.1134023, 0.1114027, 0.1095039
126 };
127
128 const static double t[31] =
129 {
130 7.673828e-4, 0.002306870, 0.003860618, 0.005438454,
131 0.007050699, 0.008708396, 0.010423570, 0.012209530,
132 0.014081250, 0.016055790, 0.018152900, 0.020395730,
133 0.022811770, 0.025434070, 0.028302960, 0.031468220,
134 0.034992330, 0.038954830, 0.043458780, 0.048640350,
135 0.054683340, 0.061842220, 0.070479830, 0.081131950,
136 0.094624440, 0.112300100, 0.136498000, 0.171688600,
137 0.227624100, 0.330498000, 0.584703100
138 };
139
140 const static double h[31] =
141 {
142 0.03920617, 0.03932705, 0.03950999, 0.03975703,
143 0.04007093, 0.04045533, 0.04091481, 0.04145507,
144 0.04208311, 0.04280748, 0.04363863, 0.04458932,
145 0.04567523, 0.04691571, 0.04833487, 0.04996298,
146 0.05183859, 0.05401138, 0.05654656, 0.05953130,
147 0.06308489, 0.06737503, 0.07264544, 0.07926471,
148 0.08781922, 0.09930398, 0.11555990, 0.14043440,
149 0.18361420, 0.27900160, 0.70104740
150 };
151
152 /*----------- Constants and definitions for Kinderman - Ramage --- */
153 /*
154 * REFERENCE
155 *
156 * Kinderman A. J. and Ramage J. G. (1976).
157 * Computer generation of normal random variables.
158 * JASA 71, 893-896.
159 */
160
161 #define C1 0.398942280401433
162 #define C2 0.180025191068563
163 #define g(x) (C1*exp(-x*x/2.0)-C2*(A-x))
164
165 const static double A = 2.216035867166471;
166
167 double s, u1, w, y, u2, u3, aa, tt, theta, R;
168 int i;
169
170 switch(_N01_kind) {
171
172 case AHRENS_DIETER: /* see Reference above */
173
174 u1 = uniform();
175 s = 0.0;
176 if (u1 > 0.5)
177 s = 1.0;
178 u1 = u1 + u1 - s;
179 u1 *= 32.0;
180 i = (int) u1;
181 if (i == 32)
182 i = 31;
183 if (i != 0) {
184 u2 = u1 - i;
185 aa = a[i - 1];
186 while (u2 <= t[i - 1]) {
187 u1 = uniform();
188 w = u1 * (a[i] - aa);
189 tt = (w * 0.5 + aa) * w;
190 repeat {
191 if (u2 > tt)
192 goto deliver;
193 u1 = uniform();
194 if (u2 < u1)
195 break;
196 tt = u1;
197 u2 = uniform();
198 }
199 u2 = uniform();
200 }
201 w = (u2 - t[i - 1]) * h[i - 1];
202 }
203 else {
204 i = 6;
205 aa = a[31];
206 repeat {
207 u1 = u1 + u1;
208 if (u1 >= 1.0)
209 break;
210 aa = aa + d[i - 1];
211 i = i + 1;
212 }
213 u1 = u1 - 1.0;
214 repeat {
215 w = u1 * d[i - 1];
216 tt = (w * 0.5 + aa) * w;
217 repeat {
218 u2 = uniform();
219 if (u2 > tt)
220 goto jump;
221 u1 = uniform();
222 if (u2 < u1)
223 break;
224 tt = u1;
225 }
226 u1 = uniform();
227 }
228 jump:;
229 }
230
231 deliver:
232 y = aa + w;
233 return (s == 1.0) ? -y : y;
234
235 case BOX_MULLER:
236 if(_BM_norm_keep != 0.0) { /* An exact test is intentional */
237 s = _BM_norm_keep;
238 _BM_norm_keep = 0.0;
239 return s;
240 } else {
241 theta = 2 * PI * uniform();
242 R = sqrt(-2 * log(uniform())) + 10*DBL_MIN; /* ensure non-zero */
243 _BM_norm_keep = R * sin(theta);
244 return R * cos(theta);
245 }
246
247 case KINDERMAN_RAMAGE: /* see Reference above */
248 /* corrected version from Josef Leydold
249 * */
250 u1 = uniform();
251 if(u1 < 0.884070402298758) {
252 u2 = uniform();
253 return A*(1.131131635444180*u1+u2-1);
254 }
255
256 if(u1 >= 0.973310954173898) { /* tail: */
257 repeat {
258 u2 = uniform();
259 u3 = uniform();
260 tt = (A*A-2*log(u3));
261 if( u2*u2<(A*A)/tt )
262 return (u1 < 0.986655477086949) ? sqrt(tt) : -sqrt(tt);
263 }
264 }
265
266 if(u1 >= 0.958720824790463) { /* region3: */
267 repeat {
268 u2 = uniform();
269 u3 = uniform();
270 tt = A - 0.630834801921960* fmin2(u2,u3);
271 if(fmax2(u2,u3) <= 0.755591531667601)
272 return (u2<u3) ? tt : -tt;
273 if(0.034240503750111*fabs(u2-u3) <= g(tt))
274 return (u2<u3) ? tt : -tt;
275 }
276 }
277
278 if(u1 >= 0.911312780288703) { /* region2: */
279 repeat {
280 u2 = uniform();
281 u3 = uniform();
282 tt = 0.479727404222441+1.105473661022070*fmin2(u2,u3);
283 if( fmax2(u2,u3)<=0.872834976671790 )
284 return (u2<u3) ? tt : -tt;
285 if( 0.049264496373128*fabs(u2-u3)<=g(tt) )
286 return (u2<u3) ? tt : -tt;
287 }
288 }
289
290 /* ELSE region1: */
291 repeat {
292 u2 = uniform();
293 u3 = uniform();
294 tt = 0.479727404222441-0.595507138015940*fmin2(u2,u3);
295 if (tt < 0.) continue;
296 if(fmax2(u2,u3) <= 0.805577924423817)
297 return (u2<u3) ? tt : -tt;
298 if(0.053377549506886*fabs(u2-u3) <= g(tt))
299 return (u2<u3) ? tt : -tt;
300 }
301
302 }/*switch*/
303
304 // Not reached, but an exit statement is required for -Wall
305 throw logic_error("Bad exit from RmathRNG::normal");
306 return 0;
307 }
308
309 } //namespace jags
310