Added Mersenne Twister random number generator.

3 files changed, 78 insertions, 1 deletions

diff --git a/src/stable/randommersenne.cpp b/src/stable/randommersenne.cpp
new file mode 100644
index 0000000..2d9e85b
--- /dev/null
+++ b/src/stable/randommersenne.cpp
@@ -0,0 +1,44 @@
+#include "bu/randommersenne.h"
+
+Bu::RandomMersenne::RandomMersenne( int32_t iSeed ) :
+    iIndex( 0 )
+{
+    seed( iSeed );
+}
+
+Bu::RandomMersenne::~RandomMersenne()
+{
+}
+
+void Bu::RandomMersenne::seed( int32_t iSeed )
+{
+     iMt[0] = iSeed;
+     for( int j = 1; j < 624; j++ )
+     {
+         iMt[j] = 0x6c078965 * (iMt[j-1]^(iMt[j-1]>>30)) + j;
+     }
+}
+
+int32_t Bu::RandomMersenne::rand()
+{
+    if( iIndex == 0 )
+    {
+        for( int j = 0; j < 624; j++ )
+        {
+            uint32_t y = (iMt[j] & 0x80000000) + (iMt[(j+1)%624] & 0x7fffffff);
+            iMt[j] = iMt[(j+397)%624] ^ (y>>1);
+            if( (y%2) )
+                iMt[j] = iMt[j]^0x9908b0df;
+        }
+    }
+
+    uint32_t y = iMt[iIndex];
+    y = y ^ (y>>11);
+    y = y ^ ((y<<7)&0x9d2c5680);
+    y = y ^ ((y<<15)&0xefc60000);
+    y = y ^ (y>>18);
+
+    iIndex = (iIndex + 1) % 624;
+    return y;
+}
+
diff --git a/src/stable/randommersenne.h b/src/stable/randommersenne.h
new file mode 100644
index 0000000..27c100e
--- /dev/null
+++ b/src/stable/randommersenne.h
@@ -0,0 +1,31 @@
+#ifndef BU_RANDOM_MERSENNE_H
+#define BU_RANDOM_MERSENNE_H
+
+#include "bu/randombase.h"
+
+namespace Bu
+{
+    /**
+     * An implementation of Mersenne Twister (MT19937) algorithm.  This as an
+     * algorithm with an excellent reputation for very good random number
+     * generation and a very large period.  It is however, relatively slow and
+     * complex compared to, say the Complementary Multiply With Carry
+     * (Bu::RandomCmwc), and may not actually be a *better* random number
+     * generator.
+     */
+    class RandomMersenne : public Bu::RandomBase
+    {
+    public:
+        RandomMersenne( int32_t iSeed=0 );
+        virtual ~RandomMersenne();
+
+        virtual void seed( int32_t iSeed );
+        virtual int32_t rand();
+
+    private:
+        int32_t iIndex;
+        uint32_t iMt[624];
+    };
+};
+
+#endif
diff --git a/src/tests/random.cpp b/src/tests/random.cpp
index 95cdb0a..d7b94ce 100644
--- a/src/tests/random.cpp
+++ b/src/tests/random.cpp
@@ -1,6 +1,7 @@
 #include <bu/randombasic.h>
 #include <bu/randomcmwc.h>
 #include <bu/randomsystem.h>
+#include <bu/randommersenne.h>
 #include <bu/sio.h>
 #include <time.h>
 
@@ -15,7 +16,7 @@ void coverage()
     uint32_t uBucket[78];
     memset( uBucket, 0, sizeof(uint32_t)*78 );
 
-    for( int j = 0; j < 1000; j++ )
+    for( int j = 0; j < 1000000; j++ )
     {
         uBucket[(int)(((uint32_t)rand.rand())/(double)(0xfffffffful)*78+0.5)]++;
     }
@@ -43,6 +44,7 @@ int main()
     coverage<RandomBasic>();
     coverage<RandomCmwc>();
     coverage<RandomSystem>();
+    coverage<RandomMersenne>();
 
     return 0;
 }