mt19937ar.c

Go to the documentation of this file.

#include <stdio.h>
#include "potion.h"

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

static unsigned long state[N]; /* the array for the state vector  */
static int left = 1;
static int initf = 0;
static unsigned long *next;

/* initializes state[N] with a seed */
void init_genrand(unsigned long s) {
  int j;
  state[0]= s & 0xffffffffUL;
  for (j=1; j<N; j++) {
    state[j] = (1812433253UL * (state[j-1] ^ (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       */
    state[j] &= 0xffffffffUL;  /* for >32 bit machines */
  }
  left = 1; initf = 1;
}

/* initialize by an array with array-length */
/* init_key is the array for initializing keys */
/* key_length is its length */
/* slight change for C++, 2004/2/26 */
void init_by_array(unsigned long init_key[], int key_length) {
  int i, j, k;
  init_genrand(19650218UL);
  i=1; j=0;
  k = (N>key_length ? N : key_length);
  for (; k; k--) {
    state[i] = (state[i] ^ ((state[i-1] ^ (state[i-1] >> 30)) * 1664525UL))
      + init_key[j] + j; /* non linear */
    state[i] &= 0xffffffffUL; /* for WORDSIZE > 32 machines */
    i++; j++;
    if (i>=N) { state[0] = state[N-1]; i=1; }
    if (j>=key_length) j=0;
  }
  for (k=N-1; k; k--) {
    state[i] = (state[i] ^ ((state[i-1] ^ (state[i-1] >> 30)) * 1566083941UL))
      - i; /* non linear */
    state[i] &= 0xffffffffUL; /* for WORDSIZE > 32 machines */
    i++;
    if (i>=N) { state[0] = state[N-1]; i=1; }
  }

  state[0] = 0x80000000UL; /* MSB is 1; assuring non-zero initial array */ 
  left = 1; initf = 1;
}

static void next_state(void) {
  unsigned long *p=state;
  int j;

  /* if init_genrand() has not been called, */
  /* a default initial seed is used     */
  if (initf==0) init_genrand(5489UL);

  left = N;
  next = state;
  
  for (j=N-M+1; --j; p++) {
    /*@ assert Value: mem_access: \valid_read(p+1); */
    /*@ assert Value: mem_access: \valid_read(p+397); */
    *p = p[M] ^ TWIST(p[0], p[1]);
  }

  for (j=M; --j; p++) {
    /*@ assert Value: mem_access: \valid_read(p+1); */
    /*@ assert Value: mem_access: \valid_read(p+(int)(397-624)); */
    *p = p[M-N] ^ TWIST(p[0], p[1]);
  }

  /*@ assert Value: mem_access: \valid(p); */
  /*@ assert Value: mem_access: \valid_read(p+(int)(397-624)); */
  *p = p[M-N] ^ TWIST(p[0], state[0]);
}

unsigned long potion_rand_int(void) {
  unsigned long y;

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

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

  return y;
}

double potion_rand_double(void) {
  unsigned long a=potion_rand_int()>>5, b=potion_rand_int()>>6; 
  return(a*67108864.0+b)*(1.0/9007199254740992.0); 
} 
PN potion_srand(Potion *P, PN cl, PN self, PN seed) {
  init_genrand(PN_INT(seed));
  return self;
}

PN potion_rand(Potion *P, PN cl, PN self) {
  return PN_NUM(potion_rand_int());
}

PN potion_num_rand(Potion *P, PN cl, PN self) {
  return PN_NUM(potion_rand_double());
}