Mersenne Twister Randome Number Genrator

mersenne_twister_engine.hpp

/*
 * Copyright (C) 2016 LAN Xingcan
 * All right reserved
 *
 * Permission to use, copy, modify, and/or distribute this software for any
 * purpose with or without fee is hereby granted, provided that the above
 * copyright notice and this permission notice appear in all copies.
 *
 * THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL WARRANTIES
 * WITH REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED WARRANTIES OF
 * MERCHANTABILITY AND FITNESS. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR
 * ANY SPECIAL, DIRECT, INDIRECT, OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES
 * WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN AN
 * ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF
 * OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE.
 */

#pragma once

#include <limits>

namespace lanxc
{

  template<typename T,
      std::size_t WORD_SIZE, std::size_t STATE_SIZE, std::size_t SHIFT_SIZE,
      std::size_t MASK_BIT, T XOR_MASK,
      std::size_t TEMPERING_U, T TEMPERING_D,
      std::size_t TEMPERING_S, T TEMPERING_B,
      std::size_t TEMPERING_T, T TEMPERING_C,
      std::size_t TEMPERING_L, T INIT_MULTIPLIER
  >
  class mersenne_twister_engine
  {
  public:
    mersenne_twister_engine(T seed = T(0)) throw()
        : index(0)
    {
      pool[0] = seed;
      for (std::size_t i = 1; i < STATE_SIZE; ++i)
      {
        pool[i] = INIT_MULTIPLIER * (pool[i - 1] ^ (pool[i - 1] >> (WORD_SIZE - 2))) + i;
      }
    }

    T operator ()() throw()
    {
      T j = index + 1;
      T k = index + SHIFT_SIZE;
      const T table[2] = { 0, XOR_MASK };

      if (j >= STATE_SIZE) j = 0;
      if (k >= STATE_SIZE) k -= STATE_SIZE;

      const T mask = ((~T(0)) >> MASK_BIT) << MASK_BIT;
      T y = (pool[index] & mask) | (pool[j] & ~mask);
      T r = pool[k] ^(y >> 1) ^table[y & 1];
      pool[index] = r;
      r ^= (r >> TEMPERING_U) & TEMPERING_D;
      r ^= ((r << TEMPERING_S) & TEMPERING_B);
      r ^= ((r << TEMPERING_T) & TEMPERING_C);
      r ^= (r >> TEMPERING_L);
      index = j;
      return r;
    };

  private:
    unsigned index;
    std::array<T, STATE_SIZE> pool;
  };

  using prng_mt19937 = mersenne_twister_engine<std::uint_fast32_t,
      32, 624, 397,
      31, 0x9908b0dfU,
      11, 0xffffffffU,
      7, 0x9d2c5680U,
      15, 0xefc60000U,
      18,
      1812433253U>;

  using prng_mt19937_64 = mersenne_twister_engine<std::uint_fast64_t,
      64, 312, 156,
      31, 0xb5026f5aa96619e9ULL,
      29, 0x5555555555555555ULL,
      17, 0x71d67fffeda60000ULL,
      37, 0xfff7eee000000000ULL,
      43,
      6364136223846793005ULL>;

}

Benchmark ( Time costed for 200000000 rounds)

Compiler	std::mt19937	lanxc::prng_mt19937
GNU G++ 6.0 with libstdc++ -Ofast	804286us	743639us
GNU G++ 6.0 with libstdc++ -O3	825381us	706003us
GNU G++ 6.0 with libstdc++ -O2	1069715us	565251us
Apple LLVM 7.3.0 Clang with libc++ -O3	1218068us	587607us
Apple LLVM 7.3.0 Clang with libc++ -O2	1205328us	586482us

As above shows, this implementation is pretty faster than the ones comes form libstdc++ or libc++.