rand.c

#include "mathutils.h"

#define W (32)
#define N (624)
#define M (397)
#define R (31)
#define A (0x9908B0DF)
#define U (11)
#define D (0xFFFFFFFF)
#define S (7)
#define B (0x9D2C5680)
#define T (15)
#define C (0xEFC60000)
#define L (18)
#define F (1812433253)
#define SECONDARY_SEED	(4357)	//This is the seed prescribed in the Nishimura, Matsumoto C-code implementation
#define U_MASK (((1u << W)-1) & ~(1u << R) - 1) 	//Create W-R 1's followed by R 0's
#define L_MASK (~U_MASK) 	//Inversion of the upper bitmask


int x[N];
int i = N;
int seeded = 0;

/*************************************************************
 Below is an implementation of the pseudorandom number generator 
 scheme -known as the Mersenne Twister (MT)- proposed by Matsumoto
 and Nishimura (1998). Returns 32-bit integers uniformly at random
 in (0, 2^W-1). 
 The coefficients defined above are taken straight from the 
 1998 publication.
*************************************************************/

extern void seed_MT19937(int seed){
	if(seeded)
		{return;}
	seeded = 1;
	x[0] = seed;
	int j;
	for(j = 1; j < N; j++)
		{x[j] = (((1u << W)-1) & (F*(x[j-1] ^ (x[j-1] >> (W-2))) + j));}
	return;
}

extern int rand_MT19937(){
	int y;
	if(!seeded)
		{seed_MT19937(SECONDARY_SEED);}
	if(i >= N){		//Once indexer i rolls over N, twist and generate new randoms
		int k = 0;
		for(;k < N;k++){
			y = ((x[k] & U_MASK) | (x[(k+1)] & L_MASK));
			x[k] = x[(k+M) % N] ^ (y >> 1) ^ (((y & 0x1) == 0x0) ? 0 : A);	//Multiply by twist transformation matrix {0,A}
		}
		i = 0;
	}
	y = x[i];
	y = y ^ (y >> U);
	y = y ^ ((y << S) & B);
	y = y ^ ((y << T) & C);
	y = y ^ (y >> L);
	i += 1;
	return y;
}