|
|
Boost : |
From: jeetsukumaran_at_[hidden]
Date: 1999-11-05 10:25:11
Hi!
I've written a C++ implementation of the Mersenne Twister PRNG, and I
wish to submit it to the boost library. I believe the class conforms
to all requirements ...
When going through the submission guidelines, I read that this mailing
list is often used to pre-test code. So, here I present it below, for
review, feedback, criticism etc.
The code has been formatted to to avoid line wraps, and has been
un-namespaced.
I posted a version on USENET some time ago, and this version
incorporates all the fixes that were found to be needed on that
version. The code is under the GNU General Public License, which I
assume fulfills BOOST standards.
-- jeet
//--- MersenneTwister.h ------------------------------------------//
/**
* \file MersenneTwister.h
*
* Declaration file for MersenneTwister: MT19973B Psuedo-Random
* Number Generator.
*
* This code was based largely on a C implementation by Takuji
* Nishimura. The original C source code stated:
*
* Copyright (C) 1997 Makoto Matsumoto and Takuji Nishimura. When
* you use this, send an email to: matumoto_at_[hidden] with
* an appropriate reference to your work.
*
* LICENCE CONDITIONS:
* Matsumoto and Nishimura consent to GNU General Public Licence
* for this code.
*
* This C++ implementation source code is also released under the
* GNU General Public License:
*
* (C) Copyright 1999 Jeet Sukumaran.
*
* This program is free software; you can redistribute it and/or
* modify it under the terms of the GNU General Public License as
* published by the Free Software Foundation; either version 2 of
* the License, or (at your option) any later version.
*
* This program 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
* General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with this program; if not, write to the Free Software
* Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA
* 02111-1307 USA
*
* Version History:
*
*------------------
* Takuji Nishimura
* 1997:
*
* Original C source.
*
*----------------
* Jeet Sukumaran
* October 1999:
*
* Re-implemented the algorithm -- it is now encapsulated in a C++
* class that maintains an appropriate internal state, with an
* interface that allows for a variety of random number generation
* protocols.
*
*----------------
* Jeet Sukumaran
* November 3 1999:
*
* - MAJOR bug- and behavior-fixes: the twister works as advertised
* now; also an interface overhaul
*
* @author Jeet Sukumaran (jeetsukumaran_at_[hidden])
* @date November 1999
*/
#ifndef MERSENNETWISTER_H
#define MERSENNETWISTER_H
#include <vector>
#include <limits>
/**
* A C++ implementation of Makoto Matsumoto's and Takuji
* Nishimura's MT19937B "Mersenne Twister" pseudo-random number
* generator algorithm.
*
* This is a C++ implementation of the MT19937B, a "Twisted
* Generalized Feedback Shift Sequence Register" pseudo-random
* number generator,
*
* After being suitably seeded with a non-zero integer, this
* generator will give you a period of 2^19937 - 1 (!!!), and
* optimal equidistribution properties in up to 623 dimensions.
* It has passed the stringent "diehard" suite of statistical
* tests for randomness, and, as far as I can tell, is one of
* the better PRNG's around -- not perfect, but very good. Full
* details on the algorithm can be found at the algorithm's
* homepage:
*
* http://www.math.keio.ac.jp/matumoto/emt.html.
*
* Also see:
*
* M. Matsumoto and T. Nishimura. 1998. "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.
*
* @author Jeet Sukumaran (jeetsukumaran_at_[hidden])
* @date October 1999
*/
class MersenneTwister {
public:
// TYPES
typedef double RealType;
typedef unsigned long UIntType;
typedef signed long SIntType;
// BEHAVIORS
//////////////
// lifecycle
// constructs a self-seeded generator
MersenneTwister(void);
// constructs a user-seeded generator: seed must be non-zero
MersenneTwister(UIntType seed);
// copy constructor
MersenneTwister(const MersenneTwister& rhs);
// assignment
MersenneTwister& operator=(const MersenneTwister& rhs);
// our destructor does not do much
virtual ~MersenneTwister(void);
////////////////////////
// operation interface
// returns uniform random integer deviate in [0, 2^32-1]
// i,e, 0 <= n <= 2^32-1
virtual UIntType getInteger(void);
// returns uniform random integer deviate in [0, range)
// i.e. 0 <= n < range
virtual UIntType getInteger(UIntType range);
// returns uniform random real deviate in [0.0, 1.0]
// i.e. 0.0 <= n <= 1.0
virtual RealType getReal(void);
// returns uniform random real deviate in [-1, 1]
// i.e. -1.0 <= r <= 1.0
virtual RealType getNoise(void);
// true if getReal() returns a value that is less than prob
virtual bool checkProbability(RealType prob);
// to be used in STL algorithms such as std::generate_n
// similar to calling getInteger() -- 0 <= n <= 2^32-1
virtual UIntType operator()(void);
// to be used in STL algorithms such as std::random_shuffle
// similar to calling getInteger(n)
virtual UIntType operator()(UIntType n);
/**
* (experimental) Returns a random value of a (built-in)
* type that varies uniformly across the entire range of
* that type.
*
* @return a value of ResultType, n, where min <= n <= max;
* min and max are the implementation/architecture
* dependent minimum and maximum values for
* variables of ResultType.
*/
template <typename ResultType>
ResultType generate(void) {
ResultType mn = std::numeric_limits<ResultType>::min();
ResultType mx = std:: numeric_limits<ResultType>::max();
return
ResultType(mn+
((static_cast<double>(this->getInteger())/maxRand_)
* (mx-mn)));
}
protected:
// system set-up and seeding
virtual void seedGenerator(UIntType seed);
virtual UIntType autoSeed(void);
private:
// period parameters, algorithm constants, masks etc.
// defined in the implementation file
static const UIntType N_;
static const UIntType M_;
static const UIntType upperMask_;
static const UIntType lowerMask_;
static const UIntType maskB_;
static const UIntType maskC_;
static const UIntType mag01_[2];
static const UIntType maxRand_;
static const UIntType signedDiv_;
private:
// algorithm state variables
std::vector<UIntType> state_;
UIntType count_;
};
// MersenneTwister
////////////////////////////////////////////////////////////////
#endif // #ifndef MERSENNETWISTER_H
//--- MersenneTwister.cpp ----------------------------------------//
/**
* \file MersenneTwister.cpp
*
* Implementation file for MersenneTwister: MT19973B Psuedo-Random
* Number Generator.
*
* This code was based largely on a C implementation by Takuji
* Nishimura. The original C source code stated:
*
* Copyright (C) 1997 Makoto Matsumoto and Takuji Nishimura. When
* you use this, send an email to: matumoto_at_[hidden] with
* an appropriate reference to your work.
*
* LICENCE CONDITIONS:
* Matsumoto and Nishimura consent to GNU General Public Licence
* for this code.
*
* This C++ implementation source code is also released under the
* GNU General Public License:
*
* (C) Copyright 1999 Jeet Sukumaran.
*
* This program is free software; you can redistribute it and/or
* modify it under the terms of the GNU General Public License as
* published by the Free Software Foundation; either version 2 of
* the License, or (at your option) any later version.
*
* This program 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
* General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with this program; if not, write to the Free Software
* Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA
* 02111-1307 USA
*
* Version History:
*
*------------------
* Takuji Nishimura
* 1997:
*
* Original C source.
*
*----------------
* Jeet Sukumaran
* October 1999:
*
* Re-implemented the algorithm -- it is now encapsulated in a C++
* class that maintains an appropriate internal state, with an
* interface that allows for a variety of random number generation
* protocols.
*
*----------------
* Jeet Sukumaran
* November 3 1999:
*
* - MAJOR bug- and behavior-fixes: the twister works as advertised
* now; also an interface overhaul
*
* @author Jeet Sukumaran (jeetsukumaran_at_[hidden])
* @date November 1999
*/
// Standard C Library, for auto-seed generation
#include <ctime>
#include <cstdlib>
#include "MersenneTwister.h"
////////////////////////////////////////////////////////////////////
// period parameters, magic vectors, algorithm constants and masks
const MersenneTwister::UIntType
MersenneTwister::N_ = 624;
const MersenneTwister::UIntType
MersenneTwister::M_ = 397;
const MersenneTwister::UIntType
MersenneTwister::upperMask_ = 0x80000000;
const MersenneTwister::UIntType
MersenneTwister::lowerMask_ = 0x7FFFFFFF;
const MersenneTwister::UIntType
MersenneTwister::mag01_[2] = {0x0, 0x9908B0DF};
const MersenneTwister::UIntType
MersenneTwister::maskB_ = 0x9D2C5680;
const MersenneTwister::UIntType
MersenneTwister::maskC_ = 0xEFC60000;
const MersenneTwister::UIntType
MersenneTwister::maxRand_ = 0xFFFFFFFF;
const MersenneTwister::UIntType
MersenneTwister::signedDiv_ = 0x80000000;
////////////////////////////////////////////////////////////////////
// lifecycle & assignment
/**
* Default constructor
*
* Constructs an autoseeded generator.
*/
MersenneTwister::MersenneTwister(void) {
this->seedGenerator(this->autoSeed());
} // default constructor
/**
* Seeded constructor
*
* Constructs an generator with the provided non-zero seed.
*
* @param seed non-zero integer used to seed the generator
*/
MersenneTwister::MersenneTwister(UIntType seed) {
this->seedGenerator(seed);;
} // Seeded constructor
/**
* Copy constructor.
*
* The default memberwise copy should be quite sufficient. This
* is simply for form and as a developmental stub.
*
* @param rhs MersenneTwister object to be copied
*/
MersenneTwister::MersenneTwister(const MersenneTwister& rhs)
: state_(rhs.state_),
count_(rhs.count_) {
} // copy constructor
/**
* Assignment operator.
*
* Again, the default works just as well.
*
* @param rhs MersenneTwister object to be assigned to this
* @return this object reference
*/
MersenneTwister&
MersenneTwister::operator=(const MersenneTwister& rhs) {
state_ = rhs.state_;
count_ = rhs.count_;
return (*this);
} // assignment
/**
* Destructor.
*
* Uneccessary, but provided on the "good citizen" principle.
*/
MersenneTwister::~MersenneTwister(void) {
} // destructor
////////////////////////////////////////////////////////////////////
// initialization, seeding & setup
/**
* Uses the given seed, WHICH MUST BE A NON-ZERO INTEGER,
* to set the initial state of the generator. If the seed
* given is zero, it is ignored and the program selects
* another seed by itself. This also sets the state vector
* index to the correct initial value.
*
* @param seed non-zero integer to seed the generator
*/
void MersenneTwister::seedGenerator(UIntType seed) {
// check for 0 seed, and change if neccesary
while (!seed) {
seed = this->autoSeed();
}
state_.push_back(seed & 0xFFFFFFFF); // seed[0]
for (unsigned int k = 1; k < N_; ++k) {
this->state_.push_back((69069 * this->state_[k-1])
& 0xffffffff);
}
this->count_ = 1;
} // seedGenerator
/**
* Generates a non-zero integer value to be used to
* seed the generator.
*
* Currently, this function uses the tried and tired "poor
* man's" random-number-generator -- i.e. the method I used
* until I learnt of the ills of the implementation prng.
* I can't, however, think of a demonstatably better way to
* do this, and, in any case, imagine that it (hopefully)
* won't make too much of a difference -- as only one value
* is retrieved from this method at a time for any given
* seed.
*
* @return non-zero integer value
*/
MersenneTwister::UIntType
MersenneTwister::autoSeed(void) {
// I know, I know ... but what else to do?
std::srand(std::time(NULL));
return std::rand();
} // autoSeed
////////////////////////////////////////////////////////////////////
// random number protocols
/**
* Algorithm kernel -- the Twister itself: pulls the next
* 32-bit integer from the sequence.
*
* This is the heart and raison d'etre of the class --
* the MT19937B Mersenne Twister algorithm. On each call, it
* pulls the next number from the MT19937B sequence, and
* modifies the internal state appropriately. It is invoked
* by all other random number generation request methods
* provided for in the class interface. The number returned
* is a uniform integer deviate in [0, 2^32-1].
*
* @return an integer value, i, where 0 <= i <= 2^32-1
*/
MersenneTwister::UIntType
MersenneTwister::getInteger(void) {
// generate N_ words at time
if (this->count_ == N_) {
// reload the generator
// this could (SHOULD) be a separate function ...
unsigned int kk;
UIntType w;
for (kk = 0; kk < N_ - M_; ++kk) {
w = (this->state_[kk] & upperMask_)
| (this->state_[kk+1] & lowerMask_);
this->state_[kk] = this->state_[kk + M_]
^ (w >> 1)
^ mag01_[w & 0x1];
}
for ( ; kk < N_ - 1; ++kk) {
w = (this->state_[kk] & upperMask_)
| (state_[kk+1] & lowerMask_);
this->state_[kk]= this->state_[kk + (M_ - N_)]
^ (w >> 1)
^ mag01_[w & 0x1];
}
w = (this->state_[N_-1] & upperMask_)
| (this->state_[0] & lowerMask_);
this->state_[N_-1] = this->state_[M_ - 1]
^ (w >> 1)
^ mag01_[w & 0x1];
this->count_ = 0;
// huh?
} // if count_ == N_
// tempering: hard-coded magic numbers directly
// from the original ...
UIntType y; // facilitate return value optimization
y = this->state_[this->count_++];
y ^= (y >> 11);
y ^= (y << 7) & maskB_;
y ^= (y << 15) & maskC_;
y &= 0xffffffff;
y ^= (y >> 18);
return y;
} // getInteger
/**
* Returns an integer between 0 (inclusive) and the specified
* limit (exclusive).
*
* This returns a uniform integer deviate in [0, range), i.e.
* from zero up to, but not including the value requested.
*
* @param range upper-limit for return value.
* @return an integer value, i, where 0 <= i < range
*/
MersenneTwister::UIntType
MersenneTwister::getInteger(UIntType range) {
// There is a 1 in 2^32-1 chance that getInteger() will
// return 0xffffffff, in which case our algorithm will
// return "range"; this lies outside the promised [0, range)
// bounds.
// To guard against this, we check for (and reject)
// a 0xffffffff value ...
UIntType i;
while ( (i = this->getInteger()) == maxRand_ ); // check
return UIntType((static_cast<RealType>(i)/maxRand_) * range);
} // getInteger
/**
* Returns a random real (floating-point) number between
* 0 and 1, inclusive.
*
* This uses the MT generator to provide a uniform real
* deviate in [0.0, 1.0], i.e. between 0 and 1, inclusive
* of the endpoints.
*
* @return a real value, r, where 0.0 <= r <= 1.0
*/
MersenneTwister::RealType
MersenneTwister::getReal(void) {
// Convert the integer to a real and return:
return static_cast<RealType>(this->getInteger())/maxRand_;
} // getReal
/**
* Returns a random real (floating-point) number between
* -1 and 1 (inclusive).
*
* This uses the MT generator to provide a uniform real
* deviate in [-1, 1], i.e. between -1 and 1, inclusive of the
* endpoints.
*
* @return real value, r, where -1 <= r <= 1
*/
MersenneTwister::RealType
MersenneTwister::getNoise(void) {
return
static_cast<RealType>(static_cast<SIntType>(this->getInteger()))
/signedDiv_;
} // getNoise
/**
* Evaluates to <code>true</code> or <code>false</code> with a
* probability of <code>prob</code>.
*
* Calls getReal() and compares result with <code>prob</code>.
* If the result is less than <code>prob</code>, then this
* function evaluates to <code>true</code>. Otherwise,
* <code>false</code>.
*
* @param prob probability of <code>true</code>
* @return <code>true</code> or <code>false</code>
* with a probability of <code>prob</code> and
* <code>1 - prob</code> each respectively.
*/
bool MersenneTwister::checkProbability(RealType prob) {
return this->getReal() < prob;
} // checkProbability
/**
* Returns a uniform 32-bit integer deviate in [0, 2^32-1].
*
* Basically duplicates the functionality of
* <code>getInteger()</code>, but with the specific signature
* required to be used in some STL algorithms, such as
* <code>std::generate_n</code>.
*
* @return an integer value, i, where 0 <= i <= 2^32-1
*/
MersenneTwister::UIntType
MersenneTwister::operator()(void) {
return UIntType(this->getInteger());
} // operator()
/**
* Returns a uniform 32-bit integer deviate in [0, n).
*
* Basically duplicates the functionality of
* <code>getInteger(n)</code>, but with the specific signature
* required to be used in some STL algorithms, such as
* <code>std::random_shuffle</code>.
*
* @return an integer value, i, where 0 <= i < n
*/
MersenneTwister::UIntType
MersenneTwister::operator()(UIntType n) {
// there is a 1 in 2^32-1 chance that getInteger() will
// return 0xffffffff, in which case our algorithm will
// return "n"; this lies outside the promised [0, n) range.
// To guard against this, we check for (and reject)
// a 0xffffffff value ...
UIntType i;
while ( (i = this->getInteger()) == maxRand_ ); // check
return UIntType( (static_cast<RealType>(i)/maxRand_) * n );
} // operator(n)
Boost list run by bdawes at acm.org, gregod at cs.rpi.edu, cpdaniel at pacbell.net, john at johnmaddock.co.uk