1 /* boost random/mersenne_twister.hpp header file 2 * 3 * Copyright Jens Maurer 2000-2001 4 * Distributed under the Boost Software License, Version 1.0. (See 5 * accompanying file LICENSE_1_0.txt or copy at 6 * http://www.boost.org/LICENSE_1_0.txt) 7 * 8 * See http://www.boost.org for most recent version including documentation. 9 * 10 * $Id: mersenne_twister.hpp,v 1.20 2005/07/21 22:04:31 jmaurer Exp $ 11 * 12 * Revision history 13 * 2001-02-18 moved to individual header files 14 */ 15 16 #ifndef DLIB_BOOST_RANDOM_MERSENNE_TWISTER_HPP 17 #define DLIB_BOOST_RANDOM_MERSENNE_TWISTER_HPP 18 19 #include <iostream> 20 #include <algorithm> // std::copy 21 #include <stdexcept> 22 #include "../uintn.h" 23 #include "../serialize.h" 24 25 namespace dlib 26 { 27 namespace random_helpers 28 { 29 30 // ------------------------------------------------------------------------------------ 31 32 // http://www.math.keio.ac.jp/matumoto/emt.html 33 template< 34 class UIntType, 35 int w, 36 int n, 37 int m, 38 int r, 39 UIntType a, 40 int u, 41 int s, 42 UIntType b, 43 int t, 44 UIntType c, 45 int l, 46 UIntType val 47 > 48 class mersenne_twister 49 { 50 public: 51 typedef UIntType result_type; 52 const static int word_size = w; 53 const static int state_size = n; 54 const static int shift_size = m; 55 const static int mask_bits = r; 56 const static UIntType parameter_a = a; 57 const static int output_u = u; 58 const static int output_s = s; 59 const static UIntType output_b = b; 60 const static int output_t = t; 61 const static UIntType output_c = c; 62 const static int output_l = l; 63 64 const static bool has_fixed_range = false; 65 mersenne_twister()66 mersenne_twister() { seed(); } 67 mersenne_twister(UIntType value)68 explicit mersenne_twister(UIntType value) { seed(value); } 69 seed()70 void seed () { seed(UIntType(5489)); } 71 72 // compiler-generated copy ctor and assignment operator are fine 73 seed(UIntType value)74 void seed(UIntType value) 75 { 76 // New seeding algorithm from 77 // http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/MT2002/emt19937ar.html 78 // In the previous versions, MSBs of the seed affected only MSBs of the 79 // state x[]. 80 const UIntType mask = ~0u; 81 x[0] = value & mask; 82 for (i = 1; i < n; i++) { 83 // See Knuth "The Art of Computer Programming" Vol. 2, 3rd ed., page 106 84 x[i] = (1812433253UL * (x[i-1] ^ (x[i-1] >> (w-2))) + i) & mask; 85 } 86 } 87 min()88 result_type min() const { return 0; } max()89 result_type max() const 90 { 91 // avoid "left shift count >= with of type" warning 92 result_type res = 0; 93 for(int i = 0; i < w; ++i) 94 res |= (1u << i); 95 return res; 96 } 97 98 result_type operator()(); 99 serialize(const mersenne_twister & item,std::ostream & out)100 friend void serialize( 101 const mersenne_twister& item, 102 std::ostream& out 103 ) 104 { 105 dlib::serialize(item.x, out); 106 dlib::serialize(item.i, out); 107 } 108 deserialize(mersenne_twister & item,std::istream & in)109 friend void deserialize( 110 mersenne_twister& item, 111 std::istream& in 112 ) 113 { 114 dlib::deserialize(item.x, in); 115 dlib::deserialize(item.i, in); 116 } 117 118 private: 119 120 void twist(int block); 121 122 // state representation: next output is o(x(i)) 123 // x[0] ... x[k] x[k+1] ... x[n-1] x[n] ... x[2*n-1] represents 124 // x(i-k) ... x(i) x(i+1) ... x(i-k+n-1) x(i-k-n) ... x[i(i-k-1)] 125 // The goal is to always have x(i-n) ... x(i-1) available for 126 // operator== and save/restore. 127 128 UIntType x[2*n]; 129 int i; 130 }; 131 132 // ------------------------------------------------------------------------------------ 133 134 template< 135 class UIntType, int w, int n, int m, int r, UIntType a, int u, 136 int s, UIntType b, int t, UIntType c, int l, UIntType val 137 > twist(int block)138 void mersenne_twister<UIntType,w,n,m,r,a,u,s,b,t,c,l,val>::twist( 139 int block 140 ) 141 { 142 const UIntType upper_mask = (~0u) << r; 143 const UIntType lower_mask = ~upper_mask; 144 145 if(block == 0) { 146 for(int j = n; j < 2*n; j++) { 147 UIntType y = (x[j-n] & upper_mask) | (x[j-(n-1)] & lower_mask); 148 x[j] = x[j-(n-m)] ^ (y >> 1) ^ (y&1 ? a : 0); 149 } 150 } else if (block == 1) { 151 // split loop to avoid costly modulo operations 152 { // extra scope for MSVC brokenness w.r.t. for scope 153 for(int j = 0; j < n-m; j++) { 154 UIntType y = (x[j+n] & upper_mask) | (x[j+n+1] & lower_mask); 155 x[j] = x[j+n+m] ^ (y >> 1) ^ (y&1 ? a : 0); 156 } 157 } 158 159 for(int j = n-m; j < n-1; j++) { 160 UIntType y = (x[j+n] & upper_mask) | (x[j+n+1] & lower_mask); 161 x[j] = x[j-(n-m)] ^ (y >> 1) ^ (y&1 ? a : 0); 162 } 163 // last iteration 164 UIntType y = (x[2*n-1] & upper_mask) | (x[0] & lower_mask); 165 x[n-1] = x[m-1] ^ (y >> 1) ^ (y&1 ? a : 0); 166 i = 0; 167 } 168 } 169 170 // ------------------------------------------------------------------------------------ 171 172 template< 173 class UIntType, int w, int n, int m, int r, UIntType a, int u, 174 int s, UIntType b, int t, UIntType c, int l, UIntType val 175 > 176 inline typename mersenne_twister<UIntType,w,n,m,r,a,u,s,b,t,c,l,val>::result_type operator()177 mersenne_twister<UIntType,w,n,m,r,a,u,s,b,t,c,l,val>::operator()( 178 ) 179 { 180 if(i == n) 181 twist(0); 182 else if(i >= 2*n) 183 twist(1); 184 // Step 4 185 UIntType z = x[i]; 186 ++i; 187 z ^= (z >> u); 188 z ^= ((z << s) & b); 189 z ^= ((z << t) & c); 190 z ^= (z >> l); 191 return z; 192 } 193 194 // ------------------------------------------------------------------------------------ 195 196 } // namespace random 197 198 199 typedef random_helpers::mersenne_twister<uint32,32,351,175,19,0xccab8ee7,11, 200 7,0x31b6ab00,15,0xffe50000,17, 0xa37d3c92> mt11213b; 201 202 // validation by experiment from mt19937.c 203 typedef random_helpers::mersenne_twister<uint32,32,624,397,31,0x9908b0df,11, 204 7,0x9d2c5680,15,0xefc60000,18, 3346425566U> mt19937; 205 206 } // namespace dlib 207 208 209 #endif // DLIB_BOOST_RANDOM_MERSENNE_TWISTER_HPP 210 211