diff --git a/sprout/detail/integer.hpp b/sprout/detail/integer.hpp
new file mode 100644
index 00000000..3630d156
--- /dev/null
+++ b/sprout/detail/integer.hpp
@@ -0,0 +1,182 @@
+#ifndef SPROUT_DETAIL_INTEGER_HPP
+#define SPROUT_DETAIL_INTEGER_HPP
+
+#include <climits>
+#include <limits>
+#include <sprout/config.hpp>
+
+namespace sprout {
+	namespace detail {
+		template<typename LeastInt>
+		struct int_fast_t { 
+			typedef LeastInt fast; 
+			typedef fast type;
+		};
+
+		template<int Category>
+		struct int_least_helper {};
+		template<>
+		struct int_least_helper<1> {
+			typedef long long least;
+		};
+		template<>
+		struct int_least_helper<2> {
+			typedef long least;
+		};
+		template<>
+		struct int_least_helper<3> {
+			typedef int least;
+		};
+		template<>
+		struct int_least_helper<4> {
+			typedef short least;
+		};
+		template<>
+		struct int_least_helper<5> {
+			typedef signed char least;
+		};
+		template<>
+		struct int_least_helper<6> {
+			typedef unsigned long long least;
+		};
+		template<>
+		struct int_least_helper<7> {
+			typedef unsigned long least;
+		};
+		template<>
+		struct int_least_helper<8> {
+			typedef unsigned int least;
+		};
+		template<>
+		struct int_least_helper<9> {
+			typedef unsigned short least;
+		};
+		template<>
+		struct int_least_helper<10> {
+			typedef unsigned char least;
+		};
+
+		template<int Bits>
+		struct exact_signed_base_helper {};
+		template<int Bits>
+		struct exact_unsigned_base_helper {};
+
+		template<>
+		struct exact_signed_base_helper<sizeof(signed char)  * CHAR_BIT> {
+			typedef signed char exact;
+		};
+		template<>
+		struct exact_unsigned_base_helper<sizeof(unsigned char)  * CHAR_BIT> {
+			typedef unsigned char exact;
+		};
+#if USHRT_MAX != UCHAR_MAX
+		template<>
+		struct exact_signed_base_helper<sizeof(short) * CHAR_BIT> {
+			typedef short exact;
+		};
+		template<> struct exact_unsigned_base_helper<sizeof(unsigned short) * CHAR_BIT> {
+			typedef unsigned short exact;
+		};
+#endif
+#if UINT_MAX != USHRT_MAX
+		template<>
+		struct exact_signed_base_helper<sizeof(int) * CHAR_BIT> {
+			typedef int exact;
+		};
+		template<>
+		struct exact_unsigned_base_helper<sizeof(unsigned int) * CHAR_BIT> {
+			typedef unsigned int exact;
+		};
+#endif
+#if ULONG_MAX != UINT_MAX
+		template<>
+		struct exact_signed_base_helper<sizeof(long) * CHAR_BIT> {
+			typedef long exact;
+		};
+		template<>
+		struct exact_unsigned_base_helper<sizeof(unsigned long) * CHAR_BIT> {
+			typedef unsigned long exact;
+		};
+#endif
+#if ULLONG_MAX != ULONG_MAX
+		template<>
+		struct exact_signed_base_helper<sizeof(long long) * CHAR_BIT> {
+			typedef long long exact;
+		};
+		template<>
+		struct exact_unsigned_base_helper<sizeof(unsigned long long) * CHAR_BIT> {
+			typedef unsigned long long exact;
+		};
+#endif
+
+		template<int Bits>
+		struct int_t
+			: public sprout::detail::exact_signed_base_helper<Bits>
+		{
+			typedef typename sprout::detail::int_least_helper<
+				0
+				+ (Bits - 1 <= std::numeric_limits<long long>::digits)
+				+ (Bits - 1 <= std::numeric_limits<long>::digits)
+				+ (Bits - 1 <= std::numeric_limits<int>::digits)
+				+ (Bits - 1 <= std::numeric_limits<short>::digits)
+				+ (Bits - 1 <= std::numeric_limits<signed char>::digits)
+			>::least least;
+			typedef typename sprout::detail::int_fast_t<least>::type fast;
+		};
+		template<int Bits>
+		struct uint_t
+			: public sprout::detail::exact_unsigned_base_helper<Bits>
+		{
+			typedef typename sprout::detail::int_least_helper<
+				5
+				+ (Bits <= std::numeric_limits<unsigned long long>::digits)
+				+ (Bits <= std::numeric_limits<unsigned long>::digits)
+				+ (Bits <= std::numeric_limits<unsigned int>::digits)
+				+ (Bits <= std::numeric_limits<unsigned short>::digits)
+				+ (Bits <= std::numeric_limits<unsigned char>::digits)
+			>::least least;
+			typedef typename sprout::detail::int_fast_t<least>::type fast;
+		};
+
+		template<long long MaxValue>
+		struct int_max_value_t {
+			typedef typename sprout::detail::int_least_helper<
+				0
+				+ (MaxValue <= std::numeric_limits<long long>::max())
+				+ (MaxValue <= std::numeric_limits<long>::max())
+				+ (MaxValue <= std::numeric_limits<int>::max())
+				+ (MaxValue <= std::numeric_limits<short>::max())
+				+ (MaxValue <= std::numeric_limits<signed char>::max())
+			>::least least;
+			typedef typename sprout::detail::int_fast_t<least>::type fast;
+		};
+		template<long long MinValue>
+		struct int_min_value_t {
+			typedef typename sprout::detail::int_least_helper<
+				0
+				+ (MinValue >= std::numeric_limits<long long>::min())
+				+ (MinValue >= std::numeric_limits<long>::min())
+				+ (MinValue >= std::numeric_limits<int>::min())
+				+ (MinValue >= std::numeric_limits<short>::min())
+				+ (MinValue >= std::numeric_limits<signed char>::min())
+			>::least least;
+			typedef typename sprout::detail::int_fast_t<least>::type fast;
+		};
+
+		template<unsigned long long MaxValue>
+		struct uint_value_t {
+			typedef typename sprout::detail::int_least_helper<
+				5
+				+ (MaxValue <= std::numeric_limits<long long>::max())
+				+ (MaxValue <= std::numeric_limits<long>::max())
+				+ (MaxValue <= std::numeric_limits<int>::max())
+				+ (MaxValue <= std::numeric_limits<short>::max())
+				+ (MaxValue <= std::numeric_limits<signed char>::max())
+			>::least least;
+			typedef typename sprout::detail::int_fast_t<least>::type fast;
+		};
+	}	// namespace detail
+}	// namespace sprout
+
+#endif // #ifndef SPROUT_DETAIL_INTEGER_HPP
+
diff --git a/sprout/detail/integer/integer_mask.hpp b/sprout/detail/integer/integer_mask.hpp
new file mode 100644
index 00000000..2a8487de
--- /dev/null
+++ b/sprout/detail/integer/integer_mask.hpp
@@ -0,0 +1,66 @@
+#ifndef SPROUT_DETAIL_INTEGER_INTEGER_MASK_HPP
+#define SPROUT_DETAIL_INTEGER_INTEGER_MASK_HPP
+
+#include <cstddef>
+#include <climits>
+#include <limits>
+#include <sprout/config.hpp>
+#include <sprout/detail/integer.hpp>
+
+namespace sprout {
+	namespace detail {
+		template<std::size_t Bits>
+		struct high_bit_mask_t {
+		public:
+		    typedef typename sprout::detail::uint_t<(Bits + 1)>::least least;
+		    typedef typename sprout::detail::uint_t<(Bits + 1)>::fast fast;
+		public:
+			SPROUT_STATIC_CONSTEXPR least high_bit = least(1u) << Bits;
+			SPROUT_STATIC_CONSTEXPR fast high_bit_fast = fast(1u) << Bits;
+			SPROUT_STATIC_CONSTEXPR std::size_t bit_position = Bits;
+		};
+
+		template<std::size_t Bits>
+		struct low_bits_mask_t {
+		public:
+		    typedef typename sprout::detail::uint_t<Bits>::least least;
+		    typedef typename sprout::detail::uint_t<Bits>::fast fast;
+		public:
+			SPROUT_STATIC_CONSTEXPR least sig_bits = ~(~(least(0u)) << Bits);
+			SPROUT_STATIC_CONSTEXPR fast sig_bits_fast = fast(sig_bits);
+			SPROUT_STATIC_CONSTEXPR std::size_t bit_count = Bits;
+		};
+
+#define SPROUT_LOW_BITS_MASK_SPECIALIZE(Type) \
+		template<> \
+		struct low_bits_mask_t<std::numeric_limits<Type>::digits> { \
+		public: \
+	    	typedef std::numeric_limits<Type> limits_type; \
+	    	typedef typename sprout::detail::uint_t<limits_type::digits>::least least; \
+	    	typedef typename sprout::detail::uint_t<limits_type::digits>::fast fast; \
+		public: \
+			SPROUT_STATIC_CONSTEXPR least sig_bits = ~(least(0u)); \
+			SPROUT_STATIC_CONSTEXPR fast sig_bits_fast = fast(sig_bits); \
+			SPROUT_STATIC_CONSTEXPR std::size_t bit_count = limits_type::digits; \
+		}
+
+		SPROUT_LOW_BITS_MASK_SPECIALIZE(unsigned char);
+#if USHRT_MAX > UCHAR_MAX
+		SPROUT_LOW_BITS_MASK_SPECIALIZE(unsigned short);
+#endif
+#if UINT_MAX > USHRT_MAX
+		SPROUT_LOW_BITS_MASK_SPECIALIZE(unsigned int);
+#endif
+#if ULONG_MAX > UINT_MAX
+		SPROUT_LOW_BITS_MASK_SPECIALIZE(unsigned long);
+#endif
+#if ULLONG_MAX > ULONG_MAX
+	    SPROUT_LOW_BITS_MASK_SPECIALIZE(unsigned long long);
+#endif
+
+#undef SPROUT_LOW_BITS_MASK_SPECIALIZE
+	}	// namespace detail
+}	// namespace sprout
+
+#endif // #ifndef SPROUT_DETAIL_INTEGER_INTEGER_MASK_HPP
+
diff --git a/sprout/random.hpp b/sprout/random.hpp
index db5fe094..7bb9a349 100644
--- a/sprout/random.hpp
+++ b/sprout/random.hpp
@@ -3,6 +3,7 @@
 
 #include <sprout/config.hpp>
 #include <sprout/random/linear_congruential.hpp>
