Imported Upstream version 3.0+20111124upstream/3.0+20111124upstream

1 files changed, 345 insertions, 0 deletions

diff --git a/src/common/Random.c b/src/common/Random.c
new file mode 100644
index 0000000..66b180e
--- /dev/null
+++ b/src/common/Random.c
@@ -0,0 +1,345 @@
+/*
+	Random.c
+		quasi- and pseudo-random-number generation
+		last modified 8 Jun 10 th
+*/
+
+
+/*
+	PART 1: Sobol quasi-random-number generator
+	adapted from ACM TOMS algorithm 659
+*/
+
+static void SobolGet(This *t, real *x)
+{
+  number seq = t->rng.sobol.seq++;
+  count zerobit = 0, dim;
+
+  while( seq & 1 ) {
+    ++zerobit;
+    seq >>= 1;
+  }
+
+  for( dim = 0; dim < t->ndim; ++dim ) {
+    t->rng.sobol.prev[dim] ^= t->rng.sobol.v[dim][zerobit];
+    x[dim] = t->rng.sobol.prev[dim]*t->rng.sobol.norm;
+  }
+}
+
+
+static void SobolSkip(This *t, number n)
+{
+  while( n-- ) {
+    number seq = t->rng.sobol.seq++;
+    count zerobit = 0, dim;
+
+    while( seq & 1 ) {
+      ++zerobit;
+      seq >>= 1;
+    }
+
+    for( dim = 0; dim < t->ndim; ++dim )
+      t->rng.sobol.prev[dim] ^= t->rng.sobol.v[dim][zerobit];
+  }
+}
+
+
+static inline void SobolIni(This *t)
+{
+  static number ini[9*40] = {
+      3,   1,   0,   0,   0,   0,   0,   0,   0,
+      7,   1,   1,   0,   0,   0,   0,   0,   0,
+     11,   1,   3,   7,   0,   0,   0,   0,   0,
+     13,   1,   1,   5,   0,   0,   0,   0,   0,
+     19,   1,   3,   1,   1,   0,   0,   0,   0,
+     25,   1,   1,   3,   7,   0,   0,   0,   0,
+     37,   1,   3,   3,   9,   9,   0,   0,   0,
+     59,   1,   3,   7,  13,   3,   0,   0,   0,
+     47,   1,   1,   5,  11,  27,   0,   0,   0,
+     61,   1,   3,   5,   1,  15,   0,   0,   0,
+     55,   1,   1,   7,   3,  29,   0,   0,   0,
+     41,   1,   3,   7,   7,  21,   0,   0,   0,
+     67,   1,   1,   1,   9,  23,  37,   0,   0,
+     97,   1,   3,   3,   5,  19,  33,   0,   0,
+     91,   1,   1,   3,  13,  11,   7,   0,   0,
+    109,   1,   1,   7,  13,  25,   5,   0,   0,
+    103,   1,   3,   5,  11,   7,  11,   0,   0,
+    115,   1,   1,   1,   3,  13,  39,   0,   0,
+    131,   1,   3,   1,  15,  17,  63,  13,   0,
+    193,   1,   1,   5,   5,   1,  27,  33,   0,
+    137,   1,   3,   3,   3,  25,  17, 115,   0,
+    145,   1,   1,   3,  15,  29,  15,  41,   0,
+    143,   1,   3,   1,   7,   3,  23,  79,   0,
+    241,   1,   3,   7,   9,  31,  29,  17,   0,
+    157,   1,   1,   5,  13,  11,   3,  29,   0,
+    185,   1,   3,   1,   9,   5,  21, 119,   0,
+    167,   1,   1,   3,   1,  23,  13,  75,   0,
+    229,   1,   3,   3,  11,  27,  31,  73,   0,
+    171,   1,   1,   7,   7,  19,  25, 105,   0,
+    213,   1,   3,   5,   5,  21,   9,   7,   0,
+    191,   1,   1,   1,  15,   5,  49,  59,   0,
+    253,   1,   1,   1,   1,   1,  33,  65,   0,
+    203,   1,   3,   5,  15,  17,  19,  21,   0,
+    211,   1,   1,   7,  11,  13,  29,   3,   0,
+    239,   1,   3,   7,   5,   7,  11, 113,   0,
+    247,   1,   1,   5,   3,  15,  19,  61,   0,
+    285,   1,   3,   1,   1,   9,  27,  89,   7,
+    369,   1,   1,   3,   7,  31,  15,  45,  23,
+    299,   1,   3,   3,   9,   9,  25, 107,  39 };
+
+  count dim, bit, nbits;
+  number max, *pini = ini;
+  cnumber nmax = 2*t->maxeval;
+
+  for( nbits = 0, max = 1; max <= nmax; max <<= 1 ) ++nbits;
+  t->rng.sobol.norm = 1./max;
+
+  for( bit = 0; bit < nbits; ++bit )
+    t->rng.sobol.v[0][bit] = (max >>= 1);
+
+  for( dim = 1; dim < t->ndim; ++dim ) {
+    number *pv = t->rng.sobol.v[dim], *pvv = pv;
+    number powers = *pini++, j;
+    int inibits = -1, bit;
+    for( j = powers; j; j >>= 1 ) ++inibits;
+
+    memcpy(pv, pini, inibits*sizeof(*pini));
+    pini += 8;
+
+    for( bit = inibits; bit < nbits; ++bit ) {
+      number newv = *pvv, j = powers;
+      int b;
+      for( b = 0; b < inibits; ++b ) {
+        if( j & 1 ) newv ^= pvv[b] << (inibits - b);
+        j >>= 1;
+      }
+      pvv[inibits] = newv;
+      ++pvv;
+    }
+
+    for( bit = 0; bit < nbits - 1; ++bit )
+      pv[bit] <<= nbits - bit - 1;
+  }
+
+  t->rng.sobol.seq = 0;
+  VecClear(t->rng.sobol.prev);
+
+  t->rng.getrandom = SobolGet;
+  t->rng.skiprandom = SobolSkip;
+}
+
+
+/*
+	PART 2: Mersenne Twister pseudo-random-number generator
+	adapted from T. Nishimura's and M. Matsumoto's C code at
+	http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/emt.html
+*/
+
+/* 32 or 53 random bits */
+#define RANDOM_BITS 32
+
+
+static inline state_t Twist(state_t a, state_t b)
+{
+  state_t mixbits = (a & 0x80000000) | (b & 0x7fffffff);
+  state_t matrixA = (-(b & 1)) & 0x9908b0df;
+  return (mixbits >> 1) ^ matrixA;
+}
+
+
+static inline void MersenneReload(state_t *state)
+{
+  state_t *s = state;
+  int j;
+
+  for( j = MERSENNE_N - MERSENNE_M + 1; --j; ++s )
+    *s = s[MERSENNE_M] ^ Twist(s[0], s[1]);
+  for( j = MERSENNE_M; --j; ++s )
+    *s = s[MERSENNE_M - MERSENNE_N] ^ Twist(s[0], s[1]);
+  *s = s[MERSENNE_M - MERSENNE_N] ^ Twist(s[0], state[0]);
+}
+
+
+static inline state_t MersenneInt(state_t s)
+{
+  s ^= s >> 11;
+  s ^= (s << 7) & 0x9d2c5680;
+  s ^= (s << 15) & 0xefc60000;
+  return s ^ (s >> 18);
+}
+
+
+static void MersenneGet(This *t, real *x)
+{
+  count next = t->rng.mersenne.next, dim;
+
+  for( dim = 0; dim < t->ndim; ++dim ) {
+#if RANDOM_BITS == 53
+    state_t a, b;
+#endif
+
+    if( next >= MERSENNE_N ) {
+      MersenneReload(t->rng.mersenne.state);
+      next = 0;
+    }
+
+#if RANDOM_BITS == 53
+    a = MersenneInt(t->rng.mersenne.state[next++]) >> 5;
+    b = MersenneInt(t->rng.mersenne.state[next++]) >> 6;
+    x[dim] = (67108864.*a + b)/9007199254740992.;
+#else
+    x[dim] = MersenneInt(t->rng.mersenne.state[next++])/4294967296.;
+#endif
+  }
+
+  t->rng.mersenne.next = next;
+}
+
+
+static void MersenneSkip(This *t, number n)
+{
+#if RANDOM_BITS == 53
+  n = 2*n*t->ndim + t->rng.mersenne.next;
+#else
+  n = n*t->ndim + t->rng.mersenne.next;
+#endif
+  t->rng.mersenne.next = n % MERSENNE_N;
+  n /= MERSENNE_N;
+  while( n-- ) MersenneReload(t->rng.mersenne.state);
+}
+
+
+static inline void MersenneIni(This *t)
+{
+  state_t seed = t->seed;
+  state_t *next = t->rng.mersenne.state;
+  count j;
+
+  for( j = 1; j <= MERSENNE_N; ++j ) {
+    *next++ = seed;
+    seed = 0x6c078965*(seed ^ (seed >> 30)) + j;
+    /* see Knuth TAOCP Vol 2, 3rd Ed, p. 106 for multiplier */
+  }
+
+  MersenneReload(t->rng.mersenne.state);
+  t->rng.mersenne.next = 0;
+
+  t->rng.getrandom = MersenneGet;
+  t->rng.skiprandom = MersenneSkip;
+}
+
+
+/*
+	PART 3: Ranlux subtract-and-borrow random-number generator 
+	proposed by Marsaglia and Zaman, implemented by F. James with 
+	the name RCARRY in 1991, and later improved by Martin Luescher 
+	in 1993 to produce "Luxury Pseudorandom Numbers".
+	Adapted from the CERNlib Fortran 77 code by F. James, 1993.
+
+	The available luxury levels are:
+
+	level 0  (p = 24): equivalent to the original RCARRY of Marsaglia
+	         and Zaman, very long period, but fails many tests.
+	level 1  (p = 48): considerable improvement in quality over level 0,
+	         now passes the gap test, but still fails spectral test.
+	level 2  (p = 97): passes all known tests, but theoretically still
+	         defective.
+	level 3  (p = 223): DEFAULT VALUE.  Any theoretically possible
+	         correlations have very small chance of being observed.
+	level 4  (p = 389): highest possible luxury, all 24 bits chaotic.
+*/
+
+
+static inline int RanluxInt(This *t, count n)
+{
+  int s = 0;
+
+  while( n-- ) {
+    s = t->rng.ranlux.state[t->rng.ranlux.j24] -
+        t->rng.ranlux.state[t->rng.ranlux.i24] + t->rng.ranlux.carry;
+    s += (t->rng.ranlux.carry = NegQ(s)) & (1 << 24);
+    t->rng.ranlux.state[t->rng.ranlux.i24] = s;
+    --t->rng.ranlux.i24;
+    t->rng.ranlux.i24 += NegQ(t->rng.ranlux.i24) & 24;
+    --t->rng.ranlux.j24;
+    t->rng.ranlux.j24 += NegQ(t->rng.ranlux.j24) & 24;
+  }
+
+  return s;
+}
+
+
+static void RanluxGet(This *t, real *x)
+{
+/* The Generator proper: "Subtract-with-borrow",
+   as proposed by Marsaglia and Zaman, FSU, March 1989 */
+
+  count dim;
+
+  for( dim = 0; dim < t->ndim; ++dim ) {
+    cint nskip = (--t->rng.ranlux.n24 >= 0) ? 0 :
+      (t->rng.ranlux.n24 = 24, t->rng.ranlux.nskip);
+    cint s = RanluxInt(t, 1 + nskip);
+    x[dim] = s*0x1p-24;
+/* small numbers (with less than 12 significant bits) are "padded" */
+    if( s < (1 << 12) )
+      x[dim] += t->rng.ranlux.state[t->rng.ranlux.j24]*0x1p-48;
+  }
+}
+
+
+static void RanluxSkip(This *t, cnumber n)
+{
+  RanluxInt(t, n + t->rng.ranlux.nskip*(n/24));
+  t->rng.ranlux.n24 = 24 - n % 24;
+}
+
+
+static inline void RanluxIni(This *t)
+{
+  cint skip[] = {24, 48, 97, 223, 389,
+    223, 223, 223, 223, 223, 223, 223, 223, 223, 223,
+    223, 223, 223, 223, 223, 223, 223, 223, 223, 223};
+  state_t seed = t->seed;
+  state_t level = RNG;
+  count i;
+
+  if( level < sizeof skip ) level = skip[level];
+  t->rng.ranlux.nskip = level - 24;
+
+  t->rng.ranlux.i24 = 23;
+  t->rng.ranlux.j24 = 9;
+  t->rng.ranlux.n24 = 24;
+
+  for( i = 0; i < 24; ++i ) {
+    cint k = seed/53668;
+    seed = 40014*(seed - k*53668) - k*12211;
+    seed += NegQ(seed) & 2147483563;
+    t->rng.ranlux.state[i] = seed & ((1 << 24) - 1);
+  }
+
+  t->rng.ranlux.carry = ~TrueQ(t->rng.ranlux.state[23]) & (1 << 24);
+
+  t->rng.getrandom = RanluxGet;
+  t->rng.skiprandom = RanluxSkip;
+}
+
+
+/*
+	PART 4: User routines:
+
+	- IniRandom sets up the random-number generator to produce a
+	  sequence of at least n ndim-dimensional random vectors.
+
+	- GetRandom retrieves one random vector.
+
+	- SkipRandom skips over n random vectors.
+*/
+
+static inline void IniRandom(This *t)
+{
+  if( t->seed == 0 ) SobolIni(t);
+  else if( RNG == 0 ) MersenneIni(t);
+  else RanluxIni(t);
+}
+