text: /*
  	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);
  }
  
  