1 /*
2 * R : A Computer Language for Statistical Data Analysis
3 * Copyright (C) 1995, 1996 Robert Gentleman and Ross Ihaka
4 * Copyright (C) 1997--2007 The R Development Core Team
5 *
6 * Copyright (C) 2006 - 2011 Dirk Eddelbuettel (dieharder adaptation)
7 *
8 *
9 * This program is free software; you can redistribute it and/or modify
10 * it under the terms of the GNU General Public License as published by
11 * the Free Software Foundation; either version 2 of the License, or
12 * (at your option) any later version.
13 *
14 * This program is distributed in the hope that it will be useful,
15 * but WITHOUT ANY WARRANTY; without even the implied warranty of
16 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
17 * GNU General Public License for more details.
18 *
19 * You should have received a copy of the GNU General Public License
20 * along with this program; if not, a copy is available at
21 * http://www.r-project.org/Licenses/
22 */
23
24 /*
25 * A few defines to select 'dieharder' mode (which is what we use here and
26 * 'GNU_R_MODE' (which is how we 'turn off' the features used by R but not
27 * here.
28 */
29 #ifndef DIEHARDER
30 #define DIEHARDER
31 #ifdef GNU_R_MODE
32 #undef GNU_R_MODE
33 #endif
34 #endif
35
36 #ifdef DIEHARDER
37 #include <dieharder/libdieharder.h>
38 /* we need M here for the MT, but dieharder also defines it */
39 #ifdef M
40 #undef M
41 #endif
42
43
44 /* define our function prototypes */
45 /* two generic 'getters' used by all GNU R RNGs */
46 static unsigned long int r_rng_get (void *vstate);
47 static double r_rng_get_double (void *vstate);
48 /* RNG-specific 'setters' for each GNU R RNG */
49 static void r_wichmann_hill_set (void *vstate, unsigned long int s);
50 static void r_marsaglia_mc_set (void *vstate, unsigned long int s);
51 static void r_super_duper_set (void *vstate, unsigned long int s);
52 static void r_mersenne_twister_set (void *vstate, unsigned long int s);
53 static void r_knuth_taocp_set (void *vstate, unsigned long int s);
54 static void r_knuth_taocp2_set (void *vstate, unsigned long int s);
55
56 typedef unsigned int Int32; /* in Random.h */
57
58 typedef void * (*DL_FUNC)();
59
60 /*
61 * The question is, will the patch below break the Debian
62 * build process? Had better test it as best I can in my Debian VM.
63 void error(const char *txt, ...);
64 char *_(char *txt) { return(txt); };
65 */
dieharder_error(const char * format,...)66 void dieharder_error(const char *format, ...)
67 {
68
69 va_list ap;
70 va_start(ap, format);
71 vfprintf(stderr, format, ap);
72 va_end(ap);
73
74 }
75
76 #define error dieharder_error
77 #define _(x) x
78
79
80 /* UINT_MAX from limits.h */
81 #define RANDNUM_MAX UINT_MAX
82
83 /* from R_ext/Random.h */
84 typedef enum {
85 WICHMANN_HILL,
86 MARSAGLIA_MULTICARRY,
87 SUPER_DUPER,
88 MERSENNE_TWISTER,
89 KNUTH_TAOCP,
90 USER_UNIF,
91 KNUTH_TAOCP2
92 } RNGtype;
93
94 /* Different kinds of "N(0,1)" generators :*/
95 typedef enum {
96 BUGGY_KINDERMAN_RAMAGE,
97 AHRENS_DIETER,
98 BOX_MULLER,
99 USER_NORM,
100 INVERSION,
101 KINDERMAN_RAMAGE
102 } N01type;
103 #endif
104
105 #ifdef GNU_R_MODE
106 /* <UTF8> char here is handled as a whole string */
107
108 #ifdef HAVE_CONFIG_H
109 #include <config.h>
110
111 #include "Defn.h"
112 #include <R_ext/Random.h>
113 #endif
114 #endif
115
116 /* Normal generator is not actually set here but in nmath/snorm.c */
117 #define RNG_DEFAULT MERSENNE_TWISTER
118 #define N01_DEFAULT INVERSION
119
120 #ifdef GNU_R_MODE
121 #include <R_ext/Rdynload.h>
122
123 static DL_FUNC User_unif_fun, User_unif_init, User_unif_nseed,
124 User_unif_seedloc;
125
126 DL_FUNC User_norm_fun; /* also in ../nmath/snorm.c */
127 #endif
128
129 static RNGtype RNG_kind = RNG_DEFAULT;
130 extern N01type N01_kind; /* from ../nmath/snorm.c */
131
132 /* typedef unsigned int Int32; in Random.h */
133
134 /* .Random.seed == (RNGkind, i_seed[0],i_seed[1],..,i_seed[n_seed-1])
135 * or == (RNGkind) or missing [--> Randomize]
136 */
137
138 typedef struct {
139 RNGtype kind;
140 N01type Nkind;
141 char *name; /* print name */
142 int n_seed; /* length of seed vector */
143 Int32 *i_seed;
144 } RNGTAB;
145
146 #ifdef DIEHARDER
147 typedef RNGTAB r_rng_tab_t;
148
149 #define TAOCP_MAX 0x3fffffffUL /* cf GSL sources */
150 #define TAOCP2_MAX (1L<<30)-1L /* cf GSL sources */
151 #endif
152
153 static Int32 dummyvec[625]; /* s/dummy/dummyvec/ */
154 static
155 RNGTAB RNG_Table[] =
156 {
157 /* kind Nkind name n_seed i_seed */
158 { 0, 0, "Wichmann-Hill", 3, dummyvec},
159 { 1, 0, "Marsaglia-MultiCarry", 2, dummyvec},
160 { 2, 0, "Super-Duper", 2, dummyvec},
161 { 3, 0, "Mersenne-Twister", 1+624, dummyvec},
162 { 4, 0, "Knuth-TAOCP", 1+100, dummyvec},
163 { 5, 0, "User-supplied", 0, dummyvec},
164 { 6, 0, "Knuth-TAOCP-2002", 1+100, dummyvec},
165 };
166
167
168 #define d2_32 4294967296./* = (double) */
169 #define i2_32m1 2.328306437080797e-10/* = 1/(2^32 - 1) */
170 #define KT 9.31322574615479e-10 /* = 2^-30 */
171
172 #define I1 (RNG_Table[RNG_kind].i_seed[0])
173 #define I2 (RNG_Table[RNG_kind].i_seed[1])
174 #define I3 (RNG_Table[RNG_kind].i_seed[2])
175
176 static void Randomize(RNGtype kind);
177 static double MT_genrand(void);
178 static Int32 KT_next(void);
179 static void RNG_Init_R_KT(Int32);
180 static void RNG_Init_KT2(Int32);
181 #define KT_pos (RNG_Table[KNUTH_TAOCP].i_seed[100])
182
fixup(double x)183 static double fixup(double x)
184 {
185 /* ensure 0 and 1 are never returned */
186 if(x <= 0.0) return 0.5*i2_32m1;
187 if((1.0 - x) <= 0.0) return 1.0 - 0.5*i2_32m1;
188 return x;
189 }
190
191
unif_rand(void)192 double unif_rand(void)
193 {
194 double value;
195
196 switch(RNG_kind) {
197
198 case WICHMANN_HILL:
199 I1 = I1 * 171 % 30269;
200 I2 = I2 * 172 % 30307;
201 I3 = I3 * 170 % 30323;
202 value = I1 / 30269.0 + I2 / 30307.0 + I3 / 30323.0;
203 return fixup(value - (int) value);/* in [0,1) */
204
205 case MARSAGLIA_MULTICARRY:/* 0177777(octal) == 65535(decimal)*/
206 I1= 36969*(I1 & 0177777) + (I1>>16);
207 I2= 18000*(I2 & 0177777) + (I2>>16);
208 return fixup(((I1 << 16)^(I2 & 0177777)) * i2_32m1); /* in [0,1) */
209
210 case SUPER_DUPER:
211 /* This is Reeds et al (1984) implementation;
212 * modified using __unsigned__ seeds instead of signed ones
213 */
214 I1 ^= ((I1 >> 15) & 0377777); /* Tausworthe */
215 I1 ^= I1 << 17;
216 I2 *= 69069; /* Congruential */
217 return fixup((I1^I2) * i2_32m1); /* in [0,1) */
218
219 case MERSENNE_TWISTER:
220 return fixup(MT_genrand());
221
222 case KNUTH_TAOCP:
223 case KNUTH_TAOCP2:
224 return fixup(KT_next() * KT);
225
226 #if GNU_R_MODE
227 case USER_UNIF:
228 return *((double *) User_unif_fun());
229 #endif
230
231 default:
232 error(_("unif_rand: unimplemented RNG kind %d"), RNG_kind);
233 return -1.;
234 }
235 }
236
237 /* we must mask global variable here, as I1-I3 hide RNG_kind
238 and we want the argument */
FixupSeeds(RNGtype RNG_kind,int initial)239 static void FixupSeeds(RNGtype RNG_kind, int initial)
240 {
241 /* Depending on RNG, set 0 values to non-0, etc. */
242
243 int j, notallzero = 0;
244
245 /* Set 0 to 1 :
246 for(j = 0; j <= RNG_Table[RNG_kind].n_seed - 1; j++)
247 if(!RNG_Table[RNG_kind].i_seed[j]) RNG_Table[RNG_kind].i_seed[j]++; */
248
249 switch(RNG_kind) {
250 case WICHMANN_HILL:
251 I1 = I1 % 30269; I2 = I2 % 30307; I3 = I3 % 30323;
252
253 /* map values equal to 0 mod modulus to 1. */
254 if(I1 == 0) I1 = 1;
255 if(I2 == 0) I2 = 1;
256 if(I3 == 0) I3 = 1;
257 return;
258
259 case SUPER_DUPER:
260 if(I1 == 0) I1 = 1;
261 /* I2 = Congruential: must be ODD */
262 I2 |= 1;
263 break;
264
265 case MARSAGLIA_MULTICARRY:
266 if(I1 == 0) I1 = 1;
267 if(I2 == 0) I2 = 1;
268 break;
269
270 case MERSENNE_TWISTER:
271 if(initial) I1 = 624;
272 /* No action unless user has corrupted .Random.seed */
273 if(I1 <= 0) I1 = 624;
274 /* check for all zeroes */
275 for (j = 1; j <= 624; j++)
276 if(RNG_Table[RNG_kind].i_seed[j] != 0) {
277 notallzero = 1;
278 break;
279 }
280 if(!notallzero) Randomize(RNG_kind);
281 break;
282
283 case KNUTH_TAOCP:
284 case KNUTH_TAOCP2:
285 if(KT_pos <= 0) KT_pos = 100;
286 /* check for all zeroes */
287 for (j = 0; j < 100; j++)
288 if(RNG_Table[RNG_kind].i_seed[j] != 0) {
289 notallzero = 1;
290 break;
291 }
292 if(!notallzero) Randomize(RNG_kind);
293 break;
294 case USER_UNIF:
295 break;
296 default:
297 error(_("FixupSeeds: unimplemented RNG kind %d"), RNG_kind);
298 }
299 }
300
301 #ifdef GNU_R_MODE
302 extern double BM_norm_keep; /* ../nmath/snorm.c */
303 #endif
304
RNG_Init(RNGtype kind,Int32 seed)305 static void RNG_Init(RNGtype kind, Int32 seed)
306 {
307 int j;
308
309 #ifdef GNU_R_MODE
310 BM_norm_keep = 0.0; /* zap Box-Muller history */
311 #endif
312
313 /* Initial scrambling */
314 for(j = 0; j < 50; j++)
315 seed = (69069 * seed + 1);
316 switch(kind) {
317 case WICHMANN_HILL:
318 case MARSAGLIA_MULTICARRY:
319 case SUPER_DUPER:
320 case MERSENNE_TWISTER:
321 /* i_seed[0] is mti, *but* this is needed for historical consistency */
322 for(j = 0; j < RNG_Table[kind].n_seed; j++) {
323 seed = (69069 * seed + 1);
324 RNG_Table[kind].i_seed[j] = seed;
325 }
326 FixupSeeds(kind, 1);
327 break;
328 case KNUTH_TAOCP:
329 RNG_Init_R_KT(seed);
330 break;
331 case KNUTH_TAOCP2:
332 RNG_Init_KT2(seed);
333 break;
334 case USER_UNIF:
335 #if GNU_R_MODE
336 User_unif_fun = R_FindSymbol("user_unif_rand", "", NULL);
337 if (!User_unif_fun) error(_("'user_unif_rand' not in load table"));
338 User_unif_init = R_FindSymbol("user_unif_init", "", NULL);
339 if (User_unif_init) (void) User_unif_init(seed);
340 User_unif_nseed = R_FindSymbol("user_unif_nseed", "", NULL);
341 User_unif_seedloc = R_FindSymbol("user_unif_seedloc", "", NULL);
342 if (User_unif_seedloc) {
343 int ns = 0;
344 if (!User_unif_nseed) {
345 warning(_("cannot read seeds unless 'user_unif_nseed' is supplied"));
346 break;
347 }
348 ns = *((int *) User_unif_nseed());
349 if (ns < 0 || ns > 625) {
350 warning(_("seed length must be in 0...625; ignored"));
351 break;
352 }
353 RNG_Table[kind].n_seed = ns;
354 RNG_Table[kind].i_seed = (Int32 *) User_unif_seedloc();
355 }
356 break;
357 #endif
358 default:
359 error(_("RNG_Init: unimplemented RNG kind %d"), kind);
360 }
361 }
362
363 #include <time.h>
Randomize(RNGtype kind)364 static void Randomize(RNGtype kind)
365 {
366 /* Only called by GetRNGstate() when there's no .Random.seed */
367
368 RNG_Init(kind, (Int32) time(NULL));
369 }
370
371 /* function needed by dieharder follow below */
372 #ifdef DIEHARDER
373
374 /* ------------------------------------------------ GNU R generic getters */
375 static unsigned long int
r_rng_get(void * vstate)376 r_rng_get (void *vstate)
377 {
378 unsigned long int j = r_rng_get_double(vstate) * RANDNUM_MAX;
379 return(j);
380 }
381
382 static double
r_rng_get_double(void * vstate)383 r_rng_get_double (void *vstate)
384 {
385 r_rng_tab_t *state = (r_rng_tab_t *) vstate;
386 RNG_kind = state->kind;
387 double u = unif_rand();
388 return(u);
389 }
390
391 /* -------------------------------------- GNU R generator 1: wichmann_hill */
392 static void
r_wichmann_hill_set(void * vstate,unsigned long int s)393 r_wichmann_hill_set (void *vstate, unsigned long int s)
394 {
395 r_rng_tab_t *state = (r_rng_tab_t *) vstate;
396 *state = RNG_Table[WICHMANN_HILL];
397 RNGtype kind = state->kind;
398 Int32 seed = s; /* FIXME: casting down from long seed to unsigned int */
399 RNG_Init(kind, seed);
400 return;
401 }
402
403 static const gsl_rng_type r_wichmann_hill_type =
404 {"R_wichmann_hill", /* name */
405 RANDNUM_MAX, /* RAND_MAX */
406 0, /* RAND_MIN */
407 sizeof (r_rng_tab_t),
408 &r_wichmann_hill_set,
409 &r_rng_get,
410 &r_rng_get_double};
411
412 const gsl_rng_type *gsl_rng_r_wichmann_hill = &r_wichmann_hill_type;
413
414 /* -------------------------------------- GNU R generator 2: marsaglia_mc */
415 static void
r_marsaglia_mc_set(void * vstate,unsigned long int s)416 r_marsaglia_mc_set (void *vstate, unsigned long int s)
417 {
418 r_rng_tab_t *state = (r_rng_tab_t *) vstate;
419 *state = RNG_Table[MARSAGLIA_MULTICARRY];
420 RNGtype kind = state->kind;
421 Int32 seed = s; /* FIXME: casting down from long seed to unsigned int */
422 RNG_Init(kind, seed);
423 return;
424 }
425
426 static const gsl_rng_type r_marsaglia_mc_type =
427 {"R_marsaglia_multic.", /* name */
428 RANDNUM_MAX, /* RAND_MAX */
429 0, /* RAND_MIN */
430 sizeof (r_rng_tab_t),
431 &r_marsaglia_mc_set,
432 &r_rng_get,
433 &r_rng_get_double};
434
435 const gsl_rng_type *gsl_rng_r_marsaglia_mc = &r_marsaglia_mc_type;
436
437 /* -------------------------------------- GNU R generator 3: super_duper */
438 static void
r_super_duper_set(void * vstate,unsigned long int s)439 r_super_duper_set (void *vstate, unsigned long int s)
440 {
441 r_rng_tab_t *state = (r_rng_tab_t *) vstate;
442 *state = RNG_Table[SUPER_DUPER];
443 RNGtype kind = state->kind;
444 Int32 seed = s; /* FIXME: casting down from long seed to unsigned int */
445 RNG_Init(kind, seed);
446 return;
447 }
448
449 static const gsl_rng_type r_super_duper_type =
450 {"R_super_duper", /* name */
451 RANDNUM_MAX, /* RAND_MAX */
452 0, /* RAND_MIN */
453 sizeof (r_rng_tab_t),
454 &r_super_duper_set,
455 &r_rng_get,
456 &r_rng_get_double};
457
458 const gsl_rng_type *gsl_rng_r_super_duper = &r_super_duper_type;
459
460 /* ------------------------------------- GNU R generator 4: mersenne_twister */
461 static void
r_mersenne_twister_set(void * vstate,unsigned long int s)462 r_mersenne_twister_set (void *vstate, unsigned long int s)
463 {
464 r_rng_tab_t *state = (r_rng_tab_t *) vstate;
465 *state = RNG_Table[MERSENNE_TWISTER];
466 RNGtype kind = state->kind;
467 Int32 seed = s; /* FIXME: casting down from long seed to unsigned int */
468 RNG_Init(kind, seed);
469 return;
470 }
471
472 static const gsl_rng_type r_mersenne_twister_type =
473 {"R_mersenne_twister", /* name */
474 RANDNUM_MAX, /* RAND_MAX */
475 0, /* RAND_MIN */
476 sizeof (r_rng_tab_t),
477 &r_mersenne_twister_set,
478 &r_rng_get,
479 &r_rng_get_double};
480
481 const gsl_rng_type *gsl_rng_r_mersenne_twister = &r_mersenne_twister_type;
482
483 /* ------------------------------------- GNU R generator 5: knuth_taocp */
484 static void
r_knuth_taocp_set(void * vstate,unsigned long int s)485 r_knuth_taocp_set (void *vstate, unsigned long int s)
486 {
487 r_rng_tab_t *state = (r_rng_tab_t *) vstate;
488 *state = RNG_Table[KNUTH_TAOCP];
489 RNGtype kind = state->kind;
490 Int32 seed = s; /* FIXME: casting down from long seed to unsigned int */
491 RNG_Init(kind, seed);
492 return;
493 }
494
495 static unsigned long int
r_rng_get_taocp(void * vstate)496 r_rng_get_taocp (void *vstate)
497 {
498 unsigned long int j = r_rng_get_double(vstate) * TAOCP_MAX;
499 return(j);
500 }
501
502 static const gsl_rng_type r_knuth_taocp_type =
503 {"R_knuth_taocp", /* name */
504 TAOCP_MAX, /* RAND_MAX */
505 0, /* RAND_MIN */
506 sizeof (r_rng_tab_t),
507 &r_knuth_taocp_set,
508 &r_rng_get_taocp,
509 &r_rng_get_double};
510
511 const gsl_rng_type *gsl_rng_r_knuth_taocp = &r_knuth_taocp_type;
512
513 /* ------------------------------------- GNU R generator 6: knuth_taocp2 */
514 static void
r_knuth_taocp2_set(void * vstate,unsigned long int s)515 r_knuth_taocp2_set (void *vstate, unsigned long int s)
516 {
517 r_rng_tab_t *state = (r_rng_tab_t *) vstate;
518 *state = RNG_Table[KNUTH_TAOCP2];
519 RNGtype kind = state->kind;
520 Int32 seed = s; /* FIXME: casting down from long seed to unsigned int */
521 RNG_Init(kind, seed);
522 return;
523 }
524
525 static unsigned long int
r_rng_get_taocp2(void * vstate)526 r_rng_get_taocp2 (void *vstate)
527 {
528 unsigned long int j = r_rng_get_double(vstate) * TAOCP2_MAX;
529 return(j);
530 }
531
532 static const gsl_rng_type r_knuth_taocp2_type =
533 {"R_knuth_taocp2", /* name */
534 TAOCP2_MAX, /* RAND_MAX */
535 0, /* RAND_MIN */
536 sizeof (r_rng_tab_t),
537 &r_knuth_taocp2_set,
538 &r_rng_get_taocp2,
539 &r_rng_get_double};
540
541 const gsl_rng_type *gsl_rng_r_knuth_taocp2 = &r_knuth_taocp2_type;
542
543 #endif
544
545 /* R internal functions below before original implementations follow */
546 #if GNU_R_MODE
GetRNGstate()547 void GetRNGstate()
548 {
549 /* Get .Random.seed into proper variables */
550 int len_seed, j, tmp;
551 SEXP seeds;
552 RNGtype newRNG; N01type newN01;
553
554 /* look only in the workspace */
555 seeds = findVarInFrame(R_GlobalEnv, R_SeedsSymbol);
556 if (seeds == R_UnboundValue) {
557 Randomize(RNG_kind);
558 }
559 else {
560 seeds = coerceVector(seeds, INTSXP);
561 if (seeds == R_MissingArg)
562 error(_(".Random.seed is a missing argument with no default"));
563 if (!isVector(seeds))
564 error(_(".Random.seed is not a vector"));
565 tmp = INTEGER(seeds)[0];
566 if (tmp == NA_INTEGER)
567 error(_(".Random.seed[1] is not a valid integer"));
568 newRNG = tmp % 100;
569 newN01 = tmp / 100;
570 /*if (RNG_kind > USER_UNIF || RNG_kind < 0) {
571 warning(".Random.seed was invalid: re-initializing");
572 RNG_kind = RNG_DEFAULT;
573 }*/
574 if (newN01 < 0 || newN01 > KINDERMAN_RAMAGE)
575 error(_(".Random.seed[0] is not a valid Normal type"));
576 switch(newRNG) {
577 case WICHMANN_HILL:
578 case MARSAGLIA_MULTICARRY:
579 case SUPER_DUPER:
580 case MERSENNE_TWISTER:
581 case KNUTH_TAOCP:
582 case KNUTH_TAOCP2:
583 break;
584 case USER_UNIF:
585 if(!User_unif_fun)
586 error(_(".Random.seed[1] = 5 but no user-supplied generator"));
587 break;
588 default:
589 error(_(".Random.seed[1] is not a valid RNG kind (code)"));
590 }
591 RNG_kind = newRNG; N01_kind = newN01;
592 len_seed = RNG_Table[RNG_kind].n_seed;
593 /* Not sure whether this test is needed: wrong for USER_UNIF */
594 if(LENGTH(seeds) > 1 && LENGTH(seeds) < len_seed + 1)
595 error(_(".Random.seed has wrong length"));
596 if(LENGTH(seeds) == 1 && RNG_kind != USER_UNIF)
597 Randomize(RNG_kind);
598 else {
599 for(j = 1; j <= len_seed; j++) {
600 tmp = INTEGER(seeds)[j];
601 /* Some generators can generate NA_INTEGER as a valid integer value */
602 /* if(tmp == NA_INTEGER)
603 error(".Random.seed[%d] is not a valid integer", j+1);*/
604 RNG_Table[RNG_kind].i_seed[j - 1] = tmp;
605 }
606 FixupSeeds(RNG_kind, 0);
607 }
608 }
609 }
610
PutRNGstate()611 void PutRNGstate()
612 {
613 /* Copy out seeds to .Random.seed */
614 int len_seed, j;
615 SEXP seeds;
616
617 if (RNG_kind < 0 || RNG_kind > KNUTH_TAOCP2 ||
618 N01_kind < 0 || N01_kind > KINDERMAN_RAMAGE) {
619 warning("Internal .Random.seed is corrupt: not saving");
620 return;
621 }
622
623 len_seed = RNG_Table[RNG_kind].n_seed;
624
625 PROTECT(seeds = allocVector(INTSXP, len_seed + 1));
626
627 INTEGER(seeds)[0] = RNG_kind + 100 * N01_kind;
628 for(j = 0; j < len_seed; j++)
629 INTEGER(seeds)[j+1] = RNG_Table[RNG_kind].i_seed[j];
630
631 /* assign only in the workspace */
632 defineVar(R_SeedsSymbol, seeds, R_GlobalEnv);
633 UNPROTECT(1);
634 }
635
RNGkind(RNGtype newkind)636 static void RNGkind(RNGtype newkind)
637 {
638 /* Choose a new kind of RNG.
639 * Initialize its seed by calling the old RNG's unif_rand()
640 */
641 if (newkind == -1) newkind = RNG_DEFAULT;
642 switch(newkind) {
643 case WICHMANN_HILL:
644 case MARSAGLIA_MULTICARRY:
645 case SUPER_DUPER:
646 case MERSENNE_TWISTER:
647 case KNUTH_TAOCP:
648 case KNUTH_TAOCP2:
649 case USER_UNIF:
650 break;
651 default:
652 error(_("RNGkind: unimplemented RNG kind %d"), newkind);
653 }
654
655 GetRNGstate();
656 RNG_Init(newkind, unif_rand() * UINT_MAX);
657 RNG_kind = newkind;
658 PutRNGstate();
659 }
660
Norm_kind(N01type kind)661 static void Norm_kind(N01type kind)
662 {
663 if (kind == -1) kind = N01_DEFAULT;
664 if (kind < 0 || kind > KINDERMAN_RAMAGE)
665 error(_("invalid Normal type in RNGkind"));
666 if (kind == USER_NORM) {
667 User_norm_fun = R_FindSymbol("user_norm_rand", "", NULL);
668 if (!User_norm_fun) error(_("'user_norm_rand' not in load table"));
669 }
670 GetRNGstate(); /* might not be initialized */
671 if (kind == BOX_MULLER)
672 BM_norm_keep = 0.0; /* zap Box-Muller history */
673 N01_kind = kind;
674 PutRNGstate();
675 }
676
677 /*------ .Internal interface ------------------------*/
678
do_RNGkind(SEXP call,SEXP op,SEXP args,SEXP env)679 SEXP attribute_hidden do_RNGkind (SEXP call, SEXP op, SEXP args, SEXP env)
680 {
681 SEXP ans, rng, norm;
682
683 checkArity(op,args);
684 PROTECT(ans = allocVector(INTSXP, 2));
685 INTEGER(ans)[0] = RNG_kind;
686 INTEGER(ans)[1] = N01_kind;
687 rng = CAR(args);
688 norm = CADR(args);
689 if(!isNull(rng)) { /* set a new RNG kind */
690 RNGkind(asInteger(rng));
691 }
692 if(!isNull(norm)) { /* set a new normal kind */
693 Norm_kind(asInteger(norm));
694 }
695 UNPROTECT(1);
696 return ans;
697 }
698
699
do_setseed(SEXP call,SEXP op,SEXP args,SEXP env)700 SEXP attribute_hidden do_setseed (SEXP call, SEXP op, SEXP args, SEXP env)
701 {
702 SEXP skind;
703 int seed;
704 RNGtype kind;
705
706 checkArity(op,args);
707 seed = asInteger(CAR(args));
708 if (seed == NA_INTEGER)
709 error(_("supplied seed is not a valid integer"));
710 skind = CADR(args);
711 if (!isNull(skind)) {
712 kind = asInteger(skind);
713 RNGkind(kind);
714 } else
715 kind = RNG_kind;
716 RNG_Init(RNG_kind, (Int32) seed);
717 PutRNGstate();
718 return R_NilValue;
719 }
720
721
722 /* S COMPATIBILITY */
723
724 /* The following entry points provide compatibility with S. */
725 /* These entry points should not be used by new R code. */
726
seed_in(long * ignored)727 void seed_in(long *ignored)
728 {
729 GetRNGstate();
730 }
731
seed_out(long * ignored)732 void seed_out(long *ignored)
733 {
734 PutRNGstate();
735 }
736 #endif
737
738 /* =================== Mersenne Twister ========================== */
739 /* From http://www.math.keio.ac.jp/~matumoto/emt.html */
740
741 /* A C-program for MT19937: Real number version([0,1)-interval)
742 (1999/10/28)
743 genrand() generates one pseudorandom real number (double)
744 which is uniformly distributed on [0,1)-interval, for each
745 call. sgenrand(seed) sets initial values to the working area
746 of 624 words. Before genrand(), sgenrand(seed) must be
747 called once. (seed is any 32-bit integer.)
748 Integer generator is obtained by modifying two lines.
749 Coded by Takuji Nishimura, considering the suggestions by
750 Topher Cooper and Marc Rieffel in July-Aug. 1997.
751
752 Copyright (C) 1997, 1999 Makoto Matsumoto and Takuji Nishimura.
753 When you use this, send an email to: matumoto@math.keio.ac.jp
754 with an appropriate reference to your work.
755
756 REFERENCE
757 M. Matsumoto and T. Nishimura,
758 "Mersenne Twister: A 623-Dimensionally Equidistributed Uniform
759 Pseudo-Random Number Generator",
760 ACM Transactions on Modeling and Computer Simulation,
761 Vol. 8, No. 1, January 1998, pp 3--30.
762 */
763
764 /* Period parameters */
765 #define N 624
766 #define M 397
767 #define MATRIX_A 0x9908b0df /* constant vector a */
768 #define UPPER_MASK 0x80000000 /* most significant w-r bits */
769 #define LOWER_MASK 0x7fffffff /* least significant r bits */
770
771 /* Tempering parameters */
772 #define TEMPERING_MASK_B 0x9d2c5680
773 #define TEMPERING_MASK_C 0xefc60000
774 #define TEMPERING_SHIFT_U(y) (y >> 11)
775 #define TEMPERING_SHIFT_S(y) (y << 7)
776 #define TEMPERING_SHIFT_T(y) (y << 15)
777 #define TEMPERING_SHIFT_L(y) (y >> 18)
778
779 static Int32 *mt = dummyvec+1; /* the array for the state vector */
780 static int mti=N+1; /* mti==N+1 means mt[N] is not initialized */
781
782 /* Initializing the array with a seed */
783 static void
MT_sgenrand(Int32 seed)784 MT_sgenrand(Int32 seed)
785 {
786 int i;
787
788 for (i = 0; i < N; i++) {
789 mt[i] = seed & 0xffff0000;
790 seed = 69069 * seed + 1;
791 mt[i] |= (seed & 0xffff0000) >> 16;
792 seed = 69069 * seed + 1;
793 }
794 mti = N;
795 }
796
797 /* Initialization by "sgenrand()" is an example. Theoretically,
798 there are 2^19937-1 possible states as an intial state.
799 Essential bits in "seed_array[]" is following 19937 bits:
800 (seed_array[0]&UPPER_MASK), seed_array[1], ..., seed_array[N-1].
801 (seed_array[0]&LOWER_MASK) is discarded.
802 Theoretically,
803 (seed_array[0]&UPPER_MASK), seed_array[1], ..., seed_array[N-1]
804 can take any values except all zeros. */
805
MT_genrand(void)806 static double MT_genrand(void)
807 {
808 Int32 y;
809 static Int32 mag01[2]={0x0, MATRIX_A};
810 /* mag01[x] = x * MATRIX_A for x=0,1 */
811
812 mti = dummyvec[0];
813
814 if (mti >= N) { /* generate N words at one time */
815 int kk;
816
817 if (mti == N+1) /* if sgenrand() has not been called, */
818 MT_sgenrand(4357); /* a default initial seed is used */
819
820 for (kk = 0; kk < N - M; kk++) {
821 y = (mt[kk] & UPPER_MASK) | (mt[kk+1] & LOWER_MASK);
822 mt[kk] = mt[kk+M] ^ (y >> 1) ^ mag01[y & 0x1];
823 }
824 for (; kk < N - 1; kk++) {
825 y = (mt[kk] & UPPER_MASK) | (mt[kk+1] & LOWER_MASK);
826 mt[kk] = mt[kk+(M-N)] ^ (y >> 1) ^ mag01[y & 0x1];
827 }
828 y = (mt[N-1] & UPPER_MASK) | (mt[0] & LOWER_MASK);
829 mt[N-1] = mt[M-1] ^ (y >> 1) ^ mag01[y & 0x1];
830
831 mti = 0;
832 }
833
834 y = mt[mti++];
835 y ^= TEMPERING_SHIFT_U(y);
836 y ^= TEMPERING_SHIFT_S(y) & TEMPERING_MASK_B;
837 y ^= TEMPERING_SHIFT_T(y) & TEMPERING_MASK_C;
838 y ^= TEMPERING_SHIFT_L(y);
839 dummyvec[0] = mti;
840
841 return ( (double)y * 2.3283064365386963e-10 ); /* reals: [0,1)-interval */
842 }
843
844 /*
845 The following code was taken from earlier versions of
846 http://www-cs-faculty.stanford.edu/~knuth/programs/rng.c-old
847 http://www-cs-faculty.stanford.edu/~knuth/programs/rng.c
848 */
849
850
851 #define long Int32
852 #define ran_arr_buf R_KT_ran_arr_buf
853 #define ran_arr_cycle R_KT_ran_arr_cycle
854 #define ran_arr_ptr R_KT_ran_arr_ptr
855 #define ran_arr_sentinel R_KT_ran_arr_sentinel
856 #define ran_x dummyvec
857
858 #define KK 100 /* the long lag */
859 #define LL 37 /* the short lag */
860 #define MM (1L<<30) /* the modulus */
861 #define TT 70 /* guaranteed separation between streams */
862 #define mod_diff(x,y) (((x)-(y))&(MM-1)) /* subtraction mod MM */
863 #define is_odd(x) ((x)&1) /* units bit of x */
ran_array(long aa[],int n)864 static void ran_array(long aa[],int n) /* put n new random numbers in aa */
865 {
866 register int i,j;
867 for (j=0;j<KK;j++) aa[j]=ran_x[j];
868 for (;j<n;j++) aa[j]=mod_diff(aa[j-KK],aa[j-LL]);
869 for (i=0;i<LL;i++,j++) ran_x[i]=mod_diff(aa[j-KK],aa[j-LL]);
870 for (;i<KK;i++,j++) ran_x[i]=mod_diff(aa[j-KK],ran_x[i-LL]);
871 }
872 #define QUALITY 1009 /* recommended quality level for high-res use */
873 static long ran_arr_buf[QUALITY];
874 static long ran_arr_sentinel=(long)-1;
875 static long *ran_arr_ptr=&ran_arr_sentinel; /* the next random number, or -1 */
876
ran_arr_cycle(void)877 static long ran_arr_cycle(void)
878 {
879 ran_array(ran_arr_buf,QUALITY);
880 ran_arr_buf[KK]=-1;
881 ran_arr_ptr=ran_arr_buf+1;
882 return ran_arr_buf[0];
883 }
884
885 /* =================== Knuth TAOCP 2002 ========================== */
886
887 /* This program by D E Knuth is in the public domain and freely copyable.
888 * It is explained in Seminumerical Algorithms, 3rd edition, Section 3.6
889 * (or in the errata to the 2nd edition --- see
890 * http://www-cs-faculty.stanford.edu/~knuth/taocp.html
891 * in the changes to Volume 2 on pages 171 and following). */
892
893 /* N.B. The MODIFICATIONS introduced in the 9th printing (2002) are
894 included here; there's no backwards compatibility with the original. */
895
896
897 #ifdef __STDC__
ran_start(long seed)898 void ran_start(long seed)
899 #else
900 void ran_start(seed) /* do this before using ran_array */
901 long seed; /* selector for different streams */
902 #endif
903 {
904 register int t,j;
905 long x[KK+KK-1]; /* the preparation buffer */
906 register long ss=(seed+2)&(MM-2);
907 for (j=0;j<KK;j++) {
908 x[j]=ss; /* bootstrap the buffer */
909 ss<<=1; if (ss>=MM) ss-=MM-2; /* cyclic shift 29 bits */
910 }
911 x[1]++; /* make x[1] (and only x[1]) odd */
912 for (ss=seed&(MM-1),t=TT-1; t; ) {
913 for (j=KK-1;j>0;j--) x[j+j]=x[j], x[j+j-1]=0; /* "square" */
914 for (j=KK+KK-2;j>=KK;j--)
915 x[j-(KK-LL)]=mod_diff(x[j-(KK-LL)],x[j]),
916 x[j-KK]=mod_diff(x[j-KK],x[j]);
917 if (is_odd(ss)) { /* "multiply by z" */
918 for (j=KK;j>0;j--) x[j]=x[j-1];
919 x[0]=x[KK]; /* shift the buffer cyclically */
920 x[LL]=mod_diff(x[LL],x[KK]);
921 }
922 if (ss) ss>>=1; else t--;
923 }
924 for (j=0;j<LL;j++) ran_x[j+KK-LL]=x[j];
925 for (;j<KK;j++) ran_x[j-LL]=x[j];
926 for (j=0;j<10;j++) ran_array(x,KK+KK-1); /* warm things up */
927 ran_arr_ptr=&ran_arr_sentinel;
928 }
929 /* ===================== end of Knuth's code ====================== */
930
RNG_Init_KT2(Int32 seed)931 static void RNG_Init_KT2(Int32 seed)
932 {
933 ran_start(seed % 1073741821);
934 KT_pos = 100;
935 }
936
KT_next(void)937 static Int32 KT_next(void)
938 {
939 if(KT_pos >= 100) {
940 ran_arr_cycle();
941 KT_pos = 0;
942 }
943 return ran_x[(KT_pos)++];
944 }
945
RNG_Init_R_KT(Int32 seed)946 static void RNG_Init_R_KT(Int32 seed)
947 {
948 ran_start(seed % 1073741821);
949 KT_pos = 100;
950 }
951