Source

enzo-3.0 / src / enzo / utilities / MersenneTwister.C

The active_particles branch has multiple heads

/* 
   This is an implementation of the Mersenne Twister which was written by  Michael
   Brundage. The original source code, which was copied and pasted into this file,
   can be found at http://qbrundage.com/michaelb/pubs/essays/random_number_generation,
   and a description of the Mersenne Twister random number generator can be found here:
   http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/emt.html
   It has a claimed periodicity of 2^{19937}-1, which is pretty incredible.  The algorithm
   as-is is mainly intended for Monte Carlo realizations, and is not intended to be used for
   cryptography.  See the website for more information.

   This code was placed into the public domain by Michael Brundage, and was put into Enzo
   by Brian W. O'Shea on 11 December 2007.  It uses the system random number generator as
   an initial seed, and the user must specify a seed for the system random number generator.

*/
#include "preincludes.h"
#include "macros_and_parameters.h"

#define MT_LEN       624

int mt_index;
unsigned_long_int mt_buffer[MT_LEN];

void mt_init(unsigned_int seed) {
  int i;

  srand(seed);

  for (i = 0; i < MT_LEN; i++)
    mt_buffer[i] = ((unsigned_long_int) rand() );
  mt_index = 0;
}

#define MT_IA           397
#define MT_IB           (MT_LEN - MT_IA)
#define UPPER_MASK      0x80000000
#define LOWER_MASK      0x7FFFFFFF
#define MATRIX_A        0x9908B0DF
#define TWIST(b,i,j)    ((b)[i] & UPPER_MASK) | ((b)[j] & LOWER_MASK)
#define MAGIC(s)        (((s)&1)*MATRIX_A)

unsigned_long_int mt_random() {
  unsigned_long_int * b = mt_buffer;
  int idx = mt_index;
  unsigned_long_int s;
  int i;

  if (idx == MT_LEN*sizeof(unsigned_long_int))
    {
      idx = 0;
      i = 0;
      for (; i < MT_IB; i++) {
	s = TWIST(b, i, i+1);
	b[i] = b[i + MT_IA] ^ (s >> 1) ^ MAGIC(s);
      }
      for (; i < MT_LEN-1; i++) {
	s = TWIST(b, i, i+1);
	b[i] = b[i - MT_IB] ^ (s >> 1) ^ MAGIC(s);
      }
      
      s = TWIST(b, MT_LEN-1, 0);
      b[MT_LEN-1] = b[MT_IA-1] ^ (s >> 1) ^ MAGIC(s);
    }
  mt_index = idx + sizeof(unsigned_long_int);
  return *(unsigned_long_int *)((unsigned char *)b + idx);
  /*
    Matsumoto and Nishimura additionally confound the bits returned to the caller
    but this doesn't increase the randomness, and slows down the generator by
    as much as 25%.  So I omit these operations here.
    
    r ^= (r >> 11);
    r ^= (r << 7) & 0x9D2C5680;
    r ^= (r << 15) & 0xEFC60000;
    r ^= (r >> 18);
  */
}
Tip: Filter by directory path e.g. /media app.js to search for public/media/app.js.
Tip: Use camelCasing e.g. ProjME to search for ProjectModifiedEvent.java.
Tip: Filter by extension type e.g. /repo .js to search for all .js files in the /repo directory.
Tip: Separate your search with spaces e.g. /ssh pom.xml to search for src/ssh/pom.xml.
Tip: Use ↑ and ↓ arrow keys to navigate and return to view the file.
Tip: You can also navigate files with Ctrl+j (next) and Ctrl+k (previous) and view the file with Ctrl+o.
Tip: You can also navigate files with Alt+j (next) and Alt+k (previous) and view the file with Alt+o.