The Mersenne Twister implementation is now simpler, shorter and more
similar to the original C code.

Follows a standalone version of the original C code:

```c
#include <limits.h>
#include <sys/types.h>

typedef unsigned int uint32_t;

/* Period parameters */
#define N 624
#define M 397
#define MATRIX_A 0x9908b0dfU	/* constant vector a */
#define UMASK 0x80000000U	/* most significant w-r bits */
#define LMASK 0x7fffffffU	/* least significant r bits */
#define MIXBITS(u,v) ( ((u) & UMASK) | ((v) & LMASK) )
#define TWIST(u,v) ((MIXBITS((u),(v)) >> 1) ^ ((v)&1U ? MATRIX_A : 0U))

struct MT {
  /* assume int is enough to store 32bits */
  uint32_t state[N]; /* the array for the state vector  */
  uint32_t *next;
  int left;
};

/* initializes state[N] with a seed */
static void init_genrand(struct MT *mt, unsigned int s) {
  int j;
  mt->state[0] = s & 0xffffffffU;
  for (j=1; j<N; j++) {
    mt->state[j] = (1812433253U * (mt->state[j-1] ^ (mt->state[j-1] >> 30)) + j);
    /* See Knuth TAOCP Vol2. 3rd Ed. P.106 for multiplier. */
    /* In the previous versions, MSBs of the seed affect   */
    /* only MSBs of the array state[].                     */
    /* 2002/01/09 modified by Makoto Matsumoto             */
    mt->state[j] &= 0xffffffff;  /* for >32 bit machines */
  }
  mt->left = 1;
  mt->next = mt->state + N;
}

static void
  next_state(struct MT *mt)
{
  uint32_t *p = mt->state;
  int j;

  mt->left = N;
  mt->next = mt->state;

  for (j=N-M+1; --j; p++)
    *p = p[M] ^ TWIST(p[0], p[1]);

  for (j=M; --j; p++)
    *p = p[M-N] ^ TWIST(p[0], p[1]);

  *p = p[M-N] ^ TWIST(p[0], mt->state[0]);
}

/* generates a random number on [0,0xffffffff]-interval */
static unsigned int
  genrand_int32(struct MT *mt)
{
  /* mt must be initialized */
  unsigned int y;

  if (--mt->left <= 0) next_state(mt);
  y = *mt->next++;

  /* Tempering */
  y ^= (y >> 11);
  y ^= (y << 7) & 0x9d2c5680;
  y ^= (y << 15) & 0xefc60000;
  y ^= (y >> 18);

  return y;
}

static double
  int_pair_to_real_exclusive(uint32_t a, uint32_t b)
{
  a >>= 5;
  b >>= 6;
  return(a*67108864.0+b)*(1.0/9007199254740992.0);
}

static double
  genrand_real(struct MT *mt)
{
  /* mt must be initialized */
  unsigned int a = genrand_int32(mt), b = genrand_int32(mt);
  return int_pair_to_real_exclusive(a, b);
}


#include <stdio.h>

int main (int argc, char const *argv[]) {
  struct MT mt = { NULL, (uint32_t*)NULL, (int)NULL };
  init_genrand(&mt, 1234);
  int j = 10;
  while (j--) {
    printf("%.0f, ", (genrand_real(&mt) * 1001.0));
  }
  printf("-----\n");

  return 0;
}
```