Commit e119c87b authored by Stephen M. Webb's avatar Stephen M. Webb

added extra rng test case and fix to mt

parent 85026cce
......@@ -36,11 +36,26 @@ Author: Stephen M. Webb <stephen.webb@bregmasoft.ca>
//
-// 1. Redistributions of source code must retain the above copyright
-// notice, this list of conditions and the following disclaimer.
-//
+// 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.,
+// 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
//
-// 2. Redistributions in binary form must reproduce the above copyright
-// notice, this list of conditions and the following disclaimer in the
-// documentation and/or other materials provided with the distribution.
-//
+// For information about WorldForge and its authors, please contact
+// the Worldforge Web Site at http://www.worldforge.org.
//
-// 3. The names of its contributors may not be used to endorse or promote
-// products derived from this software without specific prior written
-// permission.
......@@ -73,23 +88,6 @@ Author: Stephen M. Webb <stephen.webb@bregmasoft.ca>
-
-// Not thread safe (unless auto-initialization is avoided and each thread has
-// its own MTRand object)
+// 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.,
+// 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+//
+// For information about WorldForge and its authors, please contact
+// the Worldforge Web Site at http://www.worldforge.org.
+//
+#ifndef WFMATH_MERSENNE_TWISTER_H_
+#define WFMATH_MERSENNE_TWISTER_H_
......@@ -248,10 +246,11 @@ Author: Stephen M. Webb <stephen.webb@bregmasoft.ca>
-
-inline double MTRand::randExc()
- { return double(randInt()) * (1.0/4294967296.0); }
-
+{ return rand() * n; }
-inline double MTRand::randExc( const double& n )
- { return randExc() * n; }
-
-inline double MTRand::randDblExc()
- { return ( double(randInt()) + 0.5 ) * (1.0/4294967296.0); }
-
......@@ -259,11 +258,18 @@ Author: Stephen M. Webb <stephen.webb@bregmasoft.ca>
- { return randDblExc() * n; }
-
-inline double MTRand::rand53()
-{
+inline MTRand::uint32 MTRand::randInt(uint32 n)
{
- uint32 a = randInt() >> 5, b = randInt() >> 6;
- return ( a * 67108864.0 + b ) * (1.0/9007199254740992.0); // by Isaku Wada
-}
-
+ uint32 used = n;
+ used |= used >> 1;
+ used |= used >> 2;
+ used |= used >> 4;
+ used |= used >> 8;
+ used |= used >> 16;
-inline double MTRand::randNorm( const double& mean, const double& variance )
-{
- // Return a real number from a normal (Gaussian) distribution with given
......@@ -271,8 +277,13 @@ Author: Stephen M. Webb <stephen.webb@bregmasoft.ca>
- double r = std::sqrt( -2.0 * std::log( 1.0-randDblExc()) ) * variance;
- double phi = 2.0 * 3.14159265358979323846264338328 * randExc();
- return mean + r * std::cos(phi);
-}
-
+ uint32 i;
+ do
+ i = randInt() & used;
+ while( i > n );
+ return i;
}
-inline MTRand::uint32 MTRand::randInt()
-{
- // Pull a 32-bit integer from the generator state
......@@ -288,9 +299,11 @@ Author: Stephen M. Webb <stephen.webb@bregmasoft.ca>
- s1 ^= (s1 << 15) & 0xefc60000UL;
- return ( s1 ^ (s1 >> 18) );
-}
-
-inline MTRand::uint32 MTRand::randInt( const uint32& n )
-{
+#if 0
+inline void MTRand::save(uint32* saveArray) const
{
- // Find which bits are used in n
- // Optimized by Magnus Jonsson (magnus@smartelectronix.com)
- uint32 used = n;
......@@ -306,21 +319,25 @@ Author: Stephen M. Webb <stephen.webb@bregmasoft.ca>
- i = randInt() & used; // toss unused bits to shorten search
- while( i > n );
- return i;
-}
-
-
+ register uint32 *sa = saveArray;
+ register const uint32 *s = state;
+ register int i = state_size;
+ for( ; i--; *sa++ = *s++ ) {}
+ *sa = left;
}
-inline void MTRand::seed( const uint32 oneSeed )
-{
+inline void MTRand::load(uint32 *const loadArray)
{
- // Seed the generator with a simple uint32
- initialize(oneSeed);
- reload();
-}
+{ return rand() * n; }
-
-
-inline void MTRand::seed( uint32 *const bigSeed, const uint32 seedLength )
+inline MTRand::uint32 MTRand::randInt(uint32 n)
{
-{
- // Seed the generator with an array of uint32's
- // There are 2^19937-1 possible initial states. This function allows
- // all of those to be accessed by providing at least 19937 bits (with a
......@@ -352,25 +369,11 @@ Author: Stephen M. Webb <stephen.webb@bregmasoft.ca>
- }
- state[0] = 0x80000000UL; // MSB is 1, assuring non-zero initial array
- reload();
+ uint32 used = n;
+ used |= used >> 1;
+ used |= used >> 2;
+ used |= used >> 4;
+ used |= used >> 8;
+ used |= used >> 16;
+
+ uint32 i;
+ do
+ i = randInt() & used;
+ while( i > n );
+ return i;
}
-}
-
-
-inline void MTRand::initialize( const uint32 seed )
+#if 0
+inline void MTRand::save(uint32* saveArray) const
{
-{
- // Initialize generator state with seed
- // See Knuth TAOCP Vol 2, 3rd Ed, p.106 for multiplier.
- // In previous versions, most significant bits (MSBs) of the seed affect
......@@ -384,17 +387,11 @@ Author: Stephen M. Webb <stephen.webb@bregmasoft.ca>
- *s++ = ( 1812433253UL * ( *r ^ (*r >> 30) ) + i ) & 0xffffffffUL;
- r++;
- }
+ register uint32 *sa = saveArray;
+ register const uint32 *s = state;
+ register int i = state_size;
+ for( ; i--; *sa++ = *s++ ) {}
+ *sa = left;
}
-}
-
-
-inline void MTRand::reload()
+inline void MTRand::load(uint32 *const loadArray)
{
-{
- // Generate N new values in state
- // Made clearer and faster by Matthew Bellew (matthew.bellew@home.com)
- register uint32 *p = state;
......@@ -660,9 +657,9 @@ Author: Stephen M. Webb <stephen.webb@bregmasoft.ca>
+void MTRand::seed(uint32 s)
+{
+ state[0] = s;
+ for (uint32 i = 1; i < state_size; ++i)
+ for (index = 1; index < state_size; ++index)
+ {
+ state[i] = (1812433253UL * (state[i-1] ^ (state[i-1] >> 30)) + i);
+ state[index] = (1812433253UL * (state[index-1] ^ (state[index-1] >> 30)) + index);
+ }
+}
+
......@@ -700,7 +697,7 @@ Author: Stephen M. Webb <stephen.webb@bregmasoft.ca>
+ /* generate state_size words at one time */
+ if (index >= state_size)
+ {
+ int kk;
+ uint32 kk;
+ for (kk=0; kk < state_size - period; kk++)
+ {
+ y = (state[kk]&UPPER_MASK) | (state[kk+1]&LOWER_MASK);
......@@ -708,10 +705,10 @@ Author: Stephen M. Webb <stephen.webb@bregmasoft.ca>
+ }
+ for (; kk < state_size-1; kk++)
+ {
+ y = (state[kk]&UPPER_MASK)|(state[kk+1]&LOWER_MASK);
+ y = (state[kk]&UPPER_MASK) | (state[kk+1]&LOWER_MASK);
+ state[kk] = state[kk+(period - state_size)] ^ (y >> 1) ^ mag01[y & 0x01];
+ }
+ y = (state[state_size-1]&UPPER_MASK)|(state[0]&LOWER_MASK);
+ y = (state[state_size-1]&UPPER_MASK) | (state[0]&LOWER_MASK);
+ state[state_size-1] = state[period-1] ^ (y >> 1) ^ mag01[y & 0x1UL];
+
+ index = 0;
......@@ -750,7 +747,60 @@ Author: Stephen M. Webb <stephen.webb@bregmasoft.ca>
}
--- a/wfmath/randgen_test.cpp
+++ b/wfmath/randgen_test.cpp
@@ -41,13 +41,8 @@
@@ -23,8 +23,10 @@
// Author: Alistair Riddoch
// Created: 2011-2-14
+#include <assert.h>
#include "randgen.h"
#include <cstdio>
+#include <iostream>
#ifdef HAVE_CONFIG_H
#include "config.h"
@@ -32,7 +34,40 @@
using WFMath::MTRand;
-int main()
+bool test_known_sequence()
+{
+ static WFMath::MTRand::uint32 expected_results[] = {
+ 2221777491u,
+ 2873750246u,
+ 4067173416u,
+ 794519497u,
+ 3287624630u,
+ 3357287912u,
+ 1212880927u,
+ 2464917741u,
+ 949382604u,
+ 1898004827u
+ };
+ WFMath::MTRand rng;
+ bool result = true;
+
+ rng.seed(23);
+
+ for (int i = 0; i < 10; ++i)
+ {
+ WFMath::MTRand::uint32 rnd = rng.randInt();
+ if (rnd != expected_results[i])
+ {
+ std::cerr << "Mismatch between MTRand and known result sequuence\n"
+ << rnd << " != " << expected_results[i] << std::endl;
+ result = false;
+ }
+ assert(rnd == expected_results[i]);
+ }
+ return result;
+}
+
+bool test_generator_instances()
{
MTRand one(23);
@@ -41,13 +76,14 @@
MTRand two(23);
......@@ -763,9 +813,14 @@ Author: Stephen M. Webb <stephen.webb@bregmasoft.ca>
-
- printf("%.16f %.16f\n", thr.rand<float>(), thr.rand<float>());
- float thrres = thr.rand<float>();
-
- printf("%.16f %.16f %.16f\n", oneres, twores, thrres);
+ printf("%.16f %.16f\n", oneres, twores);
+ return oneres == twores;
+}
- printf("%.16f %.16f %.16f\n", oneres, twores, thrres);
+int main()
+{
+ return !(test_known_sequence() && test_generator_instances());
}
--- a/wfmath/stream.cpp
+++ b/wfmath/stream.cpp
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment