1 /*
2 * Mathlib : A C Library of Special Functions
3 * Copyright (C) 1998 Ross Ihaka
4 * Copyright (C) 2000-9 The R Core Team
5 *
6 * This program is free software; you can redistribute it and/or modify
7 * it under the terms of the GNU General Public License as published by
8 * the Free Software Foundation; either version 2 of the License, or
9 * (at your option) any later version.
10 *
11 * This program is distributed in the hope that it will be useful,
12 * but WITHOUT ANY WARRANTY; without even the implied warranty of
13 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14 * GNU General Public License for more details.
15 *
16 * You should have received a copy of the GNU General Public License
17 * along with this program; if not, a copy is available at
18 * http://www.r-project.org/Licenses/
19 *
20 * SYNOPSIS
21 *
22 * #include <Rmath.h>
23 * double norm_rand(void);
24 *
25 * DESCRIPTION
26 *
27 * Random variates from the STANDARD normal distribution N(0,1).
28 *
29 * Is called from rnorm(..), but also rt(), rf(), rgamma(), ...
30 */
31
32 #include <R_ext/Random.h>
33 #include "nmath.h"
34
35 #define repeat for(;;)
36
37 #ifdef MATHLIB_STANDALONE
38 static
39 #else
40 attribute_hidden
41 #endif
42 double BM_norm_keep = 0.0;
43
44 N01type N01_kind = INVERSION;
45
46 #ifndef MATHLIB_STANDALONE
47 typedef void * (*DL_FUNC)();
48 extern DL_FUNC User_norm_fun; /* declared and set in ../main/RNG.c */
49 #endif
50
51 /*
52 * REFERENCE
53 *
54 * Ahrens, J.H. and Dieter, U.
55 * Extensions of Forsythe's method for random sampling from
56 * the normal distribution.
57 * Math. Comput. 27, 927-937.
58 *
59 * The definitions of the constants a[k], d[k], t[k] and
60 * h[k] are according to the abovementioned article
61 */
norm_rand(void)62 double norm_rand(void)
63 {
64
65 const static double a[32] =
66 {
67 0.0000000, 0.03917609, 0.07841241, 0.1177699,
68 0.1573107, 0.19709910, 0.23720210, 0.2776904,
69 0.3186394, 0.36012990, 0.40225010, 0.4450965,
70 0.4887764, 0.53340970, 0.57913220, 0.6260990,
71 0.6744898, 0.72451440, 0.77642180, 0.8305109,
72 0.8871466, 0.94678180, 1.00999000, 1.0775160,
73 1.1503490, 1.22985900, 1.31801100, 1.4177970,
74 1.5341210, 1.67594000, 1.86273200, 2.1538750
75 };
76
77 const static double d[31] =
78 {
79 0.0000000, 0.0000000, 0.0000000, 0.0000000,
80 0.0000000, 0.2636843, 0.2425085, 0.2255674,
81 0.2116342, 0.1999243, 0.1899108, 0.1812252,
82 0.1736014, 0.1668419, 0.1607967, 0.1553497,
83 0.1504094, 0.1459026, 0.1417700, 0.1379632,
84 0.1344418, 0.1311722, 0.1281260, 0.1252791,
85 0.1226109, 0.1201036, 0.1177417, 0.1155119,
86 0.1134023, 0.1114027, 0.1095039
87 };
88
89 const static double t[31] =
90 {
91 7.673828e-4, 0.002306870, 0.003860618, 0.005438454,
92 0.007050699, 0.008708396, 0.010423570, 0.012209530,
93 0.014081250, 0.016055790, 0.018152900, 0.020395730,
94 0.022811770, 0.025434070, 0.028302960, 0.031468220,
95 0.034992330, 0.038954830, 0.043458780, 0.048640350,
96 0.054683340, 0.061842220, 0.070479830, 0.081131950,
97 0.094624440, 0.112300100, 0.136498000, 0.171688600,
98 0.227624100, 0.330498000, 0.584703100
99 };
100
101 const static double h[31] =
102 {
103 0.03920617, 0.03932705, 0.03950999, 0.03975703,
104 0.04007093, 0.04045533, 0.04091481, 0.04145507,
105 0.04208311, 0.04280748, 0.04363863, 0.04458932,
106 0.04567523, 0.04691571, 0.04833487, 0.04996298,
107 0.05183859, 0.05401138, 0.05654656, 0.05953130,
108 0.06308489, 0.06737503, 0.07264544, 0.07926471,
109 0.08781922, 0.09930398, 0.11555990, 0.14043440,
110 0.18361420, 0.27900160, 0.70104740
111 };
112
113 /*----------- Constants and definitions for Kinderman - Ramage --- */
114 /*
115 * REFERENCE
116 *
117 * Kinderman A. J. and Ramage J. G. (1976).
118 * Computer generation of normal random variables.
119 * JASA 71, 893-896.
120 */
121
122 #define C1 0.398942280401433
123 #define C2 0.180025191068563
124 #define g(x) (C1*exp(-x*x/2.0)-C2*(A-x))
125
126 const static double A = 2.216035867166471;
127
128 double s, u1, w, y, u2, u3, aa, tt, theta, R;
129 int i;
130
131 switch(N01_kind) {
132
133 case AHRENS_DIETER: /* see Reference above */
134
135 u1 = unif_rand();
136 s = 0.0;
137 if (u1 > 0.5)
138 s = 1.0;
139 u1 = u1 + u1 - s;
140 u1 *= 32.0;
141 i = (int) u1;
142 if (i == 32)
143 i = 31;
144 if (i != 0) {
145 u2 = u1 - i;
146 aa = a[i - 1];
147 while (u2 <= t[i - 1]) {
148 u1 = unif_rand();
149 w = u1 * (a[i] - aa);
150 tt = (w * 0.5 + aa) * w;
151 repeat {
152 if (u2 > tt)
153 goto deliver;
154 u1 = unif_rand();
155 if (u2 < u1)
156 break;
157 tt = u1;
158 u2 = unif_rand();
159 }
160 u2 = unif_rand();
161 }
162 w = (u2 - t[i - 1]) * h[i - 1];
163 }
164 else {
165 i = 6;
166 aa = a[31];
167 repeat {
168 u1 = u1 + u1;
169 if (u1 >= 1.0)
170 break;
171 aa = aa + d[i - 1];
172 i = i + 1;
173 }
174 u1 = u1 - 1.0;
175 repeat {
176 w = u1 * d[i - 1];
177 tt = (w * 0.5 + aa) * w;
178 repeat {
179 u2 = unif_rand();
180 if (u2 > tt)
181 goto jump;
182 u1 = unif_rand();
183 if (u2 < u1)
184 break;
185 tt = u1;
186 }
187 u1 = unif_rand();
188 }
189 jump:;
190 }
191
192 deliver:
193 y = aa + w;
194 return (s == 1.0) ? -y : y;
195
196 /*-----------------------------------------------------------*/
197
198 case BUGGY_KINDERMAN_RAMAGE: /* see Reference above */
199 /* note: this has problems, but is retained for
200 * reproducibility of older codes, with the same
201 * numeric code */
202 u1 = unif_rand();
203 if(u1 < 0.884070402298758) {
204 u2 = unif_rand();
205 return A*(1.13113163544180*u1+u2-1);
206 }
207
208 if(u1 >= 0.973310954173898) { /* tail: */
209 repeat {
210 u2 = unif_rand();
211 u3 = unif_rand();
212 tt = (A*A-2*log(u3));
213 if( u2*u2<(A*A)/tt )
214 return (u1 < 0.986655477086949) ? sqrt(tt) : -sqrt(tt);
215 }
216 }
217
218 if(u1 >= 0.958720824790463) { /* region3: */
219 repeat {
220 u2 = unif_rand();
221 u3 = unif_rand();
222 tt = A - 0.630834801921960* fmin2(u2,u3);
223 if(fmax2(u2,u3) <= 0.755591531667601)
224 return (u2<u3) ? tt : -tt;
225 if(0.034240503750111*fabs(u2-u3) <= g(tt))
226 return (u2<u3) ? tt : -tt;
227 }
228 }
229
230 if(u1 >= 0.911312780288703) { /* region2: */
231 repeat {
232 u2 = unif_rand();
233 u3 = unif_rand();
234 tt = 0.479727404222441+1.105473661022070*fmin2(u2,u3);
235 if( fmax2(u2,u3)<=0.872834976671790 )
236 return (u2<u3) ? tt : -tt;
237 if( 0.049264496373128*fabs(u2-u3)<=g(tt) )
238 return (u2<u3) ? tt : -tt;
239 }
240 }
241
242 /* ELSE region1: */
243 repeat {
244 u2 = unif_rand();
245 u3 = unif_rand();
246 tt = 0.479727404222441-0.595507138015940*fmin2(u2,u3);
247 if(fmax2(u2,u3) <= 0.805577924423817)
248 return (u2<u3) ? tt : -tt;
249 }
250 case BOX_MULLER:
251 if(BM_norm_keep != 0.0) { /* An exact test is intentional */
252 s = BM_norm_keep;
253 BM_norm_keep = 0.0;
254 return s;
255 } else {
256 theta = 2 * M_PI * unif_rand();
257 R = sqrt(-2 * log(unif_rand())) + 10*DBL_MIN; /* ensure non-zero */
258 BM_norm_keep = R * sin(theta);
259 return R * cos(theta);
260 }
261 #ifndef MATHLIB_STANDALONE
262 case USER_NORM:
263 return *((double *) User_norm_fun());
264 #endif
265 case INVERSION:
266 #define BIG 134217728 /* 2^27 */
267 /* unif_rand() alone is not of high enough precision */
268 u1 = unif_rand();
269 u1 = (int)(BIG*u1) + unif_rand();
270 return qnorm5(u1/BIG, 0.0, 1.0, 1, 0);
271 case KINDERMAN_RAMAGE: /* see Reference above */
272 /* corrected version from Josef Leydold
273 * */
274 u1 = unif_rand();
275 if(u1 < 0.884070402298758) {
276 u2 = unif_rand();
277 return A*(1.131131635444180*u1+u2-1);
278 }
279
280 if(u1 >= 0.973310954173898) { /* tail: */
281 repeat {
282 u2 = unif_rand();
283 u3 = unif_rand();
284 tt = (A*A-2*log(u3));
285 if( u2*u2<(A*A)/tt )
286 return (u1 < 0.986655477086949) ? sqrt(tt) : -sqrt(tt);
287 }
288 }
289
290 if(u1 >= 0.958720824790463) { /* region3: */
291 repeat {
292 u2 = unif_rand();
293 u3 = unif_rand();
294 tt = A - 0.630834801921960* fmin2(u2,u3);
295 if(fmax2(u2,u3) <= 0.755591531667601)
296 return (u2<u3) ? tt : -tt;
297 if(0.034240503750111*fabs(u2-u3) <= g(tt))
298 return (u2<u3) ? tt : -tt;
299 }
300 }
301
302 if(u1 >= 0.911312780288703) { /* region2: */
303 repeat {
304 u2 = unif_rand();
305 u3 = unif_rand();
306 tt = 0.479727404222441+1.105473661022070*fmin2(u2,u3);
307 if( fmax2(u2,u3)<=0.872834976671790 )
308 return (u2<u3) ? tt : -tt;
309 if( 0.049264496373128*fabs(u2-u3)<=g(tt) )
310 return (u2<u3) ? tt : -tt;
311 }
312 }
313
314 /* ELSE region1: */
315 repeat {
316 u2 = unif_rand();
317 u3 = unif_rand();
318 tt = 0.479727404222441-0.595507138015940*fmin2(u2,u3);
319 if (tt < 0.) continue;
320 if(fmax2(u2,u3) <= 0.805577924423817)
321 return (u2<u3) ? tt : -tt;
322 if(0.053377549506886*fabs(u2-u3) <= g(tt))
323 return (u2<u3) ? tt : -tt;
324 }
325 default:
326 MATHLIB_ERROR(_("norm_rand(): invalid N01_kind: %d\n"), N01_kind)
327 return 0.0;/*- -Wall */
328 }/*switch*/
329 }
330