MersenneTwister.hpp

Go to the documentation of this file.

#ifndef MTTOOLBOX_MERSENNETWISTER_HPP
#define MTTOOLBOX_MERSENNETWISTER_HPP

#include <stdint.h>
#include <inttypes.h>
#include <string>
#include <MTToolBox/AbstractGenerator.hpp>

namespace MTToolBox {
    class MersenneTwister : public AbstractGenerator<uint32_t> {
    public:
        MersenneTwister() {
            mt = new uint32_t[LARGE_N];
            seed(5489);
        }

        MersenneTwister(uint32_t value) {
            mt = new uint32_t[LARGE_N];
            seed(value);
        }

        MersenneTwister(const std::string& value) {
            mt = new uint32_t[LARGE_N];
            seed(value);
        }

        MersenneTwister(const uint32_t *value, int size) {
            mt = new uint32_t[LARGE_N];
            seed(value, size);
        }

        ~MersenneTwister() {
            delete[] mt;
        }

        void seed(uint32_t value) {
            mt[0] = value;
            for (mti = 1; mti < N; mti++) {
                mt[mti] =
                    (UINT32_C(1812433253)
                     * (mt[mti - 1] ^ (mt[mti - 1] >> 30)) + mti);
            }
        }

        void seed(const std::string& value) {
            seed_array<char>(value.c_str(), static_cast<int>(value.size()));
        }

        void seed(const uint32_t *value, int key_length) {
            seed_array<uint32_t>(value, key_length);
        }

        template<class T> void seed_array(const T *value, int key_length) {
            int i, j, k;
            seed(UINT32_C(19650218));
            i = 1;
            j = 0;
            if (N > key_length) {
                k = N;
            } else {
                k = key_length;
            }
            for (; k > 0; k--) {
                mt[i] = (mt[i] ^ ((mt[i-1] ^ (mt[i-1] >> 30))
                                  * UINT32_C(1664525)))
                    + static_cast<uint32_t>(value[j]) + j; /* non linear */
                i++;
                j++;
                if (i>=N) {
                    mt[0] = mt[N-1];
                    i=1;
                }
                if (j>=key_length) {
                    j=0;
                }
            }
            for (k=N-1; k; k--) {
                mt[i] = (mt[i] ^ ((mt[i-1] ^ (mt[i-1] >> 30))
                                  * UINT32_C(1566083941)))
                    - i; /* non linear */
                i++;
                if (i>=N) {
                    mt[0] = mt[N-1];
                    i=1;
                }
            }
            /* MSB is 1; assuring non-zero initial array */
            mt[0] = UINT32_C(0x80000000);
            mti = N;
        }

        uint32_t generate() {
            return next();
        }

        uint32_t next() {
            using namespace std;

            const uint32_t LARGE_MASK = UINT32_C(0x3ff);
            const uint32_t UPPER_MASK = UINT32_C(0x80000000);
            const uint32_t LOWER_MASK = UINT32_C(0x7fffffff);
            const uint32_t mag01[2] = {0x0UL, UINT32_C(0x9908b0df)};
            uint32_t y;

            y = (mt[(LARGE_N - N + mti) & LARGE_MASK] & UPPER_MASK)
                | (mt[(LARGE_N - N + mti + 1) & LARGE_MASK] & LOWER_MASK);
            mt[mti] = mt[(LARGE_N - N + M + mti) & LARGE_MASK]
                ^ (y >> 1) ^ mag01[y & 0x1UL];
            y = temper(mt[mti]);
            mti = (mti + 1) & LARGE_MASK;
            return y;
        }

        int bitSize() const {
            return 19937;
        }
    private:
        enum {LARGE_N = 1024, N = 624, M = 397};
        uint32_t *mt;    /* the array for the state vector  */
        int mti;
        uint32_t temper(uint32_t y) {
            y ^= (y >> 11);
            y ^= (y << 7) & UINT32_C(0x9d2c5680);
            y ^= (y << 15) & UINT32_C(0xefc60000);
            y ^= (y >> 18);
            return y;
        }
    };
}
#endif // MTTOOLBOX_MERSENNETWISTER_HPP