2 files changed, 75 insertions, 0 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