ViewVC Help

View File | Revision Log | Show Annotations | View Changeset | Root Listing

root/OpenMD/branches/development/src/math/MersenneTwister.hpp

Comparing:
trunk/src/math/MersenneTwister.hpp (file contents), Revision 385 by tim, Tue Mar 1 20:10:14 2005 UTC vs.
branches/development/src/math/MersenneTwister.hpp (file contents), Revision 1767 by gezelter, Fri Jul 6 22:01:58 2012 UTC

#	Line 67 | Line 67
67		#include <time.h>
68		#include <math.h>
69		#include <vector>
70	<	namespace oopse {
70	>	namespace OpenMD {
71
72	<	class MTRand {
73	<	// Data
74	<	public:
75	<	typedef unsigned long uint32;  // unsigned integer type, at least 32 bits
72	>	class MTRand {
73	>	// Data
74	>	public:
75	>	typedef unsigned long uint32;  // unsigned integer type, at least 32 bits
76
77	<	enum { N = 624 };       // length of state vector
78	<	enum { SAVE = N + 1 };  // length of array for save()
77	>	enum { N = 624 };       // length of state vector
78	>	enum { SAVE = N + 1 };  // length of array for save()
79
80	<	private:
81	<	enum { M = 397 };  // period parameter
80	>	private:
81	>	enum { M = 397 };  // period parameter
82
83	<	uint32 state[N];   // internal state
84	<	uint32 *pNext;     // next value to get from state
85	<	int left;          // number of values left before reload needed
86	<	int nstrides_;
87	<	int stride_;
83	>	uint32 state[N];   // internal state
84	>	uint32 *pNext;     // next value to get from state
85	>	int left;          // number of values left before reload needed
86	>	int nstrides_;
87	>	int stride_;
88
89	<	//Methods
90	<	public:
91	<	MTRand( const uint32& oneSeed, int nstrides = 1, int stride = 0);  // initialize with a simple uint32
92	<	MTRand( uint32 *const bigSeed, uint32 const seedLength = N, int nstrides = 1, int stride = 0);  // or an array
93	<	MTRand(int nstrides = 1, int stride = 0);  // auto-initialize with /dev/urandom or time() and clock()
89	>	//Methods
90	>	public:
91	>	MTRand( const uint32& oneSeed, int nstrides, int stride);  // initialize with a simple uint32
92	>	MTRand( uint32 *const bigSeed, uint32 const seedLength, int nstrides, int stride);  // or an array
93	>	MTRand(int nstrides, int stride);  // auto-initialize with /dev/urandom or time() and clock()
94
95	<	// Do NOT use for CRYPTOGRAPHY without securely hashing several returned
96	<	// values together, otherwise the generator state can be learned after
97	<	// reading 624 consecutive values.
95	>	// Do NOT use for CRYPTOGRAPHY without securely hashing several returned
96	>	// values together, otherwise the generator state can be learned after
97	>	// reading 624 consecutive values.
98
99	<	// Access to 32-bit random numbers
100	<	double rand();                          // real number in [0,1]
101	<	double rand( const double& n );         // real number in [0,n]
102	<	double randExc();                       // real number in [0,1)
103	<	double randExc( const double& n );      // real number in [0,n)
104	<	double randDblExc();                    // real number in (0,1)
105	<	double randDblExc( const double& n );   // real number in (0,n)
106	<	uint32 randInt();                       // integer in [0,2^32-1] (modified for striding)
107	<	uint32 rawRandInt();                    // original randInt
108	<	uint32 randInt( const uint32& n );      // integer in [0,n] for n < 2^32
109	<	double operator()() { return rand(); }  // same as rand()
99	>	// Access to 32-bit random numbers
100	>	RealType rand();                          // real number in [0,1]
101	>	RealType rand( const RealType& n );         // real number in [0,n]
102	>	RealType randExc();                       // real number in [0,1)
103	>	RealType randExc( const RealType& n );      // real number in [0,n)
104	>	RealType randDblExc();                    // real number in (0,1)
105	>	RealType randDblExc( const RealType& n );   // real number in (0,n)
106	>	uint32 randInt();                       // integer in [0,2^32-1] (modified for striding)
107	>	uint32 rawRandInt();                    // original randInt
108	>	uint32 randInt( const uint32& n );      // integer in [0,n] for n < 2^32
109	>	RealType operator()() { return rand(); }  // same as rand()
110
111	<	// Access to 53-bit random numbers (capacity of IEEE double precision)
112	<	double rand53();  // real number in [0,1)
111	>	// Access to 53-bit random numbers (capacity of IEEE RealType precision)
112	>	RealType rand53();  // real number in [0,1)
113
114	<	// Access to nonuniform random number distributions
115	<	double randNorm( const double& mean = 0.0, const double& variance = 0.0 );
114	>	// Access to nonuniform random number distributions
115	>	RealType randNorm( const RealType mean = 0.0, const RealType variance = 0.0 );
116
117	<	// Re-seeding functions with same behavior as initializers
118	<	void seed( const uint32 oneSeed );
119	<	void seed( uint32 *const bigSeed, const uint32 seedLength = N );
120	<	void seed();
117	>	// Re-seeding functions with same behavior as initializers
118	>	void seed( const uint32 oneSeed );
119	>	void seed( uint32 *const bigSeed, const uint32 seedLength = N );
120	>	void seed();
121
122	<	std::vector<uint32>generateSeeds();
122	>	std::vector<uint32>generateSeeds();
123
124	<	// Saving and loading generator state
125	<	void save( uint32* saveArray ) const;  // to array of size SAVE
126	<	void load( uint32 *const loadArray );  // from such array
127	<	friend std::ostream& operator<<( std::ostream& os, const MTRand& mtrand );
128	<	friend std::istream& operator>>( std::istream& is, MTRand& mtrand );
124	>	// Saving and loading generator state
125	>	void save( uint32* saveArray ) const;  // to array of size SAVE
126	>	void load( uint32 *const loadArray );  // from such array
127	>	friend std::ostream& operator<<( std::ostream& os, const MTRand& mtrand );
128	>	friend std::istream& operator>>( std::istream& is, MTRand& mtrand );
129
130	<	protected:
131	<	void initialize( const uint32 oneSeed );
132	<	void reload();
133	<	uint32 hiBit( const uint32& u ) const { return u & 0x80000000UL; }
134	<	uint32 loBit( const uint32& u ) const { return u & 0x00000001UL; }
135	<	uint32 loBits( const uint32& u ) const { return u & 0x7fffffffUL; }
136	<	uint32 mixBits( const uint32& u, const uint32& v ) const
137	<	{ return hiBit(u) | loBits(v); }
138	<	uint32 twist( const uint32& m, const uint32& s0, const uint32& s1 ) const
139	<	{ return m ^ (mixBits(s0,s1)>>1) ^ (-loBit(s1) & 0x9908b0dfUL); }
140	<	static uint32 hash( time_t t, clock_t c );
141	<	};
130	>	protected:
131	>	void initialize( const uint32 oneSeed );
132	>	void reload();
133	>	uint32 hiBit( const uint32& u ) const { return u & 0x80000000UL; }
134	>	uint32 loBit( const uint32& u ) const { return u & 0x00000001UL; }
135	>	uint32 loBits( const uint32& u ) const { return u & 0x7fffffffUL; }
136	>	uint32 mixBits( const uint32& u, const uint32& v ) const
137	>	{ return hiBit(u) | loBits(v); }
138	>	#ifdef _MSC_VER
139	>	#pragma warning( push ) // save current warning settings
140	>	#pragma warning( disable : 4146 ) // warning C4146: unary minus operator applied to unsigned type, result still unsigned
141	>	#endif
142	>	uint32 twist( const uint32& m, const uint32& s0, const uint32& s1 ) const
143	>	{ return m ^ (mixBits(s0,s1)>>1) ^ (-loBit(s1) & 0x9908b0dfUL); }
144	>	#ifdef _MSC_VER
145	>	#pragma warning( pop ) // return warning settings to what they were
146	>	#endif
147
148	+	static uint32 hash( time_t t, clock_t c );
149	+	};
150
151	<	inline MTRand::MTRand( const uint32& oneSeed, int nstrides, int stride) : nstrides_(nstrides), stride_(stride) {
151	>
152	>	inline MTRand::MTRand( const uint32& oneSeed, int nstrides, int stride) : nstrides_(nstrides), stride_(stride) {
153		assert(stride_ < nstrides_ && stride_ >= 0);
154		seed(oneSeed);
155	<	}
155	>	}
156
157	<	inline MTRand::MTRand( uint32 *const bigSeed, const uint32 seedLength, int nstrides, int stride) : nstrides_(nstrides), stride_(stride) {
157	>	inline MTRand::MTRand( uint32 *const bigSeed, const uint32 seedLength, int nstrides, int stride) : nstrides_(nstrides), stride_(stride) {
158		assert(stride_ < nstrides_ && stride_ >= 0);
159		seed(bigSeed,seedLength);
160	<	}
160	>	}
161
162	<	inline MTRand::MTRand(int nstrides, int stride) : nstrides_(nstrides), stride_(stride){
162	>	inline MTRand::MTRand(int nstrides, int stride)       : nstrides_(nstrides), stride_(stride){
163		assert(stride_ < nstrides_ && stride_ >= 0);
164		seed();
165	<	}
165	>	}
166
167	<	inline double MTRand::rand()
168	<	{ return double(randInt()) * (1.0/4294967295.0); }
167	>	inline RealType MTRand::rand()
168	>	{ return RealType(randInt()) * (1.0/4294967295.0); }
169
170	<	inline double MTRand::rand( const double& n )
171	<	{ return rand() * n; }
170	>	inline RealType MTRand::rand( const RealType& n )
171	>	{ return rand() * n; }
172
173	<	inline double MTRand::randExc()
174	<	{ return double(randInt()) * (1.0/4294967296.0); }
173	>	inline RealType MTRand::randExc()
174	>	{ return RealType(randInt()) * (1.0/4294967296.0); }
175
176	<	inline double MTRand::randExc( const double& n )
177	<	{ return randExc() * n; }
176	>	inline RealType MTRand::randExc( const RealType& n )
177	>	{ return randExc() * n; }
178
179	<	inline double MTRand::randDblExc()
180	<	{ return ( double(randInt()) + 0.5 ) * (1.0/4294967296.0); }
179	>	inline RealType MTRand::randDblExc()
180	>	{ return ( RealType(randInt()) + 0.5 ) * (1.0/4294967296.0); }
181
182	<	inline double MTRand::randDblExc( const double& n )
183	<	{ return randDblExc() * n; }
176	<
177	<	inline double MTRand::rand53()
178	<	{
179	<	uint32 a = randInt() >> 5, b = randInt() >> 6;
180	<	return ( a * 67108864.0 + b ) * (1.0/9007199254740992.0);  // by Isaku Wada
181	<	}
182	>	inline RealType MTRand::randDblExc( const RealType& n )
183	>	{ return randDblExc() * n; }
184
185	<	inline double MTRand::randNorm( const double& mean, const double& variance )
186	<	{
187	<	// Return a real number from a normal (Gaussian) distribution with given
188	<	// mean and variance by Box-Muller method
189	<	double r = sqrt( -2.0 * log( 1.0-randDblExc()) ) * variance;
188	<	double phi = 2.0 * 3.14159265358979323846264338328 * randExc();
189	<	return mean + r * cos(phi);
190	<	}
185	>	inline RealType MTRand::rand53()
186	>	{
187	>	uint32 a = randInt() >> 5, b = randInt() >> 6;
188	>	return ( a * 67108864.0 + b ) * (1.0/9007199254740992.0);  // by Isaku Wada
189	>	}
190
191	<	/**
192	<	* This function is modified from the original to allow for random
193	<	* streams on parallel jobs.  It now takes numbers from by striding
194	<	* through the random stream and picking up only one of the random
195	<	* numbers per nstrides_.  The number it picks is the stride_'th
196	<	* number in the stride sequence.
197	<	*/
198	<	inline MTRand::uint32 MTRand::randInt() {
200	<
201	<	std::vector<uint32> ranNums(nstrides_);
202	<
203	<	for (int i = 0; i < nstrides_; ++i) {
204	<	ranNums[i] = rawRandInt();
191	>	inline RealType MTRand::randNorm( const RealType mean, const RealType variance )
192	>	{
193	>	// Return a real number from a normal (Gaussian) distribution with given
194	>	// mean and variance by Box-Muller method
195	>	assert(variance > 0);
196	>	RealType r = sqrt( -2.0 * log( 1.0-randDblExc()) * variance);
197	>	RealType phi = 2.0 * 3.14159265358979323846264338328 * randExc();
198	>	return mean + r * cos(phi);
199		}
206	–
207	–	return ranNums[stride_];
208	–	}
200
201	<	/**
202	<	* This is the original randInt function which implements the mersenne
203	<	* twister.
204	<	*/
205	<	inline MTRand::uint32 MTRand::rawRandInt()
206	<	{
207	<	// Pull a 32-bit integer from the generator state
208	<	// Every other access function simply transforms the numbers extracted here
218	<
219	<	if( left == 0 ) reload();
220	<	--left;
221	<
222	<	register uint32 s1;
223	<	s1 = *pNext++;
224	<	s1 ^= (s1 >> 11);
225	<	s1 ^= (s1 <<  7) & 0x9d2c5680UL;
226	<	s1 ^= (s1 << 15) & 0xefc60000UL;
227	<	return ( s1 ^ (s1 >> 18) );
228	<	}
201	>	/**
202	>	* This function is modified from the original to allow for random
203	>	* streams on parallel jobs.  It now takes numbers from by striding
204	>	* through the random stream and picking up only one of the random
205	>	* numbers per nstrides_.  The number it picks is the stride_'th
206	>	* number in the stride sequence.
207	>	*/
208	>	inline MTRand::uint32 MTRand::randInt() {
209
210	<	inline MTRand::uint32 MTRand::randInt( const uint32& n )
211	<	{
212	<	// Find which bits are used in n
213	<	// Optimized by Magnus Jonsson (magnus@smartelectronix.com)
214	<	uint32 used = n;
215	<	used |= used >> 1;
216	<	used |= used >> 2;
237	<	used |= used >> 4;
238	<	used |= used >> 8;
239	<	used |= used >> 16;
240	<
241	<	// Draw numbers until one is found in [0,n]
242	<	uint32 i;
243	<	do
244	<	i = randInt() & used;  // toss unused bits to shorten search
245	<	while( i > n );
246	<	return i;
247	<	}
248	<
249	<
250	<	inline void MTRand::seed( const uint32 oneSeed )
251	<	{
252	<	// Seed the generator with a simple uint32
253	<	initialize(oneSeed);
254	<	reload();
255	<	}
256	<
257	<
258	<	inline void MTRand::seed( uint32 *const bigSeed, const uint32 seedLength )
259	<	{
260	<	// Seed the generator with an array of uint32's
261	<	// There are 2^19937-1 possible initial states.  This function allows
262	<	// all of those to be accessed by providing at least 19937 bits (with a
263	<	// default seed length of N = 624 uint32's).  Any bits above the lower 32
264	<	// in each element are discarded.
265	<	// Just call seed() if you want to get array from /dev/urandom
266	<	initialize(19650218UL);
267	<	register int i = 1;
268	<	register uint32 j = 0;
269	<	register int k = ( N > seedLength ? N : seedLength );
270	<	for( ; k; --k )
271	<	{
272	<	state[i] =
273	<	state[i] ^ ( (state[i-1] ^ (state[i-1] >> 30)) * 1664525UL );
274	<	state[i] += ( bigSeed[j] & 0xffffffffUL ) + j;
275	<	state[i] &= 0xffffffffUL;
276	<	++i;  ++j;
277	<	if( i >= N ) { state[0] = state[N-1];  i = 1; }
278	<	if( j >= seedLength ) j = 0;
279	<	}
280	<	for( k = N - 1; k; --k )
281	<	{
282	<	state[i] =
283	<	state[i] ^ ( (state[i-1] ^ (state[i-1] >> 30)) * 1566083941UL );
284	<	state[i] -= i;
285	<	state[i] &= 0xffffffffUL;
286	<	++i;
287	<	if( i >= N ) { state[0] = state[N-1];  i = 1; }
288	<	}
289	<	state[0] = 0x80000000UL;  // MSB is 1, assuring non-zero initial array
290	<	reload();
291	<	}
292	<
293	<
294	<	inline void MTRand::seed()
295	<	{
296	<	std::vector<uint32> seeds;
297	<
298	<	seeds = generateSeeds();
299	<
300	<	if (seeds.size() == 1) {
301	<	seed( seeds[0] );
302	<	} else {
303	<	seed( &seeds[0], seeds.size() );
210	>	std::vector<uint32> ranNums(nstrides_);
211	>
212	>	for (int i = 0; i < nstrides_; ++i) {
213	>	ranNums[i] = rawRandInt();
214	>	}
215	>
216	>	return ranNums[stride_];
217		}
305	–	}
218
219	<
220	<	inline std::vector<MTRand::uint32> MTRand::generateSeeds() {
221	<	// Seed the generator with an array from /dev/urandom if available
222	<	// Otherwise use a hash of time() and clock() values
223	<
224	<	std::vector<uint32> bigSeed;
219	>	/**
220	>	* This is the original randInt function which implements the mersenne
221	>	* twister.
222	>	*/
223	>	inline MTRand::uint32 MTRand::rawRandInt()
224	>	{
225	>	// Pull a 32-bit integer from the generator state
226	>	// Every other access function simply transforms the numbers extracted here
227	>
228	>	if( left == 0 ) reload();
229	>	--left;
230	>
231	>	register uint32 s1;
232	>	s1 = *pNext++;
233	>	s1 ^= (s1 >> 11);
234	>	s1 ^= (s1 <<  7) & 0x9d2c5680UL;
235	>	s1 ^= (s1 << 15) & 0xefc60000UL;
236	>	return ( s1 ^ (s1 >> 18) );
237	>	}
238
239	<	// First try getting an array from /dev/urandom
240	<	FILE* urandom = fopen( "/dev/urandom", "rb" );
241	<	if( urandom )
242	<	{
243	<	bigSeed.resize(N);
244	<	register uint32 *s = &bigSeed[0];
245	<	register int i = N;
246	<	register bool success = true;
247	<	while( success && i-- )
248	<	success = fread( s++, sizeof(uint32), 1, urandom );
249	<	fclose(urandom);
250	<	if( success ) { return bigSeed; }
239	>	inline MTRand::uint32 MTRand::randInt( const uint32& n )
240	>	{
241	>	// Find which bits are used in n
242	>	// Optimized by Magnus Jonsson (magnus@smartelectronix.com)
243	>	uint32 used = n;
244	>	used |= used >> 1;
245	>	used |= used >> 2;
246	>	used |= used >> 4;
247	>	used |= used >> 8;
248	>	used |= used >> 16;
249	>
250	>	// Draw numbers until one is found in [0,n]
251	>	uint32 i;
252	>	do
253	>	i = randInt() & used;  // toss unused bits to shorten search
254	>	while( i > n );
255	>	return i;
256	>	}
257	>
258	>
259	>	inline void MTRand::seed( const uint32 oneSeed )
260	>	{
261	>	// Seed the generator with a simple uint32
262	>	initialize(oneSeed);
263	>	reload();
264	>	}
265	>
266	>
267	>	inline void MTRand::seed( uint32 *const bigSeed, const uint32 seedLength )
268	>	{
269	>	// Seed the generator with an array of uint32's
270	>	// There are 2^19937-1 possible initial states.  This function allows
271	>	// all of those to be accessed by providing at least 19937 bits (with a
272	>	// default seed length of N = 624 uint32's).  Any bits above the lower 32
273	>	// in each element are discarded.
274	>	// Just call seed() if you want to get array from /dev/urandom
275	>	initialize(19650218UL);
276	>	register int i = 1;
277	>	register uint32 j = 0;
278	>	register int k = ( N > seedLength ? N : seedLength );
279	>	for( ; k; --k )
280	>	{
281	>	state[i] =
282	>	state[i] ^ ( (state[i-1] ^ (state[i-1] >> 30)) * 1664525UL );
283	>	state[i] += ( bigSeed[j] & 0xffffffffUL ) + j;
284	>	state[i] &= 0xffffffffUL;
285	>	++i;  ++j;
286	>	if( i >= N ) { state[0] = state[N-1];  i = 1; }
287	>	if( j >= seedLength ) j = 0;
288	>	}
289	>	for( k = N - 1; k; --k )
290	>	{
291	>	state[i] =
292	>	state[i] ^ ( (state[i-1] ^ (state[i-1] >> 30)) * 1566083941UL );
293	>	state[i] -= i;
294	>	state[i] &= 0xffffffffUL;
295	>	++i;
296	>	if( i >= N ) { state[0] = state[N-1];  i = 1; }
297	>	}
298	>	state[0] = 0x80000000UL;  // MSB is 1, assuring non-zero initial array
299	>	reload();
300	>	}
301	>
302	>
303	>	inline void MTRand::seed()
304	>	{
305	>	std::vector<uint32> seeds;
306	>
307	>	seeds = generateSeeds();
308	>
309	>	if (seeds.size() == 1) {
310	>	seed( seeds[0] );
311	>	} else {
312	>	seed( &seeds[0], seeds.size() );
313		}
314	+	}
315	+
316	+
317	+	inline std::vector<MTRand::uint32> MTRand::generateSeeds() {
318	+	// Seed the generator with an array from /dev/urandom if available
319	+	// Otherwise use a hash of time() and clock() values
320	+
321	+	std::vector<uint32> bigSeed;
322	+
323	+	// First try getting an array from /dev/urandom
324	+	FILE* urandom = fopen( "/dev/urandom", "rb" );
325	+	if( urandom )
326	+	{
327	+	bigSeed.resize(N);
328	+	register uint32 *s = &bigSeed[0];
329	+	register int i = N;
330	+	register bool success = true;
331	+	while( success && i-- )
332	+	success = (fread( s++, sizeof(uint32), 1, urandom ) == 0);
333	+	fclose(urandom);
334	+	if( success ) { return bigSeed; }
335	+	}
336
337	<	// Was not successful, so use time() and clock() instead
337	>	// Was not successful, so use time() and clock() instead
338
339	<	bigSeed.push_back(hash( time(NULL), clock()));
340	<	return bigSeed;
341	<	}
339	>	bigSeed.push_back(hash( time(NULL), clock()));
340	>	return bigSeed;
341	>	}
342
343
344	<	inline void MTRand::initialize( const uint32 seed )
345	<	{
346	<	// Initialize generator state with seed
347	<	// See Knuth TAOCP Vol 2, 3rd Ed, p.106 for multiplier.
348	<	// In previous versions, most significant bits (MSBs) of the seed affect
349	<	// only MSBs of the state array.  Modified 9 Jan 2002 by Makoto Matsumoto.
350	<	register uint32 *s = state;
351	<	register uint32 *r = state;
352	<	register int i = 1;
353	<	*s++ = seed & 0xffffffffUL;
354	<	for( ; i < N; ++i )
355	<	{
356	<	*s++ = ( 1812433253UL * ( *r ^ (*r >> 30) ) + i ) & 0xffffffffUL;
357	<	r++;
358	<	}
359	<	}
344	>	inline void MTRand::initialize( const uint32 seed )
345	>	{
346	>	// Initialize generator state with seed
347	>	// See Knuth TAOCP Vol 2, 3rd Ed, p.106 for multiplier.
348	>	// In previous versions, most significant bits (MSBs) of the seed affect
349	>	// only MSBs of the state array.  Modified 9 Jan 2002 by Makoto Matsumoto.
350	>	register uint32 *s = state;
351	>	register uint32 *r = state;
352	>	register int i = 1;
353	>	*s++ = seed & 0xffffffffUL;
354	>	for( ; i < N; ++i )
355	>	{
356	>	*s++ = ( 1812433253UL * ( *r ^ (*r >> 30) ) + i ) & 0xffffffffUL;
357	>	r++;
358	>	}
359	>	}
360
361
362	<	inline void MTRand::reload()
363	<	{
364	<	// Generate N new values in state
365	<	// Made clearer and faster by Matthew Bellew (matthew.bellew@home.com)
366	<	register uint32 *p = state;
367	<	register int i;
368	<	for( i = N - M; i--; ++p )
369	<	*p = twist( p[M], p[0], p[1] );
370	<	for( i = M; --i; ++p )
371	<	*p = twist( p[M-N], p[0], p[1] );
372	<	*p = twist( p[M-N], p[0], state[0] );
362	>	inline void MTRand::reload()
363	>	{
364	>	// Generate N new values in state
365	>	// Made clearer and faster by Matthew Bellew (matthew.bellew@home.com)
366	>	register uint32 *p = state;
367	>	register int i;
368	>	for( i = N - M; i--; ++p )
369	>	*p = twist( p[M], p[0], p[1] );
370	>	for( i = M; --i; ++p )
371	>	*p = twist( p[M-N], p[0], p[1] );
372	>	*p = twist( p[M-N], p[0], state[0] );
373
374	<	left = N, pNext = state;
375	<	}
374	>	left = N, pNext = state;
375	>	}
376
377
378	<	inline MTRand::uint32 MTRand::hash( time_t t, clock_t c )
379	<	{
380	<	// Get a uint32 from t and c
381	<	// Better than uint32(x) in case x is floating point in [0,1]
382	<	// Based on code by Lawrence Kirby (fred@genesis.demon.co.uk)
378	>	inline MTRand::uint32 MTRand::hash( time_t t, clock_t c )
379	>	{
380	>	// Get a uint32 from t and c
381	>	// Better than uint32(x) in case x is floating point in [0,1]
382	>	// Based on code by Lawrence Kirby (fred@genesis.demon.co.uk)
383
384	<	static uint32 differ = 0;  // guarantee time-based seeds will change
384	>	static uint32 differ = 0;  // guarantee time-based seeds will change
385
386	<	uint32 h1 = 0;
387	<	unsigned char *p = (unsigned char *) &t;
388	<	for( size_t i = 0; i < sizeof(t); ++i )
389	<	{
390	<	h1 *= UCHAR_MAX + 2U;
391	<	h1 += p[i];
392	<	}
393	<	uint32 h2 = 0;
394	<	p = (unsigned char *) &c;
395	<	for( size_t j = 0; j < sizeof(c); ++j )
396	<	{
397	<	h2 *= UCHAR_MAX + 2U;
398	<	h2 += p[j];
399	<	}
400	<	return ( h1 + differ++ ) ^ h2;
401	<	}
386	>	uint32 h1 = 0;
387	>	unsigned char *p = (unsigned char *) &t;
388	>	for( size_t i = 0; i < sizeof(t); ++i )
389	>	{
390	>	h1 *= UCHAR_MAX + 2U;
391	>	h1 += p[i];
392	>	}
393	>	uint32 h2 = 0;
394	>	p = (unsigned char *) &c;
395	>	for( size_t j = 0; j < sizeof(c); ++j )
396	>	{
397	>	h2 *= UCHAR_MAX + 2U;
398	>	h2 += p[j];
399	>	}
400	>	return ( h1 + differ++ ) ^ h2;
401	>	}
402
403
404	<	inline void MTRand::save( uint32* saveArray ) const
405	<	{
406	<	register uint32 *sa = saveArray;
407	<	register const uint32 *s = state;
408	<	register int i = N;
409	<	for( ; i--; *sa++ = *s++ ) {}
410	<	*sa = left;
411	<	}
404	>	inline void MTRand::save( uint32* saveArray ) const
405	>	{
406	>	register uint32 *sa = saveArray;
407	>	register const uint32 *s = state;
408	>	register int i = N;
409	>	for( ; i--; *sa++ = *s++ ) {}
410	>	*sa = left;
411	>	}
412
413
414	<	inline void MTRand::load( uint32 *const loadArray )
415	<	{
416	<	register uint32 *s = state;
417	<	register uint32 *la = loadArray;
418	<	register int i = N;
419	<	for( ; i--; *s++ = *la++ ) {}
420	<	left = *la;
421	<	pNext = &state[N-left];
422	<	}
414	>	inline void MTRand::load( uint32 *const loadArray )
415	>	{
416	>	register uint32 *s = state;
417	>	register uint32 *la = loadArray;
418	>	register int i = N;
419	>	for( ; i--; *s++ = *la++ ) {}
420	>	left = *la;
421	>	pNext = &state[N-left];
422	>	}
423
424
425	<	inline std::ostream& operator<<( std::ostream& os, const MTRand& mtrand )
426	<	{
427	<	register const MTRand::uint32 *s = mtrand.state;
428	<	register int i = mtrand.N;
429	<	for( ; i--; os << *s++ << "\t" ) {}
430	<	return os << mtrand.left;
431	<	}
425	>	inline std::ostream& operator<<( std::ostream& os, const MTRand& mtrand )
426	>	{
427	>	register const MTRand::uint32 *s = mtrand.state;
428	>	register int i = mtrand.N;
429	>	for( ; i--; os << *s++ << "\t" ) {}
430	>	return os << mtrand.left;
431	>	}
432
433
434	<	inline std::istream& operator>>( std::istream& is, MTRand& mtrand )
435	<	{
436	<	register MTRand::uint32 *s = mtrand.state;
437	<	register int i = mtrand.N;
438	<	for( ; i--; is >> *s++ ) {}
439	<	is >> mtrand.left;
440	<	mtrand.pNext = &mtrand.state[mtrand.N-mtrand.left];
441	<	return is;
442	<	}
434	>	inline std::istream& operator>>( std::istream& is, MTRand& mtrand )
435	>	{
436	>	register MTRand::uint32 *s = mtrand.state;
437	>	register int i = mtrand.N;
438	>	for( ; i--; is >> *s++ ) {}
439	>	is >> mtrand.left;
440	>	mtrand.pNext = &mtrand.state[mtrand.N-mtrand.left];
441	>	return is;
442	>	}
443
444		}
445		#endif  // MERSENNETWISTER_H

Comparing:
trunk/src/math/MersenneTwister.hpp (property svn:keywords), Revision 385 by tim, Tue Mar 1 20:10:14 2005 UTC vs.
branches/development/src/math/MersenneTwister.hpp (property svn:keywords), Revision 1767 by gezelter, Fri Jul 6 22:01:58 2012 UTC

#	Line 0 | Line 1
1	+	Author Id Revision Date

Diff Legend

–	Removed lines
+	Added lines
<	Changed lines
>	Changed lines