text: /*
  	Random.c
  		quasi- and pseudo-random-number generation
  		last modified 5 Jan 06 th
  */
  
  
  /*
  	PART 1: Sobol quasi-random-number generator
  	adapted from ACM TOMS algorithm 659
  */
  
  #define SOBOL_MINDIM 1
  #define SOBOL_MAXDIM 40
  
  static struct {
    real norm;
    number v[SOBOL_MAXDIM][30], prev[SOBOL_MAXDIM];
    number seq;
  } sobol_;
  
  
  static inline void SobolIni(cnumber n)
  {
    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;
  
    for( nbits = 0, max = 1; max <= n; max <<= 1 ) ++nbits;
    sobol_.norm = 1./max;
  
    for( bit = 0; bit < nbits; ++bit )
      sobol_.v[0][bit] = (max >>= 1);
  
    for( dim = 1; dim < ndim_; ++dim ) {
      number *pv = 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;
    }
  
    sobol_.seq = 0;
    VecClear(sobol_.prev);
  }
  
  
  static inline void SobolGet(real *x)
  {
    number seq = sobol_.seq++;
    count zerobit = 0, dim;
  
    while( seq & 1 ) {
      ++zerobit;
      seq >>= 1;
    }
  
    for( dim = 0; dim < ndim_; ++dim ) {
      sobol_.prev[dim] ^= sobol_.v[dim][zerobit];
      x[dim] = sobol_.prev[dim]*sobol_.norm;
    }
  }
  
  
  static inline void SobolSkip(number n)
  {
    while( n-- ) {
      number seq = sobol_.seq++;
      count zerobit = 0, dim;
  
      while( seq & 1 ) {
        ++zerobit;
        seq >>= 1;
      }
  
      for( dim = 0; dim < ndim_; ++dim )
        sobol_.prev[dim] ^= sobol_.v[dim][zerobit];
    }
  }
  
  
  /*
  	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
  */
  
  /* length of state vector */
  #define MERSENNE_N 624
  
  /* period parameter */
  #define MERSENNE_M 397
  
  #define DEFAULT_SEED 5489
  
  /* 32 or 53 random bits */
  #define RANDOM_BITS 32
  
  typedef unsigned int state_t;
  
  static struct {
    state_t state[MERSENNE_N];
    count next;
  } mersenne_;
  
  
  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 *s = mersenne_.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], mersenne_.state[0]);
  }
  
  
  static inline void MersenneIni(state_t seed)
  {
    state_t *next = mersenne_.state;
    int 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();
    mersenne_.next = 0;
  }
  
  
  static inline state_t MersenneInt(count next)
  {
    state_t s = mersenne_.state[next];
    s ^= s >> 11;
    s ^= (s << 7) & 0x9d2c5680;
    s ^= (s << 15) & 0xefc60000;
    return s ^ (s >> 18);
  }
  
  
  static inline void MersenneGet(real *x)
  {
    count next = mersenne_.next, dim;
  
    for( dim = 0; dim < ndim_; ++dim ) {
  #if RANDOM_BITS == 53
      state_t a, b;
  #endif
  
      if( next >= MERSENNE_N ) {
        MersenneReload();
        next = 0;
      }
  
  #if RANDOM_BITS == 53
      a = MersenneInt(next++) >> 5;
      b = MersenneInt(next++) >> 6;
      x[dim] = (67108864.*a + b)/9007199254740992.;
  #else
      x[dim] = MersenneInt(next++)/4294967295.;
  #endif
    }
  
    mersenne_.next = next;
  }
  
  
  static inline void MersenneSkip(number n)
  {
  #if RANDOM_BITS == 53
    n = 2*n*ndim_ + mersenne_.next;
  #else
    n = n*ndim_ + mersenne_.next;
  #endif
    mersenne_.next = n % MERSENNE_N;
    n /= MERSENNE_N;
    while( n-- ) MersenneReload();
  }
  
  
  /*
  	PART 3: 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 void IniRandom(cnumber n, cint flags)
  {
    if( PSEUDORNG ) {
      sobol_.seq = -1;
      MersenneIni(DEFAULT_SEED);
    }
    else SobolIni(n);
  }
  
  static inline void GetRandom(real *x)
  {
    if( sobol_.seq == -1 ) MersenneGet(x);
    else SobolGet(x);
  }
  
  static inline void SkipRandom(cnumber n)
  {
    if( sobol_.seq == -1 ) MersenneSkip(n);
    else SobolSkip(n);
  }
  
  