
#include "rand_mersenne.h"


/******************************************************************************/
/*                             VARIABLES GLOBALES                             */
/******************************************************************************/

unsigned long mt[MERS_N];    /* vecteur d'etat         */
int mti;                     /* indice dans le vecteur */
architecture_t architecture; /* type d'architecture    */


/******************************************************************************/
/*                   INITIALISATION DU GENERATEUR MERSENNE                    */
/******************************************************************************/

void rand_init(unsigned long graine) {

  /* variables locales */

  conversion_t conversion;


  /* graine du generateur */

  mt[0]=graine;
  for (mti=1;mti<MERS_N;mti++) {
    mt[mti]=(1812433253UL*(mt[mti-1]^(mt[mti-1]>>30))+mti);
  }
  
  
  /* determination du type d'architecture */

  conversion.f=1.;
  if      (conversion.i[1]==0x3FF00000) architecture=LITTLEENDIAN;
  else if (conversion.i[0]==0x3FF00000) architecture=BIGENDIAN;
  else                                  architecture=NONIEEE;
}


/******************************************************************************/
/*          TIRAGE D'UN DOUBLE DANS [0,1[ PAR LE GENERATEUR MERSENNE          */
/******************************************************************************/

double randu(void) {

  /* variables locales */

  conversion_t conversion;
  unsigned long y;
  static unsigned long mag01[2]={0,MERS_A};
  int kk;


  /* tirage aleatoire de 32 bits */

  if (mti>=MERS_N) {
    for (kk=0;kk<MERS_N-MERS_M;kk++) {
      y=(mt[kk]&UPPER_MASK)|(mt[kk+1]&LOWER_MASK);
      mt[kk]=mt[kk+MERS_M]^(y>>1)^mag01[y&1];
    }
    for (;kk<MERS_N-1;kk++) {
      y=(mt[kk]&UPPER_MASK)|(mt[kk+1]&LOWER_MASK);
      mt[kk]=mt[kk+(MERS_M-MERS_N)]^(y>>1)^mag01[y&1];
    }
    y=(mt[MERS_N-1]&UPPER_MASK)|(mt[0]&LOWER_MASK);
    mt[MERS_N-1]=mt[MERS_M-1]^(y>>1)^mag01[y&1];
    mti=0;
  }
  y=mt[mti++];
  y^= y>>MERS_U;
  y^=(y<<MERS_S)&MERS_B;
  y^=(y<<MERS_T)&MERS_C;
  y^= y>>MERS_L;


  /* conversion des bits generes en double dans [0,1[ */

  switch (architecture) {
    case LITTLEENDIAN:
      conversion.i[0]= y<<20;
      conversion.i[1]=(y>>12)|0x3FF00000;
      return conversion.f-1.;
    case BIGENDIAN:
      conversion.i[1]= y<<20;
      conversion.i[0]=(y>>12)|0x3FF00000;
      return conversion.f-1.;
    case NONIEEE:
    default:
      break;
  }
  return (double)y*(1./((double)(unsigned long)(-1L)+1.));
}