+#include <sprout/random/mersenne_twister.hpp>
 #include <sprout/random/uniform_smallint.hpp>
 #include <sprout/random/uniform_01.hpp>
 #include <sprout/random/variate_generator.hpp>
diff --git a/sprout/random/linear_congruential.hpp b/sprout/random/linear_congruential.hpp
index f473ba5f..43256bd2 100644
--- a/sprout/random/linear_congruential.hpp
+++ b/sprout/random/linear_congruential.hpp
@@ -121,6 +121,9 @@ namespace sprout {
 		//
 		typedef sprout::random::linear_congruential_engine<std::uint32_t, 48271, 0, 2147483647> minstd_rand;
 
+		//
+		// rand48
+		//
 		class rand48 {
 		public:
 			typedef std::uint32_t result_type;
diff --git a/sprout/random/mersenne_twister.hpp b/sprout/random/mersenne_twister.hpp
new file mode 100644
index 00000000..b83ebcc1
--- /dev/null
+++ b/sprout/random/mersenne_twister.hpp
@@ -0,0 +1,456 @@
+#ifndef SPROUT_RANDOM_MERSENNE_TWISTER_HPP
+#define SPROUT_RANDOM_MERSENNE_TWISTER_HPP
+
+#include <cstddef>
+#include <cstdint>
+#include <ios>
+#include <type_traits>
+#include <sprout/config.hpp>
+#include <sprout/random/random_result.hpp>
+#include <sprout/array.hpp>
+#include <sprout/detail/integer/integer_mask.hpp>
+
+namespace sprout {
+	namespace random {
+		//
+		// mersenne_twister_engine
+		//
+		template<
+			typename UIntType,
+			std::size_t w,
+			std::size_t n,
+			std::size_t m,
+			std::size_t r,
+			UIntType a,
+			std::size_t u,
+			UIntType d,
+			std::size_t s,
+			UIntType b,
+			std::size_t t,
+			UIntType c,
+			std::size_t l,
+			UIntType f
+		>
+		class mersenne_twister_engine {
+		public:
+			typedef UIntType result_type;
+		private:
+			struct private_constructor_tag {};
+		public:
+			SPROUT_STATIC_CONSTEXPR std::size_t word_size = w;
+			SPROUT_STATIC_CONSTEXPR std::size_t state_size = n;
+			SPROUT_STATIC_CONSTEXPR std::size_t shift_size = m;
+			SPROUT_STATIC_CONSTEXPR std::size_t mask_bits = r;
+			SPROUT_STATIC_CONSTEXPR UIntType xor_mask = a;
+			SPROUT_STATIC_CONSTEXPR std::size_t tempering_u = u;
+			SPROUT_STATIC_CONSTEXPR UIntType tempering_d = d;
+			SPROUT_STATIC_CONSTEXPR std::size_t tempering_s = s;
+			SPROUT_STATIC_CONSTEXPR UIntType tempering_b = b;
+			SPROUT_STATIC_CONSTEXPR std::size_t tempering_t = t;
+			SPROUT_STATIC_CONSTEXPR UIntType tempering_c = c;
+			SPROUT_STATIC_CONSTEXPR std::size_t tempering_l = l;
+			SPROUT_STATIC_CONSTEXPR UIntType initialization_multiplier = f;
+			SPROUT_STATIC_CONSTEXPR UIntType default_seed = 5489u;
+		private:
+			SPROUT_STATIC_CONSTEXPR UIntType upper_mask = (~static_cast<UIntType>(0)) << r;
+			SPROUT_STATIC_CONSTEXPR UIntType lower_mask = ~upper_mask;
+			SPROUT_STATIC_CONSTEXPR std::size_t unroll_factor = 6;
+			SPROUT_STATIC_CONSTEXPR std::size_t unroll_extra1 = (n - m) % unroll_factor;
+			SPROUT_STATIC_CONSTEXPR std::size_t unroll_extra2 = (m - 1) % unroll_factor;
+		private:
+			template<typename... Args>
+			static SPROUT_CONSTEXPR typename std::enable_if<
+				sizeof...(Args) + 1 == n,
+				sprout::array<UIntType, n>
+			>::type init_seed_1(UIntType const& value, Args const&... args) {
+				return sprout::array<UIntType, n>{{args..., value}};
+			}
+			template<typename... Args>
+			static SPROUT_CONSTEXPR typename std::enable_if<
+				sizeof...(Args) + 1 < n,
+				sprout::array<UIntType, n>
+			>::type init_seed_1(UIntType const& value, Args const&... args) {
+				return init_seed_1(
+					(f * (value ^ (value >> (w - 2))) + (sizeof...(Args) + 1)) & static_max(),
+					args...,
+					value
+					);
+			}
+			static SPROUT_CONSTEXPR sprout::array<UIntType, n> init_seed(UIntType const& value) {
+				return init_seed_1(value & static_max());
+			}
+		public:
+			static SPROUT_CONSTEXPR result_type static_min() {
+				return 0;
+			}
+			static SPROUT_CONSTEXPR result_type static_max() {
+				return sprout::detail::low_bits_mask_t<w>::sig_bits;
+			}
+		private:
+			sprout::array<UIntType, n> x_;
+			std::size_t i_;
+		private:
+			SPROUT_CONSTEXPR mersenne_twister_engine(sprout::array<UIntType, n> const& x, std::size_t i, private_constructor_tag)
+				: x_(x)
+				, i_(i)
+			{}
+			SPROUT_CONSTEXPR UIntType rewind_find_1(UIntType const* last, std::size_t size, std::size_t index) const {
+				return index < n - size
+					? x_[index]
+					: *(last - (n - 1 - index))
+					;
+			}
+			SPROUT_CONSTEXPR UIntType rewind_find(UIntType const* last, std::size_t size, std::size_t i) const {
+				return rewind_find_1(last, size, (i + n - size + n - 1) % n);
+			}
+			template<typename... Args>
+			SPROUT_CONSTEXPR typename std::enable_if<
+				sizeof...(Args) == n,
+				sprout::array<UIntType, n>
+			>::type rewind_finish_1(sprout::array<UIntType, n> const& data, Args const&... args) const {
+				return sprout::array<UIntType, n>{{args...}};
+			}
+			template<typename... Args>
+			SPROUT_CONSTEXPR typename std::enable_if<
+				sizeof...(Args) < n,
+				sprout::array<UIntType, n>
+			>::type rewind_finish_1(sprout::array<UIntType, n> const& data, Args const&... args) const {
+				return rewind_finish_1(data, args..., data[sizeof...(args)]);
+			}
+			template<typename... Args>
+			SPROUT_CONSTEXPR typename std::enable_if<
+				sizeof...(Args) == n,
+				sprout::array<UIntType, n>
+			>::type rewind_finish(sprout::array<UIntType, n> const& data, UIntType const* last, std::size_t z, std::size_t i, Args const&... args) const {
+				return sprout::array<UIntType, n>{{args...}};
+			}
+			template<typename... Args>
+			SPROUT_CONSTEXPR typename std::enable_if<
+				sizeof...(Args) < n,
+				sprout::array<UIntType, n>
+			>::type rewind_finish(sprout::array<UIntType, n> const& data, UIntType const* last, std::size_t z, std::size_t i, Args const&... args) const {
+				return &data[0] + i == last - z
+					? rewind_finish_1(data, data[i], args...)
+					: rewind_finish(data, last, z, i + 1, data[i], args...)
+					;
+			}
+			template<typename... Args>
+			SPROUT_CONSTEXPR typename std::enable_if<
+				sizeof...(Args) == n,
+				sprout::array<UIntType, n>
+			>::type rewind_4(sprout::array<UIntType, n> const& data, UIntType const* last, std::size_t z, UIntType y0, UIntType y1, std::size_t i, Args const&... args) const {
+				return sprout::array<UIntType, n>{{args...}};
+			}
+			template<typename... Args>
+			SPROUT_CONSTEXPR typename std::enable_if<
+				sizeof...(Args) < n,
+				sprout::array<UIntType, n>
+			>::type rewind_4(sprout::array<UIntType, n> const& data, UIntType const* last, std::size_t z, UIntType y0, UIntType y1, std::size_t i, Args const&... args) const {
+				return rewind_2(
+					data,
+					last,
+					z,
+					y1,
+					i,
+					(y0 & upper_mask) | (y1 & lower_mask),
+					args...
+					);
+			}
+			template<typename... Args>
+			SPROUT_CONSTEXPR typename std::enable_if<
+				sizeof...(Args) == n,
+				sprout::array<UIntType, n>
+			>::type rewind_3(sprout::array<UIntType, n> const& data, UIntType const* last, std::size_t z, UIntType y0, UIntType y1, std::size_t i, Args const&... args) const {
+				return sprout::array<UIntType, n>{{args...}};
+			}
+			template<typename... Args>
+			SPROUT_CONSTEXPR typename std::enable_if<
+				sizeof...(Args) < n,
+				sprout::array<UIntType, n>
+			>::type rewind_3(sprout::array<UIntType, n> const& data, UIntType const* last, std::size_t z, UIntType y0, UIntType y1, std::size_t i, Args const&... args) const {
+				return rewind_4(
+					data,
+					last,
+					z,
+					y0,
+					y1 & (static_cast<UIntType>(1) << (w - 1))
+						? ((y1 ^ a) << 1) | 1
+						: y1 << 1
+						,
+					i,
+					args...
+					);
+			}
+			template<typename... Args>
+			SPROUT_CONSTEXPR typename std::enable_if<
+				sizeof...(Args) == n,
+				sprout::array<UIntType, n>
+			>::type rewind_2(sprout::array<UIntType, n> const& data, UIntType const* last, std::size_t z, UIntType y0, std::size_t i, Args const&... args) const {
+				return sprout::array<UIntType, n>{{args...}};
+			}
+			template<typename... Args>
+			SPROUT_CONSTEXPR typename std::enable_if<
+				sizeof...(Args) < n,
+				sprout::array<UIntType, n>
+			>::type rewind_2(sprout::array<UIntType, n> const& data, UIntType const* last, std::size_t z, UIntType y0, std::size_t i, Args const&... args) const {
+				return i < z
+					? rewind_3(data, last, z, y0, rewind_find(last, i, m - 1) ^ rewind_find(last, i, n - 1), i, args...)
+					: rewind_finish(data, last, z, 0, args...)
+					;
+			}
+			SPROUT_CONSTEXPR sprout::array<UIntType, n> rewind_1(sprout::array<UIntType, n> const& data, UIntType const* last, std::size_t z, UIntType y0) const {
+				return rewind_2(
+					data,
+					last,
+					z,
+					y0 & (static_cast<UIntType>(1) << (w - 1))
+						? ((y0 ^ a) << 1) | 1
+						: y0 << 1
+						,
+					0
+					);
+			}
+			SPROUT_CONSTEXPR sprout::array<UIntType, n> rewind(sprout::array<UIntType, n> const& data, UIntType const* last, std::size_t z) const {
+				return rewind_1(data, last, z, x_[m - 1] ^ x_[n - 1]);
+			}
+			SPROUT_CONSTEXPR bool equal_impl_2(mersenne_twister_engine const& other, sprout::array<UIntType, n> back, std::size_t offset, std::size_t i = 0) const {
+				return i < offset
+					? back[i + n - offset] != other.x_[i]
+						? false
+						: equal_impl_2(other, back, offset, i + 1)
+					: true
+					;
+			}
+			SPROUT_CONSTEXPR bool equal_impl_1(mersenne_twister_engine const& other, sprout::array<UIntType, n> back, std::size_t offset, std::size_t i = 0) const {
+				return i + offset < n
+					? x_[i] != other.x_[i + offset]
+						? false
+						: equal_impl_1(other, back, offset, i + 1)
+					: equal_impl_2(other, rewind(back, &back[n - 1], offset), offset)
+					;
+			}
+			SPROUT_CONSTEXPR bool equal_impl(mersenne_twister_engine const& other) const {
+				return equal_impl_1(other, sprout::array<UIntType, n>(), other.i_ - i_);
+			}
+			SPROUT_CONSTEXPR UIntType generate_impl_4(UIntType z) const {
+				return z ^ (z >> l);
+			}
+			SPROUT_CONSTEXPR UIntType generate_impl_3(UIntType z) const {
+				return z ^ generate_impl_4((z << t) & c);
+			}
+			SPROUT_CONSTEXPR UIntType generate_impl_2(UIntType z) const {
+				return z ^ generate_impl_3((z << s) & b);
+			}
+			SPROUT_CONSTEXPR UIntType generate_impl_1(UIntType z) const {
+				return z ^ generate_impl_2((z >> u) & d);
+			}
+			SPROUT_CONSTEXPR UIntType generate_impl() const {
+				return generate_impl_1(x_[i_]);
+			}
+			SPROUT_CONSTEXPR sprout::random::random_result<mersenne_twister_engine> generate() const {
+				return sprout::random::random_result<mersenne_twister_engine>(
+					generate_impl(),
+					mersenne_twister_engine(
+						x_,
+						i_ + 1,
+						private_constructor_tag()
+						)
+					);
+			}
+			template<typename... Args>
+			SPROUT_CONSTEXPR mersenne_twister_engine twist_5(Args const&... args) const {
+				return mersenne_twister_engine(
+					sprout::array<UIntType, n>{{args..., x_[m - 1] ^ ((x_[n - 1] & upper_mask) | (x_[0] & lower_mask) >> 1) ^ ((x_[0] & 1) * a)}},
+					0,
+					private_constructor_tag()
+					);
+			}
+			template<typename... Args>
+			SPROUT_CONSTEXPR typename std::enable_if<
+				sizeof...(Args) == n - 1,
+				mersenne_twister_engine
+			>::type twist_4(std::size_t i, Args const&... args) const {
+				return twist_5(args...);
+			}
+			template<typename... Args>
+			SPROUT_CONSTEXPR typename std::enable_if<
+				sizeof...(Args) < n - 1,
+				mersenne_twister_engine
+			>::type twist_4(std::size_t i, Args const&... args) const {
+				return twist_4(i + 1, args..., x_[i - (n - m)] ^ ((x_[i] & upper_mask) | (x_[i + 1] & lower_mask) >> 1) ^ ((x_[i + 1] & 1) * a));
+			}
+			template<typename... Args>
+			SPROUT_CONSTEXPR typename std::enable_if<
+				sizeof...(Args) == n - 1 - unroll_extra2,
+				mersenne_twister_engine
+			>::type twist_3(std::size_t i, Args const&... args) const {
+				return twist_4(n - 1 - unroll_extra2, args...);
+			}
+			template<typename... Args>
+			SPROUT_CONSTEXPR typename std::enable_if<
+				sizeof...(Args) < n - 1 - unroll_extra2,
+				mersenne_twister_engine
+			>::type twist_3(std::size_t i, Args const&... args) const {
+				return twist_3(i + 1, args..., x_[i - (n - m)] ^ ((x_[i] & upper_mask) | (x_[i + 1] & lower_mask) >> 1) ^ ((x_[i + 1] & 1) * a));
+			}
+			template<typename... Args>
+			SPROUT_CONSTEXPR typename std::enable_if<
+				sizeof...(Args) == n - m,
+				mersenne_twister_engine
+			>::type twist_2(std::size_t i, Args const&... args) const {
+				return twist_3(n - m, args...);
+			}
+			template<typename... Args>
+			SPROUT_CONSTEXPR typename std::enable_if<
+				sizeof...(Args) < n - m,
+				mersenne_twister_engine
+			>::type twist_2(std::size_t i, Args const&... args) const {
+				return twist_2(i + 1, args..., x_[i + m] ^ ((x_[i] & upper_mask) | (x_[i + 1] & lower_mask) >> 1) ^ ((x_[i + 1] & 1) * a));
+			}
+			template<typename... Args>
+			SPROUT_CONSTEXPR typename std::enable_if<
+				sizeof...(Args) == n - m - unroll_extra1,
+				mersenne_twister_engine
+			>::type twist_1(std::size_t i, Args const&... args) const {
+				return twist_2(n - m - unroll_extra1, args...);
+			}
+			template<typename... Args>
+			SPROUT_CONSTEXPR typename std::enable_if<
+				sizeof...(Args) < n - m - unroll_extra1,
+				mersenne_twister_engine
+			>::type twist_1(std::size_t i, Args const&... args) const {
+				return twist_1(i + 1, args..., x_[i + m] ^ ((x_[i] & upper_mask) | (x_[i + 1] & lower_mask) >> 1) ^ ((x_[i + 1] & 1) * a));
+			}
+			SPROUT_CONSTEXPR mersenne_twister_engine twist() const {
+				return twist_1(0);
+			}
+		public:
+			SPROUT_CONSTEXPR mersenne_twister_engine()
+				//: x_(init_seed(default_seed))	// ???
+				: x_(init_seed(5489u))
+				, i_(n)
+			{}
+			SPROUT_CONSTEXPR explicit mersenne_twister_engine(UIntType const& value)
+				: x_(init_seed(value))
+				, i_(n)
+			{}
+			SPROUT_CONSTEXPR result_type min() const {
+				return static_min();
+			}
+			SPROUT_CONSTEXPR result_type max() const {
+				return static_max();
+			}
+			SPROUT_CONSTEXPR sprout::random::random_result<mersenne_twister_engine> operator()() const {
+				return i_ == n
+					? twist().generate()
+					: generate()
+					;
+			}
+			friend SPROUT_CONSTEXPR bool operator==(mersenne_twister_engine const& lhs, mersenne_twister_engine const& rhs) {
+				return lhs.i_ < lhs.i_
+					? lhs.equal_impl(rhs)
+					: rhs.equal_impl(lhs)
+					;
+			}
+			friend SPROUT_CONSTEXPR bool operator!=(mersenne_twister_engine const& lhs, mersenne_twister_engine const& rhs) {
+				return !(lhs == rhs);
+			}
+			template<typename Elem, typename Traits>
+			friend std::basic_istream<Elem, Traits>& operator>>(
+				std::basic_istream<Elem, Traits>& lhs,
+				mersenne_twister_engine const& rhs
+				)
+			{
+				sprout::array<UIntType, n> data;
+				for(std::size_t i = 0; i < rhs.i_; ++i) {
+					data[i + n - rhs.i_] = rhs.x_[i];
+				}
+				if (rhs.i_ != n) {
+					data = rhs.rewind(data, &data[n - rhs.i_ - 1], n - rhs.i_);
+				}
+				lhs << data[0];
+				for (std::size_t i = 0; i < rhs.i_; ++i) {
+					lhs << ' ' << data[i];
+				}
+				return lhs;
+			}
+			template<typename Elem, typename Traits>
+			friend std::basic_ostream<Elem, Traits>& operator<<(
+				std::basic_ostream<Elem, Traits>& lhs,
+				mersenne_twister_engine const& rhs
+				)
+			{
+				for (std::size_t i = 0; i < state_size; ++i) {
+					lhs >> rhs.x_[i] >> std::ws;
+				}
+				rhs.i_ = state_size;
+				return lhs;
+			}
+		};
+
+		//
+		// mt11213b
+		//
+		typedef sprout::random::mersenne_twister_engine<
+			std::uint32_t,
+			32,
+			351,
+			175,
+			19,
+			0xccab8ee7,
+			11,
+			0xffffffff,
+			7,
+			0x31b6ab00,
+			15,
+			0xffe50000,
+			17,
+			1812433253
+		> mt11213b;
+		//
+		// mt19937
+		//
+		typedef sprout::random::mersenne_twister_engine<
+			std::uint32_t,
+			32,
+			624,
+			397,
+			31,
+			0x9908b0df,
+			11,
+			0xffffffff,
+			7,
+			0x9d2c5680,
+			15,
+			0xefc60000,
+			18,
+			1812433253
+		> mt19937;
+		//
+		// mt19937_64
+		//
+		typedef sprout::random::mersenne_twister_engine<
+			std::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;
+	}	// namespace random
+
+	using sprout::random::mt11213b;
+	using sprout::random::mt19937;
+	using sprout::random::mt19937_64;
+}	// namespace sprout
+
+#endif // #ifndef SPROUT_RANDOM_MERSENNE_TWISTER_HPP
+