| /* |
| * R : A Computer Language for Statistical Data Analysis |
| * Copyright (C) 1995, 1996 Robert Gentleman and Ross Ihaka |
| * Copyright (C) 1997--2019 The R Core Team |
| * |
| * This program is free software; you can redistribute it and/or modify |
| * it under the terms of the GNU General Public License as published by |
| * the Free Software Foundation; either version 2 of the License, or |
| * (at your option) any later version. |
| * |
| * This program is distributed in the hope that it will be useful, |
| * but WITHOUT ANY WARRANTY; without even the implied warranty of |
| * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the |
| * GNU General Public License for more details. |
| * |
| * You should have received a copy of the GNU General Public License |
| * along with this program; if not, a copy is available at |
| * https://www.R-project.org/Licenses/ |
| */ |
| |
| |
| #ifdef HAVE_CONFIG_H |
| #include <config.h> |
| #endif |
| |
| #include <Defn.h> |
| #include <Internal.h> |
| #include <R_ext/Random.h> |
| |
| /* Normal generator is not actually set here but in ../nmath/snorm.c */ |
| #define RNG_DEFAULT MERSENNE_TWISTER |
| #define N01_DEFAULT INVERSION |
| #define Sample_DEFAULT REJECTION |
| |
| |
| #include <R_ext/Rdynload.h> |
| |
| static DL_FUNC User_unif_fun, User_unif_nseed, |
| User_unif_seedloc; |
| typedef void (*UnifInitFun)(Int32); |
| |
| UnifInitFun User_unif_init = NULL; /* some picky compilers */ |
| |
| DL_FUNC User_norm_fun = NULL; /* also in ../nmath/snorm.c */ |
| |
| #include "nmath2.h" |
| static RNGtype RNG_kind = RNG_DEFAULT; |
| //extern N01type N01_kind; /* from ../nmath/snorm.c */ |
| //extern double BM_norm_keep; /* ../nmath/snorm.c */ |
| static Sampletype Sample_kind = REJECTION; |
| |
| /* typedef unsigned int Int32; in Random.h */ |
| |
| /* .Random.seed == (RNGkind, i_seed[0],i_seed[1],..,i_seed[n_seed-1]) |
| * or == (RNGkind) or missing [--> Randomize] |
| * where RNGkind := RNG_kind + 100 * N01_kind + 10000 * Sample_kind |
| * currently in outer(outer(0:7, 100*(0:5), "+"), 10000*(0:1), "+") |
| */ |
| |
| typedef struct { |
| RNGtype kind; |
| N01type Nkind; |
| char *name; /* print name */ |
| int n_seed; /* length of seed vector */ |
| Int32 *i_seed; |
| } RNGTAB; |
| |
| |
| static Int32 dummy[625]; |
| static |
| RNGTAB RNG_Table[] = |
| { |
| /* kind Nkind name n_seed i_seed */ |
| { WICHMANN_HILL, BUGGY_KINDERMAN_RAMAGE, "Wichmann-Hill", 3, dummy}, |
| { MARSAGLIA_MULTICARRY, BUGGY_KINDERMAN_RAMAGE, "Marsaglia-MultiCarry", 2, dummy}, |
| { SUPER_DUPER, BUGGY_KINDERMAN_RAMAGE, "Super-Duper", 2, dummy}, |
| { MERSENNE_TWISTER, BUGGY_KINDERMAN_RAMAGE, "Mersenne-Twister", 1+624, dummy}, |
| { KNUTH_TAOCP, BUGGY_KINDERMAN_RAMAGE, "Knuth-TAOCP", 1+100, dummy}, |
| { USER_UNIF, BUGGY_KINDERMAN_RAMAGE, "User-supplied", 0, dummy}, |
| { KNUTH_TAOCP2, BUGGY_KINDERMAN_RAMAGE, "Knuth-TAOCP-2002", 1+100, dummy}, |
| { LECUYER_CMRG, BUGGY_KINDERMAN_RAMAGE, "L'Ecuyer-CMRG", 6, dummy}, |
| }; |
| |
| |
| #define d2_32 4294967296./* = (double) */ |
| #define i2_32m1 2.328306437080797e-10/* = 1/(2^32 - 1) */ |
| #define KT 9.31322574615479e-10 /* = 2^-30 */ |
| |
| #define I1 (RNG_Table[RNG_kind].i_seed[0]) |
| #define I2 (RNG_Table[RNG_kind].i_seed[1]) |
| #define I3 (RNG_Table[RNG_kind].i_seed[2]) |
| |
| static void Randomize(RNGtype kind); |
| static double MT_genrand(void); |
| static Int32 KT_next(void); |
| static void RNG_Init_R_KT(Int32); |
| static void RNG_Init_KT2(Int32); |
| #define KT_pos (RNG_Table[KNUTH_TAOCP].i_seed[100]) |
| |
| static double fixup(double x) |
| { |
| /* ensure 0 and 1 are never returned */ |
| if(x <= 0.0) return 0.5*i2_32m1; |
| if((1.0 - x) <= 0.0) return 1.0 - 0.5*i2_32m1; |
| return x; |
| } |
| |
| |
| double unif_rand(void) |
| { |
| double value; |
| |
| switch(RNG_kind) { |
| |
| case WICHMANN_HILL: |
| I1 = I1 * 171 % 30269; |
| I2 = I2 * 172 % 30307; |
| I3 = I3 * 170 % 30323; |
| value = I1 / 30269.0 + I2 / 30307.0 + I3 / 30323.0; |
| return fixup(value - (int) value);/* in [0,1) */ |
| |
| case MARSAGLIA_MULTICARRY:/* 0177777(octal) == 65535(decimal)*/ |
| I1= 36969*(I1 & 0177777) + (I1>>16); |
| I2= 18000*(I2 & 0177777) + (I2>>16); |
| return fixup(((I1 << 16)^(I2 & 0177777)) * i2_32m1); /* in [0,1) */ |
| |
| case SUPER_DUPER: |
| /* This is Reeds et al (1984) implementation; |
| * modified using __unsigned__ seeds instead of signed ones |
| */ |
| I1 ^= ((I1 >> 15) & 0377777); /* Tausworthe */ |
| I1 ^= I1 << 17; |
| I2 *= 69069; /* Congruential */ |
| return fixup((I1^I2) * i2_32m1); /* in [0,1) */ |
| |
| case MERSENNE_TWISTER: |
| return fixup(MT_genrand()); |
| |
| case KNUTH_TAOCP: |
| case KNUTH_TAOCP2: |
| return fixup(KT_next() * KT); |
| |
| case USER_UNIF: |
| return *((double *) User_unif_fun()); |
| |
| case LECUYER_CMRG: |
| { |
| /* Based loosely on the GPL-ed version of |
| http://www.iro.umontreal.ca/~lecuyer/myftp/streams00/c2010/RngStream.c |
| but using int_least64_t, which C99 guarantees. |
| */ |
| int k; |
| int_least64_t p1, p2; |
| |
| #define II(i) (RNG_Table[RNG_kind].i_seed[i]) |
| #define m1 4294967087 |
| #define m2 4294944443 |
| #define normc 2.328306549295727688e-10 |
| #define a12 (int_least64_t)1403580 |
| #define a13n (int_least64_t)810728 |
| #define a21 (int_least64_t)527612 |
| #define a23n (int_least64_t)1370589 |
| |
| p1 = a12 * (unsigned int)II(1) - a13n * (unsigned int)II(0); |
| /* p1 % m1 would surely do */ |
| k = (int) (p1 / m1); |
| p1 -= k * m1; |
| if (p1 < 0.0) p1 += m1; |
| II(0) = II(1); II(1) = II(2); II(2) = (int) p1; |
| |
| p2 = a21 * (unsigned int)II(5) - a23n * (unsigned int)II(3); |
| k = (int) (p2 / m2); |
| p2 -= k * m2; |
| if (p2 < 0.0) p2 += m2; |
| II(3) = II(4); II(4) = II(5); II(5) = (int) p2; |
| |
| return (double)((p1 > p2) ? (p1 - p2) : (p1 - p2 + m1)) * normc; |
| } |
| default: |
| error(_("unif_rand: unimplemented RNG kind %d"), RNG_kind); |
| return -1.; |
| } |
| } |
| |
| /* we must mask global variable here, as I1-I3 hide RNG_kind |
| and we want the argument */ |
| static void FixupSeeds(RNGtype RNG_kind, int initial) |
| { |
| /* Depending on RNG, set 0 values to non-0, etc. */ |
| |
| int j, notallzero = 0; |
| |
| /* Set 0 to 1 : |
| for(j = 0; j <= RNG_Table[RNG_kind].n_seed - 1; j++) |
| if(!RNG_Table[RNG_kind].i_seed[j]) RNG_Table[RNG_kind].i_seed[j]++; */ |
| |
| switch(RNG_kind) { |
| case WICHMANN_HILL: |
| I1 = I1 % 30269; I2 = I2 % 30307; I3 = I3 % 30323; |
| |
| /* map values equal to 0 mod modulus to 1. */ |
| if(I1 == 0) I1 = 1; |
| if(I2 == 0) I2 = 1; |
| if(I3 == 0) I3 = 1; |
| return; |
| |
| case SUPER_DUPER: |
| if(I1 == 0) I1 = 1; |
| /* I2 = Congruential: must be ODD */ |
| I2 |= 1; |
| break; |
| |
| case MARSAGLIA_MULTICARRY: |
| if(I1 == 0) I1 = 1; |
| if(I2 == 0) I2 = 1; |
| break; |
| |
| case MERSENNE_TWISTER: |
| if(initial) I1 = 624; |
| /* No action unless user has corrupted .Random.seed */ |
| if(I1 <= 0) I1 = 624; |
| /* check for all zeroes */ |
| for (j = 1; j <= 624; j++) |
| if(RNG_Table[RNG_kind].i_seed[j] != 0) { |
| notallzero = 1; |
| break; |
| } |
| if(!notallzero) Randomize(RNG_kind); |
| break; |
| |
| case KNUTH_TAOCP: |
| case KNUTH_TAOCP2: |
| if(KT_pos <= 0) KT_pos = 100; |
| /* check for all zeroes */ |
| for (j = 0; j < 100; j++) |
| if(RNG_Table[RNG_kind].i_seed[j] != 0) { |
| notallzero = 1; |
| break; |
| } |
| if(!notallzero) Randomize(RNG_kind); |
| break; |
| case USER_UNIF: |
| break; |
| case LECUYER_CMRG: |
| /* first set: not all zero, in [0, m1) |
| second set: not all zero, in [0, m2) */ |
| { |
| unsigned int tmp; |
| int allOK = 1; |
| for (j = 0; j < 3; j++) { |
| tmp = RNG_Table[RNG_kind].i_seed[j]; |
| if(tmp != 0) notallzero = 1; |
| if (tmp >= m1) allOK = 0; |
| } |
| if(!notallzero || !allOK) Randomize(RNG_kind); |
| for (j = 3; j < 6; j++) { |
| tmp = RNG_Table[RNG_kind].i_seed[j]; |
| if(tmp != 0) notallzero = 1; |
| if (tmp >= m2) allOK = 0; |
| } |
| if(!notallzero || !allOK) Randomize(RNG_kind); |
| } |
| break; |
| default: |
| error(_("FixupSeeds: unimplemented RNG kind %d"), RNG_kind); |
| } |
| } |
| |
| static void RNG_Init(RNGtype kind, Int32 seed) |
| { |
| int j; |
| |
| BM_norm_keep = 0.0; /* zap Box-Muller history */ |
| |
| /* Initial scrambling */ |
| for(j = 0; j < 50; j++) |
| seed = (69069 * seed + 1); |
| switch(kind) { |
| case WICHMANN_HILL: |
| case MARSAGLIA_MULTICARRY: |
| case SUPER_DUPER: |
| case MERSENNE_TWISTER: |
| /* i_seed[0] is mti, *but* this is needed for historical consistency */ |
| for(j = 0; j < RNG_Table[kind].n_seed; j++) { |
| seed = (69069 * seed + 1); |
| RNG_Table[kind].i_seed[j] = seed; |
| } |
| FixupSeeds(kind, 1); |
| break; |
| case KNUTH_TAOCP: |
| RNG_Init_R_KT(seed); |
| break; |
| case KNUTH_TAOCP2: |
| RNG_Init_KT2(seed); |
| break; |
| case LECUYER_CMRG: |
| for(j = 0; j < RNG_Table[kind].n_seed; j++) { |
| seed = (69069 * seed + 1); |
| while(seed >= m2) seed = (69069 * seed + 1); |
| RNG_Table[kind].i_seed[j] = seed; |
| } |
| break; |
| case USER_UNIF: |
| User_unif_fun = R_FindSymbol("user_unif_rand", "", NULL); |
| if (!User_unif_fun) error(_("'user_unif_rand' not in load table")); |
| User_unif_init = (UnifInitFun) R_FindSymbol("user_unif_init", "", NULL); |
| if (User_unif_init) (void) User_unif_init(seed); |
| User_unif_nseed = R_FindSymbol("user_unif_nseed", "", NULL); |
| User_unif_seedloc = R_FindSymbol("user_unif_seedloc", "", NULL); |
| if (User_unif_seedloc) { |
| int ns = 0; |
| if (!User_unif_nseed) { |
| warning(_("cannot read seeds unless 'user_unif_nseed' is supplied")); |
| break; |
| } |
| ns = *((int *) User_unif_nseed()); |
| if (ns < 0 || ns > 625) { |
| warning(_("seed length must be in 0...625; ignored")); |
| break; |
| } |
| RNG_Table[kind].n_seed = ns; |
| RNG_Table[kind].i_seed = (Int32 *) User_unif_seedloc(); |
| } |
| break; |
| default: |
| error(_("RNG_Init: unimplemented RNG kind %d"), kind); |
| } |
| } |
| |
| static SEXP GetSeedsFromVar(void) |
| { |
| SEXP seeds = findVarInFrame(R_GlobalEnv, R_SeedsSymbol); |
| if (TYPEOF(seeds) == PROMSXP) |
| seeds = eval(R_SeedsSymbol, R_GlobalEnv); |
| return seeds; |
| } |
| |
| static void Randomize(RNGtype kind) |
| { |
| /* Only called by GetRNGstate() when there is no .Random.seed */ |
| RNG_Init(kind, TimeToSeed()); |
| } |
| |
| static Rboolean GetRNGkind(SEXP seeds) |
| { |
| /* Load RNG_kind, N01_kind Sample_kind from .Random.seed if present */ |
| int tmp, *is; |
| RNGtype newRNG; N01type newN01; Sampletype newSample; |
| |
| if (isNull(seeds)) |
| seeds = GetSeedsFromVar(); |
| if (seeds == R_UnboundValue) return TRUE; |
| if (!isInteger(seeds)) { |
| if (seeds == R_MissingArg) /* How can this happen? */ |
| error(_("'.Random.seed' is a missing argument with no default")); |
| warning(_("'.Random.seed' is not an integer vector but of type '%s', so ignored"), |
| type2char(TYPEOF(seeds))); |
| goto invalid; |
| } |
| is = INTEGER(seeds); |
| tmp = is[0]; |
| /* avoid overflow here: max current value is 10705 */ |
| if (tmp == NA_INTEGER || tmp < 0 || tmp > 11000) { |
| warning(_("'.Random.seed[1]' is not a valid integer, so ignored")); |
| goto invalid; |
| } |
| newRNG = (RNGtype) (tmp % 100); |
| newN01 = (N01type) (tmp % 10000 / 100); |
| newSample = (Sampletype) (tmp / 10000); |
| if (newN01 > KINDERMAN_RAMAGE || newSample > REJECTION) { |
| warning(_("'.Random.seed[1]' is not a valid Normal type, so ignored")); |
| goto invalid; |
| } |
| switch(newRNG) { |
| case WICHMANN_HILL: |
| case MARSAGLIA_MULTICARRY: |
| case SUPER_DUPER: |
| case MERSENNE_TWISTER: |
| case KNUTH_TAOCP: |
| case KNUTH_TAOCP2: |
| case LECUYER_CMRG: |
| break; |
| case USER_UNIF: |
| if(!User_unif_fun) { |
| warning(_("'.Random.seed[1] = 5' but no user-supplied generator, so ignored")); |
| goto invalid; |
| } |
| break; |
| default: |
| warning(_("'.Random.seed[1]' is not a valid RNG kind so ignored")); |
| goto invalid; |
| } |
| RNG_kind = newRNG; N01_kind = newN01; Sample_kind = newSample; |
| return FALSE; |
| invalid: |
| RNG_kind = RNG_DEFAULT; N01_kind = N01_DEFAULT; Sample_kind = Sample_DEFAULT; |
| |
| Randomize(RNG_kind); |
| PutRNGstate(); // write out to .Random.seed |
| return TRUE; |
| } |
| |
| |
| void GetRNGstate() |
| { |
| /* Get .Random.seed into proper variables */ |
| int len_seed; |
| SEXP seeds; |
| |
| seeds = GetSeedsFromVar(); |
| if (seeds == R_UnboundValue) { |
| Randomize(RNG_kind); |
| } else { |
| /* this might re-set the generator */ |
| if(GetRNGkind(seeds)) return; |
| len_seed = RNG_Table[RNG_kind].n_seed; |
| /* Not sure whether this test is needed: wrong for USER_UNIF */ |
| if(LENGTH(seeds) > 1 && LENGTH(seeds) < len_seed + 1) |
| error(_("'.Random.seed' has wrong length")); |
| if(LENGTH(seeds) == 1 && RNG_kind != USER_UNIF) |
| Randomize(RNG_kind); |
| else { |
| int j, *is = INTEGER(seeds); |
| for(j = 1; j <= len_seed; j++) |
| RNG_Table[RNG_kind].i_seed[j - 1] = is[j]; |
| FixupSeeds(RNG_kind, 0); |
| } |
| } |
| } |
| |
| void PutRNGstate() |
| { |
| /* Copy out seeds to .Random.seed */ |
| int len_seed, j; |
| SEXP seeds; |
| |
| if (RNG_kind > LECUYER_CMRG || N01_kind > KINDERMAN_RAMAGE || Sample_kind > REJECTION) { |
| warning("Internal .Random.seed is corrupt: not saving"); |
| return; |
| } |
| |
| len_seed = RNG_Table[RNG_kind].n_seed; |
| |
| PROTECT(seeds = allocVector(INTSXP, len_seed + 1)); |
| |
| INTEGER(seeds)[0] = RNG_kind + 100 * N01_kind + 10000 * Sample_kind; |
| for(j = 0; j < len_seed; j++) |
| INTEGER(seeds)[j+1] = RNG_Table[RNG_kind].i_seed[j]; |
| |
| /* assign only in the workspace */ |
| defineVar(R_SeedsSymbol, seeds, R_GlobalEnv); |
| UNPROTECT(1); |
| } |
| |
| static void RNGkind(RNGtype newkind) |
| { |
| /* Choose a new kind of RNG. |
| * Initialize its seed by calling the old RNG's unif_rand() |
| */ |
| if (newkind == (RNGtype)-1) newkind = RNG_DEFAULT; |
| switch(newkind) { |
| case WICHMANN_HILL: |
| case MARSAGLIA_MULTICARRY: |
| case SUPER_DUPER: |
| case MERSENNE_TWISTER: |
| case KNUTH_TAOCP: |
| case USER_UNIF: |
| case KNUTH_TAOCP2: |
| case LECUYER_CMRG: |
| break; |
| default: |
| error(_("RNGkind: unimplemented RNG kind %d"), newkind); |
| } |
| GetRNGstate(); |
| // precaution against corruption as per package randtoolbox |
| double u = unif_rand(); |
| if (u < 0.0 || u > 1.0) { |
| warning("someone corrupted the random-number generator: re-initializing"); |
| RNG_Init(newkind, TimeToSeed()); |
| } else |
| RNG_Init(newkind, (Int32) (u * UINT_MAX)); |
| RNG_kind = newkind; |
| PutRNGstate(); |
| } |
| |
| static void Norm_kind(N01type kind) |
| { |
| /* N01type is an enumeration type, so this will probably get |
| mapped to an unsigned integer type. */ |
| if (kind == (N01type)-1) kind = N01_DEFAULT; |
| if (kind > KINDERMAN_RAMAGE) |
| error(_("invalid Normal type in 'RNGkind'")); |
| if (kind == USER_NORM) { |
| User_norm_fun = R_FindSymbol("user_norm_rand", "", NULL); |
| if (!User_norm_fun) error(_("'user_norm_rand' not in load table")); |
| } |
| GetRNGstate(); /* might not be initialized */ |
| if (kind == BOX_MULLER) |
| BM_norm_keep = 0.0; /* zap Box-Muller history */ |
| N01_kind = kind; |
| PutRNGstate(); |
| } |
| |
| static void Samp_kind(Sampletype kind) |
| { |
| /* Sampletype is an enumeration type, so this will probably get |
| mapped to an unsigned integer type. */ |
| if (kind == (Sampletype)-1) kind = Sample_DEFAULT; |
| if (kind > REJECTION) |
| error(_("invalid sample type in 'RNGkind'")); |
| GetRNGstate(); /* might not be initialized */ |
| Sample_kind = kind; |
| PutRNGstate(); |
| } |
| |
| |
| /*------ .Internal interface ------------------------*/ |
| |
| SEXP attribute_hidden do_RNGkind (SEXP call, SEXP op, SEXP args, SEXP env) |
| { |
| SEXP ans, rng, norm, sample; |
| |
| checkArity(op,args); |
| GetRNGstate(); /* might not be initialized */ |
| PROTECT(ans = allocVector(INTSXP, 3)); |
| INTEGER(ans)[0] = RNG_kind; |
| INTEGER(ans)[1] = N01_kind; |
| INTEGER(ans)[2] = Sample_kind; |
| rng = CAR(args); |
| norm = CADR(args); |
| sample = CADDR(args); |
| GetRNGkind(R_NilValue); /* pull from .Random.seed if present */ |
| if(!isNull(rng)) { /* set a new RNG kind */ |
| RNGkind((RNGtype) asInteger(rng)); |
| } |
| if(!isNull(norm)) { /* set a new normal kind */ |
| Norm_kind((N01type) asInteger(norm)); |
| } |
| if(!isNull(sample)) { /* set a new sample kind */ |
| Samp_kind((Sampletype) asInteger(sample)); |
| } |
| UNPROTECT(1); |
| return ans; |
| } |
| |
| |
| SEXP attribute_hidden do_setseed (SEXP call, SEXP op, SEXP args, SEXP env) |
| { |
| SEXP skind, nkind, sampkind; |
| int seed; |
| |
| checkArity(op, args); |
| if(!isNull(CAR(args))) { |
| seed = asInteger(CAR(args)); |
| if (seed == NA_INTEGER) |
| error(_("supplied seed is not a valid integer")); |
| } else seed = TimeToSeed(); |
| skind = CADR(args); |
| nkind = CADDR(args); |
| sampkind = CADDDR(args); |
| GetRNGkind(R_NilValue); /* pull RNG_kind, N01_kind from |
| .Random.seed if present */ |
| if (!isNull(skind)) RNGkind((RNGtype) asInteger(skind)); |
| if (!isNull(nkind)) Norm_kind((N01type) asInteger(nkind)); |
| if(!isNull(sampkind)) Samp_kind((Sampletype) asInteger(sampkind)); |
| RNG_Init(RNG_kind, (Int32) seed); /* zaps BM history */ |
| PutRNGstate(); |
| return R_NilValue; |
| } |
| |
| |
| /* S COMPATIBILITY */ |
| |
| /* The following entry points provide compatibility with S. */ |
| /* These entry points should not be used by new R code. */ |
| |
| void seed_in(long *ignored) |
| { |
| GetRNGstate(); |
| } |
| |
| void seed_out(long *ignored) |
| { |
| PutRNGstate(); |
| } |
| |
| /* =================== Mersenne Twister ========================== */ |
| /* From http://www.math.keio.ac.jp/~matumoto/emt.html */ |
| /* New URL (accessed 2018-11-08): |
| http://www.math.sci.hiroshima-u.ac.jp/~m-mat/eindex.html |
| |
| The initialization method in the 1998 code and paper had a minor |
| issue that was addressed with new initialization approaches in an |
| update in 2002. R has always used a different initialization |
| approach and is not affected by that issue. |
| */ |
| |
| /* A C-program for MT19937: Real number version([0,1)-interval) |
| (1999/10/28) |
| genrand() generates one pseudorandom real number (double) |
| which is uniformly distributed on [0,1)-interval, for each |
| call. sgenrand(seed) sets initial values to the working area |
| of 624 words. Before genrand(), sgenrand(seed) must be |
| called once. (seed is any 32-bit integer.) |
| Integer generator is obtained by modifying two lines. |
| Coded by Takuji Nishimura, considering the suggestions by |
| Topher Cooper and Marc Rieffel in July-Aug. 1997. |
| |
| Copyright (C) 1997, 1999 Makoto Matsumoto and Takuji Nishimura. |
| When you use this, send an email to: matumoto@math.keio.ac.jp |
| with an appropriate reference to your work. |
| |
| REFERENCE |
| M. Matsumoto and T. Nishimura, |
| "Mersenne Twister: A 623-Dimensionally Equidistributed Uniform |
| Pseudo-Random Number Generator", |
| ACM Transactions on Modeling and Computer Simulation, |
| Vol. 8, No. 1, January 1998, pp 3--30. |
| */ |
| |
| /* Period parameters */ |
| #define N 624 |
| #define M 397 |
| #define MATRIX_A 0x9908b0df /* constant vector a */ |
| #define UPPER_MASK 0x80000000 /* most significant w-r bits */ |
| #define LOWER_MASK 0x7fffffff /* least significant r bits */ |
| |
| /* Tempering parameters */ |
| #define TEMPERING_MASK_B 0x9d2c5680 |
| #define TEMPERING_MASK_C 0xefc60000 |
| #define TEMPERING_SHIFT_U(y) (y >> 11) |
| #define TEMPERING_SHIFT_S(y) (y << 7) |
| #define TEMPERING_SHIFT_T(y) (y << 15) |
| #define TEMPERING_SHIFT_L(y) (y >> 18) |
| |
| static Int32 *mt = dummy+1; /* the array for the state vector */ |
| static int mti=N+1; /* mti==N+1 means mt[N] is not initialized */ |
| |
| /* Initializing the array with a seed */ |
| static void |
| MT_sgenrand(Int32 seed) |
| { |
| int i; |
| |
| for (i = 0; i < N; i++) { |
| mt[i] = seed & 0xffff0000; |
| seed = 69069 * seed + 1; |
| mt[i] |= (seed & 0xffff0000) >> 16; |
| seed = 69069 * seed + 1; |
| } |
| mti = N; |
| } |
| |
| /* Initialization by "sgenrand()" is an example. Theoretically, |
| there are 2^19937-1 possible states as an intial state. |
| Essential bits in "seed_array[]" is following 19937 bits: |
| (seed_array[0]&UPPER_MASK), seed_array[1], ..., seed_array[N-1]. |
| (seed_array[0]&LOWER_MASK) is discarded. |
| Theoretically, |
| (seed_array[0]&UPPER_MASK), seed_array[1], ..., seed_array[N-1] |
| can take any values except all zeros. */ |
| |
| static double MT_genrand(void) |
| { |
| Int32 y; |
| static Int32 mag01[2]={0x0, MATRIX_A}; |
| /* mag01[x] = x * MATRIX_A for x=0,1 */ |
| |
| mti = dummy[0]; |
| |
| if (mti >= N) { /* generate N words at one time */ |
| int kk; |
| |
| if (mti == N+1) /* if sgenrand() has not been called, */ |
| MT_sgenrand(4357); /* a default initial seed is used */ |
| |
| for (kk = 0; kk < N - M; kk++) { |
| y = (mt[kk] & UPPER_MASK) | (mt[kk+1] & LOWER_MASK); |
| mt[kk] = mt[kk+M] ^ (y >> 1) ^ mag01[y & 0x1]; |
| } |
| for (; kk < N - 1; kk++) { |
| y = (mt[kk] & UPPER_MASK) | (mt[kk+1] & LOWER_MASK); |
| mt[kk] = mt[kk+(M-N)] ^ (y >> 1) ^ mag01[y & 0x1]; |
| } |
| y = (mt[N-1] & UPPER_MASK) | (mt[0] & LOWER_MASK); |
| mt[N-1] = mt[M-1] ^ (y >> 1) ^ mag01[y & 0x1]; |
| |
| mti = 0; |
| } |
| |
| y = mt[mti++]; |
| y ^= TEMPERING_SHIFT_U(y); |
| y ^= TEMPERING_SHIFT_S(y) & TEMPERING_MASK_B; |
| y ^= TEMPERING_SHIFT_T(y) & TEMPERING_MASK_C; |
| y ^= TEMPERING_SHIFT_L(y); |
| dummy[0] = mti; |
| |
| return ( (double)y * 2.3283064365386963e-10 ); /* reals: [0,1)-interval */ |
| } |
| |
| /* |
| The following code was taken from earlier versions of |
| http://www-cs-faculty.stanford.edu/~knuth/programs/rng.c-old |
| http://www-cs-faculty.stanford.edu/~knuth/programs/rng.c |
| */ |
| |
| |
| /* This define may give a warning with clang, but is needed to comply |
| with the prohibition on changing the code. */ |
| #ifdef __clang__ |
| #pragma clang diagnostic push |
| #pragma clang diagnostic ignored "-Wkeyword-macro" |
| #endif |
| #define long Int32 |
| #ifdef __clang__ |
| #pragma clang diagnostic pop |
| #endif |
| |
| #define ran_arr_buf R_KT_ran_arr_buf |
| #define ran_arr_cycle R_KT_ran_arr_cycle |
| #define ran_arr_ptr R_KT_ran_arr_ptr |
| #define ran_arr_sentinel R_KT_ran_arr_sentinel |
| #define ran_x dummy |
| |
| #define KK 100 /* the long lag */ |
| #define LL 37 /* the short lag */ |
| #define MM (1L<<30) /* the modulus */ |
| #define TT 70 /* guaranteed separation between streams */ |
| #define mod_diff(x,y) (((x)-(y))&(MM-1)) /* subtraction mod MM */ |
| #define is_odd(x) ((x)&1) /* units bit of x */ |
| static void ran_array(long aa[],int n) /* put n new random numbers in aa */ |
| { |
| register int i,j; |
| for (j=0;j<KK;j++) aa[j]=ran_x[j]; |
| for (;j<n;j++) aa[j]=mod_diff(aa[j-KK],aa[j-LL]); |
| for (i=0;i<LL;i++,j++) ran_x[i]=mod_diff(aa[j-KK],aa[j-LL]); |
| for (;i<KK;i++,j++) ran_x[i]=mod_diff(aa[j-KK],ran_x[i-LL]); |
| } |
| #define QUALITY 1009 /* recommended quality level for high-res use */ |
| static long ran_arr_buf[QUALITY]; |
| static long ran_arr_sentinel=(long)-1; |
| static long *ran_arr_ptr=&ran_arr_sentinel; /* the next random number, or -1 */ |
| |
| static long ran_arr_cycle(void) |
| { |
| ran_array(ran_arr_buf,QUALITY); |
| ran_arr_buf[KK]=(long)(-1); |
| ran_arr_ptr=ran_arr_buf+1; |
| return ran_arr_buf[0]; |
| } |
| |
| /* =================== Knuth TAOCP 2002 ========================== */ |
| |
| /* This program by D E Knuth is in the public domain and freely copyable. |
| * It is explained in Seminumerical Algorithms, 3rd edition, Section 3.6 |
| * (or in the errata to the 2nd edition --- see |
| * http://www-cs-faculty.stanford.edu/~knuth/taocp.html |
| * in the changes to Volume 2 on pages 171 and following). */ |
| |
| /* N.B. The MODIFICATIONS introduced in the 9th printing (2002) are |
| included here; there's no backwards compatibility with the original. */ |
| |
| |
| static void ran_start(long seed) |
| { |
| register int t,j; |
| long x[KK+KK-1]; /* the preparation buffer */ |
| register long ss=(seed+2)&(MM-2); |
| for (j=0;j<KK;j++) { |
| x[j]=ss; /* bootstrap the buffer */ |
| ss<<=1; if (ss>=MM) ss-=MM-2; /* cyclic shift 29 bits */ |
| } |
| x[1]++; /* make x[1] (and only x[1]) odd */ |
| for (ss=seed&(MM-1),t=TT-1; t; ) { |
| for (j=KK-1;j>0;j--) x[j+j]=x[j], x[j+j-1]=0; /* "square" */ |
| for (j=KK+KK-2;j>=KK;j--) |
| x[j-(KK-LL)]=mod_diff(x[j-(KK-LL)],x[j]), |
| x[j-KK]=mod_diff(x[j-KK],x[j]); |
| if (is_odd(ss)) { /* "multiply by z" */ |
| for (j=KK;j>0;j--) x[j]=x[j-1]; |
| x[0]=x[KK]; /* shift the buffer cyclically */ |
| x[LL]=mod_diff(x[LL],x[KK]); |
| } |
| if (ss) ss>>=1; else t--; |
| } |
| for (j=0;j<LL;j++) ran_x[j+KK-LL]=x[j]; |
| for (;j<KK;j++) ran_x[j-LL]=x[j]; |
| for (j=0;j<10;j++) ran_array(x,KK+KK-1); /* warm things up */ |
| ran_arr_ptr=&ran_arr_sentinel; |
| } |
| /* ===================== end of Knuth's code ====================== */ |
| |
| static void RNG_Init_KT2(Int32 seed) |
| { |
| ran_start(seed % 1073741821); |
| KT_pos = 100; |
| } |
| |
| static Int32 KT_next(void) |
| { |
| if(KT_pos >= 100) { |
| ran_arr_cycle(); |
| KT_pos = 0; |
| } |
| return ran_x[(KT_pos)++]; |
| } |
| |
| static void RNG_Init_R_KT(Int32 seed) |
| { |
| SEXP fun, sseed, call, ans; |
| PROTECT(fun = findVar1(install(".TAOCP1997init"), R_BaseEnv, CLOSXP, FALSE)); |
| if(fun == R_UnboundValue) |
| error("function '.TAOCP1997init' is missing"); |
| PROTECT(sseed = ScalarInteger((int)(seed % 1073741821))); |
| PROTECT(call = lang2(fun, sseed)); |
| ans = eval(call, R_GlobalEnv); |
| memcpy(dummy, INTEGER(ans), 100*sizeof(int)); |
| UNPROTECT(3); |
| KT_pos = 100; |
| } |
| |
| /* Our PRNGs have at most 32 bit of precision. All generators except |
| Knuth-TAOCP, Knuth-TAOCP-2002, and possibly the user-supplied ones |
| have 31 or 32 bits bits of precision; the others are assumed to |
| have at least at least 25. */ |
| static R_INLINE double ru() |
| { |
| double U = 33554432.0; |
| return (floor(U*unif_rand()) + unif_rand())/U; |
| } |
| |
| static double R_unif_index_0(double dn) |
| { |
| double cut = INT_MAX; |
| |
| switch(RNG_kind) { |
| case KNUTH_TAOCP: |
| case USER_UNIF: |
| case KNUTH_TAOCP2: |
| cut = 33554431.0; /* 2^25 - 1 */ |
| break; |
| default: |
| break; |
| } |
| |
| double u = dn > cut ? ru() : unif_rand(); |
| return floor(dn * u); |
| } |
| |
| //generate a random non-negative integer < 2 ^ bits in 16 bit chunks |
| static double rbits(int bits) |
| { |
| int_least64_t v = 0; |
| for (int n = 0; n <= bits; n += 16) { |
| int v1 = (int) floor(unif_rand() * 65536); |
| v = 65536 * v + v1; |
| } |
| const int_least64_t one64 = 1L; |
| // mask out the bits in the result that are not needed |
| return (double) (v & ((one64 << bits) - 1)); |
| } |
| |
| double R_unif_index(double dn) |
| { |
| if (Sample_kind == ROUNDING) |
| return R_unif_index_0(dn); |
| |
| // rejection sampling from integers below the next larger power of two |
| if (dn <= 0) |
| return 0.0; |
| int bits = (int) ceil(log2(dn)); |
| double dv; |
| do { dv = rbits(bits); } while (dn <= dv); |
| return dv; |
| } |
| |
| Sampletype R_sample_kind() { return Sample_kind; } |
| |