MersenneTwister64.hpp

Go to the documentation of this file.

#ifndef MTTOOLBOX_MERSENNETWISTER64_HPP
#define MTTOOLBOX_MERSENNETWISTER64_HPP

#include <stdint.h>
#include <inttypes.h>
#include <string>
#include <MTToolBox/AbstractGenerator.hpp>
#if defined(DEBUG)
#include <iostream>
#include <iomanip>
#endif

namespace MTToolBox {
    class MersenneTwister64 : public AbstractGenerator<uint64_t> {
    public:
        MersenneTwister64() {
            mt = new uint64_t[N];
            seed(UINT64_C(19650218));
        }

        MersenneTwister64(uint64_t value) {
            mt = new uint64_t[N];
            seed(value);
        }

        MersenneTwister64(const std::string& value) {
            mt = new uint64_t[N];
            seed(value);
        }

        MersenneTwister64(const uint64_t *value, int size) {
            mt = new uint64_t[N];
            seed(value, size);
        }

        ~MersenneTwister64() {
            delete[] mt;
        }

        void seed(uint64_t value) {
            mt[0] = value;
            for (mti = 1; mti < N; mti++) {
                mt[mti] = mti
                    + UINT64_C(6364136223846793005)
                    * (mt[mti-1] ^ (mt[mti-1] >> 62));
            }
            mti = 0;
        }

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

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

        template<class T> void seed_array(const T *init_key,
                                          int key_length) {
            int i, j;
            uint64_t k;
            seed(UINT64_C(19650218));
            i = 1;
            j = 0;
            k = (N > key_length ? N : key_length);
            for (; k; k--) {
                mt[i] = (mt[i] ^ ((mt[i-1] ^ (mt[i-1] >> 62))
                                  * UINT64_C(3935559000370003845)))
                    + init_key[j] + (uint64_t)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] >> 62))
                                  * UINT64_C(2862933555777941757)))
                    - (uint64_t)i; /* non linear */
                i++;
                if (i>=N) {
                    mt[0] = mt[N-1];
                    i=1;
                }
            }
            /* MSB is 1; assuring non-zero initial array */
            mt[0] = UINT64_C(1) << 63;
            mti = 0;
        }

        uint64_t generate() {
            return next();
        }

        uint64_t next() {
            const uint64_t UPPER_MASK = UINT64_C(0xFFFFFFFF80000000);
            const uint64_t LOWER_MASK = UINT64_C(0x000000007FFFFFFF);
            const uint64_t mag01[2] = {0x0UL, UINT64_C(0xB5026F5AA96619E9)};
            uint64_t x;

            x = (mt[mti] & UPPER_MASK) | (mt[(mti + 1) % N] & LOWER_MASK);
            mt[mti] = mt[(mti + M) % N] ^ (x >> 1)
                ^ mag01[(int)(x & UINT64_C(1))];
            x = mt[mti];
            mti = (mti + 1) % N;
            return temper(x);
        }

        int bitSize() const {
            return 19937;
        }
    private:
        enum {N = 312, M = 156};
        uint64_t *mt;    /* the array for the state vector  */
        int mti;
        uint64_t temper(uint64_t y) {
            y ^= (y >> 29) & UINT64_C(0x5555555555555555);
            y ^= (y << 17) & UINT64_C(0x71D67FFFEDA60000);
            y ^= (y << 37) & UINT64_C(0xFFF7EEE000000000);
            y ^= (y >> 43);
            return y;
        }
    };
}
#endif // MTTOOLBOX_MERSENNETWISTER64_HPP