/* A C-program for MT19937: Real number version([0,1)-interval) */
/* (1999/10/28) */
/* genrand() generates one pseudorandom real number (double) */
/* which is uniformly distributed on [0,1)-interval, for each */
/* call. sgenrand(seed) sets initial values to the working area */
/* of 624 words. Before genrand(), sgenrand(seed) must be */
/* called once. (seed is any 32-bit integer.) */
/* Integer generator is obtained by modifying two lines. */
/* Coded by Takuji Nishimura, considering the suggestions by */
/* Topher Cooper and Marc Rieffel in July-Aug. 1997. */
/* This library is free software; you can redistribute it and/or */
/* modify it under the terms of the GNU Library General Public */
/* License as published by the Free Software Foundation; either */
/* version 2 of the License, or (at your option) any later */
/* version. */
/* This library is distributed in the hope that it will be useful, */
/* but WITHOUT ANY WARRANTY; without even the implied warranty of */
/* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. */
/* See the GNU Library General Public License for more details. */
/* You should have received a copy of the GNU Library General */
/* Public License along with this library; if not, write to the */
/* Free Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA */
/* 02111-1307 USA */
/* Copyright (C) 1997, 1999 Makoto Matsumoto and Takuji Nishimura. */
/* Any feedback is very welcome. For any question, comments, */
/* see http://www.math.keio.ac.jp/matumoto/emt.html or email */
/* matumoto@math.keio.ac.jp */
/* REFERENCE */
/* M. Matsumoto and T. Nishimura, */
/* "Mersenne Twister: A 623-Dimensionally Equidistributed Uniform */
/* Pseudo-Random Number Generator", */
/* ACM Transactions on Modeling and Computer Simulation, */
/* Vol. 8, No. 1, January 1998, pp 3--30. */
/* Adapted by Martin Reinsprecht and Raphael Langerhorst */
/* for the G System */
#include "PseudoRNG.h"
namespace GBE
{
namespace Util
{
// Standard Constructor
PseudoRNG::PseudoRNG()
{}
// Standard Destructor
PseudoRNG::~PseudoRNG()
{}
double PseudoRNG::getNumberDouble()
{
return this->genrandDouble();
}
unsigned long PseudoRNG::getNumberInt()
{
return this->genrandInt();
}
void PseudoRNG::setNumber( unsigned long j )
{
this->sgenrand( j );
}
/* Initializing the array with a seed */
void PseudoRNG::sgenrand( unsigned long seed )
{
int i;
for ( i = 0;i < N;i++ )
{
mt[ i ] = seed & 0xffff0000;
seed = 69069 * seed + 1;
mt[ i ] |= ( seed & 0xffff0000 ) >> 16;
seed = 69069 * seed + 1;
}
mti = N;
}
/* Initialization by "sgenrand()" is an example. Theoretically, */
/* there are 2^19937-1 possible states as an intial state. */
/* This function allows to choose any of 2^19937-1 ones. */
/* Essential bits in "seed_array[]" is following 19937 bits: */
/* (seed_array[0]&UPPER_MASK), seed_array[1], ..., seed_array[N-1]. */
/* (seed_array[0]&LOWER_MASK) is discarded. */
/* Theoretically, */
/* (seed_array[0]&UPPER_MASK), seed_array[1], ..., seed_array[N-1] */
/* can take any values except all zeros. */
void PseudoRNG::lsgenrand( unsigned long seed_array[] )
{
/* the length of seed_array[] must be at least N */
int i;
for ( i = 0;i < N;i++ )
mt[ i ] = seed_array[ i ];
mti = N;
}
/* generating reals */
// alte Form
//double PseudoRNG::genrand(unsigned long y,unsigned long mag01[2])
double PseudoRNG::genrandDouble()
{
return ( ( double ) genrandInt() * 2.3283064365386963e-10 ); /* reals: [0,1)-interval */
/* return y; */ /* for integer generation */
}
/* unsigned long */
// alte Form
//double PseudoRNG::genrand(unsigned long y,unsigned long mag01[2])
unsigned long PseudoRNG::genrandInt()
{
unsigned long mag01[ 2 ];
unsigned long y;
/* mag01[x] = x * MATRIX_A for x=0,1 */
mag01[ 0 ] = 0x0;
mag01[ 1 ] = MATRIX_A;
if ( mti >= N )
{ /* generate N words at one time */
int kk;
if ( mti == N + 1 ) /* if sgenrand() has not been called, */
sgenrand( 4357 ); /* a default initial seed is used */
for ( kk = 0;kk < N - M;kk++ )
{
y = ( mt[ kk ] & UPPER_MASK ) | ( mt[ kk + 1 ] & LOWER_MASK );
mt[ kk ] = mt[ kk + M ] ^ ( y >> 1 ) ^ mag01[ y & 0x1 ];
}
for ( ;kk < N - 1;kk++ )
{
y = ( mt[ kk ] & UPPER_MASK ) | ( mt[ kk + 1 ] & LOWER_MASK );
mt[ kk ] = mt[ kk + ( M - N ) ] ^ ( y >> 1 ) ^ mag01[ y & 0x1 ];
}
y = ( mt[ N - 1 ] & UPPER_MASK ) | ( mt[ 0 ] & LOWER_MASK );
mt[ N - 1 ] = mt[ M - 1 ] ^ ( y >> 1 ) ^ mag01[ y & 0x1 ];
mti = 0;
}
y = mt[ mti++ ];
y ^= TEMPERING_SHIFT_U( y );
y ^= TEMPERING_SHIFT_S( y ) & TEMPERING_MASK_B;
y ^= TEMPERING_SHIFT_T( y ) & TEMPERING_MASK_C;
y ^= TEMPERING_SHIFT_L( y );
return y;
}
/* This main() outputs first 1000 generated numbers. */
// we need no main here :-)
/*
main()
{
int i;
sgenrand(4357);
for (i=0; i<1000; i++) {
printf("%10.8f ", genrand());
if (i%5==4) printf("\n");
}
}*/
}
}
syntax highlighted by Code2HTML, v. 0.9.1