MT Rand

// Minimalistic Mersenne Twister Engine
// Author(s):       iFarbod <[email protected]>
//
// Copyright (c) 2014-2016 San Andreas Online
// Use of this source code is governed by the MIT license that can be
// found in the LICENSE file.

#pragma once
#include <cstdint> // for int types
#include <cstddef> // for size_t

#define PREVENT_MACRO_SUBSTITUTION

template<class UIntType,
    size_t w, size_t n, size_t m, size_t r,
    UIntType a, size_t u, UIntType d, size_t s,
    UIntType b, size_t t,
    UIntType c, size_t l, UIntType f>
class MersenneTwister
{
public:
    typedef UIntType result_type;

    static constexpr size_t word_size = w;
    static constexpr size_t state_size = n;
    static constexpr size_t shift_size = m;
    static constexpr size_t mask_bits = r;
    static constexpr UIntType xor_mask = a;
    static constexpr size_t tempering_u = u;
    static constexpr UIntType tempering_d = d;
    static constexpr size_t tempering_s = s;
    static constexpr UIntType tempering_b = b;
    static constexpr size_t tempering_t = t;
    static constexpr UIntType tempering_c = c;
    static constexpr size_t tempering_l = l;
    static constexpr UIntType initialization_multiplier = f;
    static constexpr UIntType default_seed = 5489u;

    // backwards compatibility
    static constexpr UIntType parameter_a = a;
    static constexpr size_t output_u = u;
    static constexpr size_t output_s = s;
    static constexpr UIntType output_b = b;
    static constexpr size_t output_t = t;
    static constexpr UIntType output_c = c;
    static constexpr size_t output_l = l;

    // old concept requirements
    static constexpr bool has_fixed_range = false;

    MersenneTwister(UIntType _Seed = default_seed)
    {
        seed(_Seed);
    }

    /**
    * Sets the state x(0) to v mod 2w. Then, iteratively,
    * sets x(i) to
    * (i + f * (x(i-1) xor (x(i-1) rshift w-2))) mod 2<sup>w</sup>
    * for i = 1 .. n-1. x(n) is the first value to be returned by operator().
    */
    void seed(const UIntType& value)
    {
        // New seeding algorithm from
        // http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/MT2002/emt19937ar.html
        // In the previous versions, MSBs of the seed affected only MSBs of the
        // state x[].
        const UIntType mask = (max)();
        x[0] = value & mask;
        for (i = 1; i < n; i++) {
            // See Knuth "The Art of Computer Programming"
            // Vol. 2, 3rd ed., page 106
            x[i] = (f * (x[i - 1] ^ (x[i - 1] >> (w - 2))) + i) & mask;
        }

        normalize_state();
    }

    void seed() { seed(default_seed); }

    /** Returns the smallest value that the generator can produce. */
    static result_type min PREVENT_MACRO_SUBSTITUTION()
    {
        return 0;
    }
    /** Returns the largest value that the generator can produce. */
    static result_type max PREVENT_MACRO_SUBSTITUTION()
    {
        return _WMSK;
    }

    result_type operator()()
    {
        if (i == n)
            twist();
        // Step 4
        UIntType z = x[i];
        ++i;
        z ^= ((z >> u) & d);
        z ^= ((z << s) & b);
        z ^= ((z << t) & c);
        z ^= (z >> l);
        return z;
    }

    void discard(uintmax_t z)
    {
        for (uintmax_t j = 0; j < z; ++j) {
            (*this)();
        }
    }
private:
    void twist()
    {
        const UIntType upper_mask = (~static_cast<UIntType>(0)) << r;
        const UIntType lower_mask = ~upper_mask;

        const size_t unroll_factor = 6;
        const size_t unroll_extra1 = (n - m) % unroll_factor;
        const size_t unroll_extra2 = (m - 1) % unroll_factor;

        // split loop to avoid costly modulo operations
        {  // extra scope for MSVC brokenness w.r.t. for scope
            for (size_t j = 0; j < n - m - unroll_extra1; j++) {
                UIntType y = (x[j] & upper_mask) | (x[j + 1] & lower_mask);
                x[j] = x[j + m] ^ (y >> 1) ^ ((x[j + 1] & 1) * a);
            }
        }
        {
            for (size_t j = n - m - unroll_extra1; j < n - m; j++) {
                UIntType y = (x[j] & upper_mask) | (x[j + 1] & lower_mask);
                x[j] = x[j + m] ^ (y >> 1) ^ ((x[j + 1] & 1) * a);
            }
        }
        {
            for (size_t j = n - m; j < n - 1 - unroll_extra2; j++) {
                UIntType y = (x[j] & upper_mask) | (x[j + 1] & lower_mask);
                x[j] = x[j - (n - m)] ^ (y >> 1) ^ ((x[j + 1] & 1) * a);
            }
        }
        {
            for (size_t j = n - 1 - unroll_extra2; j < n - 1; j++) {
                UIntType y = (x[j] & upper_mask) | (x[j + 1] & lower_mask);
                x[j] = x[j - (n - m)] ^ (y >> 1) ^ ((x[j + 1] & 1) * a);
            }
        }
        // last iteration
        UIntType y = (x[n - 1] & upper_mask) | (x[0] & lower_mask);
        x[n - 1] = x[m - 1] ^ (y >> 1) ^ ((x[0] & 1) * a);
        i = 0;
    }

    void normalize_state()
    {
        const UIntType upper_mask = (~static_cast<UIntType>(0)) << r;
        const UIntType lower_mask = ~upper_mask;
        UIntType y0 = x[m - 1] ^ x[n - 1];
        if (y0 & (static_cast<UIntType>(1) << (w - 1))) {
            y0 = ((y0 ^ a) << 1) | 1;
        }
        else {
            y0 = y0 << 1;
        }
        x[0] = (x[0] & upper_mask) | (y0 & lower_mask);

        // fix up the state if it's all zeroes.
        for (size_t j = 0; j < n; ++j) {
            if (x[j] != 0) return;
        }
        x[0] = static_cast<UIntType>(1) << (w - 1);
    }

    UIntType x[n];
    size_t i;

    static constexpr UIntType _WMSK = ~((~UIntType(0) << (w - 1)) << 1);
    static constexpr UIntType _HMSK = (_WMSK << r) & _WMSK;
    static constexpr UIntType _LMSK = ~_HMSK & _WMSK;
};

// Memory usage:
// mt11213b   : 352 * sizeof(uint32_t) = 352 * 4 = 1408 bytes = ~1 KB
// mt19937    : 625 * sizeof(uint32_t) = 625 * 4 = 2500 bytes = ~2 KB
// mt19937_64 : 312 * sizeof(uint64_t) = 312 * 8 = 2496 bytes = ~2 KB

typedef MersenneTwister<uint32_t, 32, 351, 175, 19, 0xccab8ee7,
    11, 0xffffffff, 7, 0x31b6ab00, 15, 0xffe50000, 17, 1812433253> MT11213b;

typedef MersenneTwister<uint32_t, 32, 624, 397, 31, 0x9908b0df,
    11, 0xffffffff, 7, 0x9d2c5680, 15, 0xefc60000, 18, 1812433253> MT19937;

typedef MersenneTwister<uint64_t, 64, 312, 156, 31,
    UINT64_C(0xb5026f5aa96619e9), 29, UINT64_C(0x5555555555555555), 17,
    UINT64_C(0x71d67fffeda60000), 37, UINT64_C(0xfff7eee000000000), 43,
    UINT64_C(6364136223846793005)> MT19937_64;