File indexing completed on 2025-01-30 09:59:12
0001
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
0016
0017
0018 #ifndef BOOST_RANDOM_MERSENNE_TWISTER_HPP
0019 #define BOOST_RANDOM_MERSENNE_TWISTER_HPP
0020
0021 #include <iosfwd>
0022 #include <istream>
0023 #include <stdexcept>
0024 #include <boost/config.hpp>
0025 #include <boost/cstdint.hpp>
0026 #include <boost/integer/integer_mask.hpp>
0027 #include <boost/random/detail/config.hpp>
0028 #include <boost/random/detail/ptr_helper.hpp>
0029 #include <boost/random/detail/seed.hpp>
0030 #include <boost/random/detail/seed_impl.hpp>
0031 #include <boost/random/detail/generator_seed_seq.hpp>
0032 #include <boost/random/detail/polynomial.hpp>
0033
0034 #include <boost/random/detail/disable_warnings.hpp>
0035
0036 namespace boost {
0037 namespace random {
0038
0039
0040
0041
0042
0043
0044
0045
0046
0047
0048
0049
0050
0051
0052
0053
0054
0055
0056
0057
0058
0059
0060
0061
0062
0063
0064
0065
0066
0067 template<class UIntType,
0068 std::size_t w, std::size_t n, std::size_t m, std::size_t r,
0069 UIntType a, std::size_t u, UIntType d, std::size_t s,
0070 UIntType b, std::size_t t,
0071 UIntType c, std::size_t l, UIntType f>
0072 class mersenne_twister_engine
0073 {
0074 public:
0075 typedef UIntType result_type;
0076 BOOST_STATIC_CONSTANT(std::size_t, word_size = w);
0077 BOOST_STATIC_CONSTANT(std::size_t, state_size = n);
0078 BOOST_STATIC_CONSTANT(std::size_t, shift_size = m);
0079 BOOST_STATIC_CONSTANT(std::size_t, mask_bits = r);
0080 BOOST_STATIC_CONSTANT(UIntType, xor_mask = a);
0081 BOOST_STATIC_CONSTANT(std::size_t, tempering_u = u);
0082 BOOST_STATIC_CONSTANT(UIntType, tempering_d = d);
0083 BOOST_STATIC_CONSTANT(std::size_t, tempering_s = s);
0084 BOOST_STATIC_CONSTANT(UIntType, tempering_b = b);
0085 BOOST_STATIC_CONSTANT(std::size_t, tempering_t = t);
0086 BOOST_STATIC_CONSTANT(UIntType, tempering_c = c);
0087 BOOST_STATIC_CONSTANT(std::size_t, tempering_l = l);
0088 BOOST_STATIC_CONSTANT(UIntType, initialization_multiplier = f);
0089 BOOST_STATIC_CONSTANT(UIntType, default_seed = 5489u);
0090
0091
0092 BOOST_STATIC_CONSTANT(UIntType, parameter_a = a);
0093 BOOST_STATIC_CONSTANT(std::size_t, output_u = u);
0094 BOOST_STATIC_CONSTANT(std::size_t, output_s = s);
0095 BOOST_STATIC_CONSTANT(UIntType, output_b = b);
0096 BOOST_STATIC_CONSTANT(std::size_t, output_t = t);
0097 BOOST_STATIC_CONSTANT(UIntType, output_c = c);
0098 BOOST_STATIC_CONSTANT(std::size_t, output_l = l);
0099
0100
0101 BOOST_STATIC_CONSTANT(bool, has_fixed_range = false);
0102
0103
0104
0105
0106
0107 mersenne_twister_engine() { seed(); }
0108
0109
0110
0111
0112 BOOST_RANDOM_DETAIL_ARITHMETIC_CONSTRUCTOR(mersenne_twister_engine,
0113 UIntType, value)
0114 { seed(value); }
0115 template<class It> mersenne_twister_engine(It& first, It last)
0116 { seed(first,last); }
0117
0118
0119
0120
0121
0122
0123
0124
0125
0126 BOOST_RANDOM_DETAIL_SEED_SEQ_CONSTRUCTOR(mersenne_twister_engine,
0127 SeedSeq, seq)
0128 { seed(seq); }
0129
0130
0131
0132
0133 void seed() { seed(default_seed); }
0134
0135
0136
0137
0138
0139
0140
0141 BOOST_RANDOM_DETAIL_ARITHMETIC_SEED(mersenne_twister_engine, UIntType, value)
0142 {
0143
0144
0145
0146
0147 const UIntType mask = (max)();
0148 x[0] = value & mask;
0149 for (i = 1; i < n; i++) {
0150
0151
0152 x[i] = (f * (x[i-1] ^ (x[i-1] >> (w-2))) + i) & mask;
0153 }
0154
0155 normalize_state();
0156 }
0157
0158
0159
0160
0161 BOOST_RANDOM_DETAIL_SEED_SEQ_SEED(mersenne_twister_engine, SeeqSeq, seq)
0162 {
0163 detail::seed_array_int<w>(seq, x);
0164 i = n;
0165
0166 normalize_state();
0167 }
0168
0169
0170 template<class It>
0171 void seed(It& first, It last)
0172 {
0173 detail::fill_array_int<w>(first, last, x);
0174 i = n;
0175
0176 normalize_state();
0177 }
0178
0179
0180 static BOOST_CONSTEXPR result_type min BOOST_PREVENT_MACRO_SUBSTITUTION ()
0181 { return 0; }
0182
0183 static BOOST_CONSTEXPR result_type max BOOST_PREVENT_MACRO_SUBSTITUTION ()
0184 { return boost::low_bits_mask_t<w>::sig_bits; }
0185
0186
0187 result_type operator()();
0188
0189
0190 template<class Iter>
0191 void generate(Iter first, Iter last)
0192 { detail::generate_from_int(*this, first, last); }
0193
0194
0195
0196
0197
0198
0199
0200
0201
0202
0203 void discard(boost::uintmax_t z)
0204 {
0205 #ifndef BOOST_RANDOM_MERSENNE_TWISTER_DISCARD_THRESHOLD
0206 #define BOOST_RANDOM_MERSENNE_TWISTER_DISCARD_THRESHOLD 10000000
0207 #endif
0208 if(z > BOOST_RANDOM_MERSENNE_TWISTER_DISCARD_THRESHOLD) {
0209 discard_many(z);
0210 } else {
0211 for(boost::uintmax_t j = 0; j < z; ++j) {
0212 (*this)();
0213 }
0214 }
0215 }
0216
0217 #ifndef BOOST_RANDOM_NO_STREAM_OPERATORS
0218
0219 template<class CharT, class Traits>
0220 friend std::basic_ostream<CharT,Traits>&
0221 operator<<(std::basic_ostream<CharT,Traits>& os,
0222 const mersenne_twister_engine& mt)
0223 {
0224 mt.print(os);
0225 return os;
0226 }
0227
0228
0229 template<class CharT, class Traits>
0230 friend std::basic_istream<CharT,Traits>&
0231 operator>>(std::basic_istream<CharT,Traits>& is,
0232 mersenne_twister_engine& mt)
0233 {
0234 for(std::size_t j = 0; j < mt.state_size; ++j)
0235 is >> mt.x[j] >> std::ws;
0236
0237
0238
0239 mt.i = mt.state_size;
0240 return is;
0241 }
0242 #endif
0243
0244
0245
0246
0247
0248 friend bool operator==(const mersenne_twister_engine& x_,
0249 const mersenne_twister_engine& y_)
0250 {
0251 if(x_.i < y_.i) return x_.equal_imp(y_);
0252 else return y_.equal_imp(x_);
0253 }
0254
0255
0256
0257
0258 friend bool operator!=(const mersenne_twister_engine& x_,
0259 const mersenne_twister_engine& y_)
0260 { return !(x_ == y_); }
0261
0262 private:
0263
0264
0265 void twist();
0266
0267
0268
0269
0270
0271
0272
0273
0274
0275 bool equal_imp(const mersenne_twister_engine& other) const
0276 {
0277 UIntType back[n];
0278 std::size_t offset = other.i - i;
0279 for(std::size_t j = 0; j + offset < n; ++j)
0280 if(x[j] != other.x[j+offset])
0281 return false;
0282 rewind(&back[n-1], offset);
0283 for(std::size_t j = 0; j < offset; ++j)
0284 if(back[j + n - offset] != other.x[j])
0285 return false;
0286 return true;
0287 }
0288
0289
0290
0291
0292
0293 template<class CharT, class Traits>
0294 void print(std::basic_ostream<CharT, Traits>& os) const
0295 {
0296 UIntType data[n];
0297 for(std::size_t j = 0; j < i; ++j) {
0298 data[j + n - i] = x[j];
0299 }
0300 if(i != n) {
0301 rewind(&data[n - i - 1], n - i);
0302 }
0303 os << data[0];
0304 for(std::size_t j = 1; j < n; ++j) {
0305 os << ' ' << data[j];
0306 }
0307 }
0308
0309
0310
0311
0312
0313 void rewind(UIntType* last, std::size_t z) const
0314 {
0315 const UIntType upper_mask = (~static_cast<UIntType>(0)) << r;
0316 const UIntType lower_mask = ~upper_mask;
0317 UIntType y0 = x[m-1] ^ x[n-1];
0318 if(y0 & (static_cast<UIntType>(1) << (w-1))) {
0319 y0 = ((y0 ^ a) << 1) | 1;
0320 } else {
0321 y0 = y0 << 1;
0322 }
0323 for(std::size_t sz = 0; sz < z; ++sz) {
0324 UIntType y1 =
0325 rewind_find(last, sz, m-1) ^ rewind_find(last, sz, n-1);
0326 if(y1 & (static_cast<UIntType>(1) << (w-1))) {
0327 y1 = ((y1 ^ a) << 1) | 1;
0328 } else {
0329 y1 = y1 << 1;
0330 }
0331 *(last - sz) = (y0 & upper_mask) | (y1 & lower_mask);
0332 y0 = y1;
0333 }
0334 }
0335
0336
0337
0338
0339
0340
0341
0342
0343
0344
0345
0346 void normalize_state()
0347 {
0348 const UIntType upper_mask = (~static_cast<UIntType>(0)) << r;
0349 const UIntType lower_mask = ~upper_mask;
0350 UIntType y0 = x[m-1] ^ x[n-1];
0351 if(y0 & (static_cast<UIntType>(1) << (w-1))) {
0352 y0 = ((y0 ^ a) << 1) | 1;
0353 } else {
0354 y0 = y0 << 1;
0355 }
0356 x[0] = (x[0] & upper_mask) | (y0 & lower_mask);
0357
0358
0359 for(std::size_t j = 0; j < n; ++j) {
0360 if(x[j] != 0) return;
0361 }
0362 x[0] = static_cast<UIntType>(1) << (w-1);
0363 }
0364
0365
0366
0367
0368
0369
0370 UIntType
0371 rewind_find(UIntType* last, std::size_t size, std::size_t j) const
0372 {
0373 std::size_t index = (j + n - size + n - 1) % n;
0374 if(index < n - size) {
0375 return x[index];
0376 } else {
0377 return *(last - (n - 1 - index));
0378 }
0379 }
0380
0381
0382
0383
0384
0385
0386
0387
0388
0389
0390 void discard_many(boost::uintmax_t z)
0391 {
0392
0393
0394
0395
0396
0397
0398
0399
0400
0401
0402 detail::polynomial phi = get_characteristic_polynomial();
0403
0404
0405 detail::polynomial g = mod_pow_x(z, phi);
0406
0407
0408 detail::polynomial h;
0409 const std::size_t num_bits = w*n - r;
0410 for(std::size_t j = 0; j < num_bits * 2; ++j) {
0411
0412
0413
0414
0415 if(i >= n) twist();
0416 h[2*num_bits - j - 1] = x[i++] & UIntType(1);
0417 }
0418
0419 detail::polynomial gh = g * h;
0420 detail::polynomial result;
0421 for(std::size_t j = 0; j <= num_bits; ++j) {
0422 result[j] = gh[2*num_bits - j - 1];
0423 }
0424 reconstruct_state(result);
0425 }
0426 static detail::polynomial get_characteristic_polynomial()
0427 {
0428 const std::size_t num_bits = w*n - r;
0429 detail::polynomial helper;
0430 helper[num_bits - 1] = 1;
0431 mersenne_twister_engine tmp;
0432 tmp.reconstruct_state(helper);
0433
0434
0435 for(std::size_t j = 0; j < num_bits; ++j) {
0436 if(tmp.i >= n) tmp.twist();
0437 if(j == num_bits - 1)
0438 assert((tmp.x[tmp.i] & 1) == 1);
0439 else
0440 assert((tmp.x[tmp.i] & 1) == 0);
0441 ++tmp.i;
0442 }
0443 detail::polynomial phi;
0444 phi[num_bits] = 1;
0445 detail::polynomial next_bits = tmp.as_polynomial(num_bits);
0446 for(std::size_t j = 0; j < num_bits; ++j) {
0447 int val = next_bits[j] ^ phi[num_bits-j-1];
0448 phi[num_bits-j-1] = val;
0449 if(val) {
0450 for(std::size_t k = j + 1; k < num_bits; ++k) {
0451 phi[num_bits-k-1] ^= next_bits[k-j-1];
0452 }
0453 }
0454 }
0455 return phi;
0456 }
0457 detail::polynomial as_polynomial(std::size_t size) {
0458 detail::polynomial result;
0459 for(std::size_t j = 0; j < size; ++j) {
0460 if(i >= n) twist();
0461 result[j] = x[i++] & UIntType(1);
0462 }
0463 return result;
0464 }
0465 void reconstruct_state(const detail::polynomial& p)
0466 {
0467 const UIntType upper_mask = (~static_cast<UIntType>(0)) << r;
0468 const UIntType lower_mask = ~upper_mask;
0469 const std::size_t num_bits = w*n - r;
0470 for(std::size_t j = num_bits - n + 1; j <= num_bits; ++j)
0471 x[j % n] = p[j];
0472
0473 UIntType y0 = 0;
0474 for(std::size_t j = num_bits + 1; j >= n - 1; --j) {
0475 UIntType y1 = x[j % n] ^ x[(j + m) % n];
0476 if(p[j - n + 1])
0477 y1 = (y1 ^ a) << UIntType(1) | UIntType(1);
0478 else
0479 y1 = y1 << UIntType(1);
0480 x[(j + 1) % n] = (y0 & upper_mask) | (y1 & lower_mask);
0481 y0 = y1;
0482 }
0483 i = 0;
0484 }
0485
0486
0487
0488
0489
0490
0491
0492 UIntType x[n];
0493 std::size_t i;
0494 };
0495
0496
0497
0498 #ifndef BOOST_NO_INCLASS_MEMBER_INITIALIZATION
0499
0500 #define BOOST_RANDOM_MT_DEFINE_CONSTANT(type, name) \
0501 template<class UIntType, std::size_t w, std::size_t n, std::size_t m, \
0502 std::size_t r, UIntType a, std::size_t u, UIntType d, std::size_t s, \
0503 UIntType b, std::size_t t, UIntType c, std::size_t l, UIntType f> \
0504 const type mersenne_twister_engine<UIntType,w,n,m,r,a,u,d,s,b,t,c,l,f>::name
0505 BOOST_RANDOM_MT_DEFINE_CONSTANT(std::size_t, word_size);
0506 BOOST_RANDOM_MT_DEFINE_CONSTANT(std::size_t, state_size);
0507 BOOST_RANDOM_MT_DEFINE_CONSTANT(std::size_t, shift_size);
0508 BOOST_RANDOM_MT_DEFINE_CONSTANT(std::size_t, mask_bits);
0509 BOOST_RANDOM_MT_DEFINE_CONSTANT(UIntType, xor_mask);
0510 BOOST_RANDOM_MT_DEFINE_CONSTANT(std::size_t, tempering_u);
0511 BOOST_RANDOM_MT_DEFINE_CONSTANT(UIntType, tempering_d);
0512 BOOST_RANDOM_MT_DEFINE_CONSTANT(std::size_t, tempering_s);
0513 BOOST_RANDOM_MT_DEFINE_CONSTANT(UIntType, tempering_b);
0514 BOOST_RANDOM_MT_DEFINE_CONSTANT(std::size_t, tempering_t);
0515 BOOST_RANDOM_MT_DEFINE_CONSTANT(UIntType, tempering_c);
0516 BOOST_RANDOM_MT_DEFINE_CONSTANT(std::size_t, tempering_l);
0517 BOOST_RANDOM_MT_DEFINE_CONSTANT(UIntType, initialization_multiplier);
0518 BOOST_RANDOM_MT_DEFINE_CONSTANT(UIntType, default_seed);
0519 BOOST_RANDOM_MT_DEFINE_CONSTANT(UIntType, parameter_a);
0520 BOOST_RANDOM_MT_DEFINE_CONSTANT(std::size_t, output_u );
0521 BOOST_RANDOM_MT_DEFINE_CONSTANT(std::size_t, output_s);
0522 BOOST_RANDOM_MT_DEFINE_CONSTANT(UIntType, output_b);
0523 BOOST_RANDOM_MT_DEFINE_CONSTANT(std::size_t, output_t);
0524 BOOST_RANDOM_MT_DEFINE_CONSTANT(UIntType, output_c);
0525 BOOST_RANDOM_MT_DEFINE_CONSTANT(std::size_t, output_l);
0526 BOOST_RANDOM_MT_DEFINE_CONSTANT(bool, has_fixed_range);
0527 #undef BOOST_RANDOM_MT_DEFINE_CONSTANT
0528 #endif
0529
0530 template<class UIntType,
0531 std::size_t w, std::size_t n, std::size_t m, std::size_t r,
0532 UIntType a, std::size_t u, UIntType d, std::size_t s,
0533 UIntType b, std::size_t t,
0534 UIntType c, std::size_t l, UIntType f>
0535 void
0536 mersenne_twister_engine<UIntType,w,n,m,r,a,u,d,s,b,t,c,l,f>::twist()
0537 {
0538 const UIntType upper_mask = (~static_cast<UIntType>(0)) << r;
0539 const UIntType lower_mask = ~upper_mask;
0540
0541 const std::size_t unroll_factor = 6;
0542 const std::size_t unroll_extra1 = (n-m) % unroll_factor;
0543 const std::size_t unroll_extra2 = (m-1) % unroll_factor;
0544
0545
0546 {
0547 for(std::size_t j = 0; j < n-m-unroll_extra1; j++) {
0548 UIntType y = (x[j] & upper_mask) | (x[j+1] & lower_mask);
0549 x[j] = x[j+m] ^ (y >> 1) ^ ((x[j+1]&1) * a);
0550 }
0551 }
0552 {
0553 for(std::size_t j = n-m-unroll_extra1; j < n-m; j++) {
0554 UIntType y = (x[j] & upper_mask) | (x[j+1] & lower_mask);
0555 x[j] = x[j+m] ^ (y >> 1) ^ ((x[j+1]&1) * a);
0556 }
0557 }
0558 {
0559 for(std::size_t j = n-m; j < n-1-unroll_extra2; j++) {
0560 UIntType y = (x[j] & upper_mask) | (x[j+1] & lower_mask);
0561 x[j] = x[j-(n-m)] ^ (y >> 1) ^ ((x[j+1]&1) * a);
0562 }
0563 }
0564 {
0565 for(std::size_t j = n-1-unroll_extra2; j < n-1; j++) {
0566 UIntType y = (x[j] & upper_mask) | (x[j+1] & lower_mask);
0567 x[j] = x[j-(n-m)] ^ (y >> 1) ^ ((x[j+1]&1) * a);
0568 }
0569 }
0570
0571 UIntType y = (x[n-1] & upper_mask) | (x[0] & lower_mask);
0572 x[n-1] = x[m-1] ^ (y >> 1) ^ ((x[0]&1) * a);
0573 i = 0;
0574 }
0575
0576
0577 template<class UIntType,
0578 std::size_t w, std::size_t n, std::size_t m, std::size_t r,
0579 UIntType a, std::size_t u, UIntType d, std::size_t s,
0580 UIntType b, std::size_t t,
0581 UIntType c, std::size_t l, UIntType f>
0582 inline typename
0583 mersenne_twister_engine<UIntType,w,n,m,r,a,u,d,s,b,t,c,l,f>::result_type
0584 mersenne_twister_engine<UIntType,w,n,m,r,a,u,d,s,b,t,c,l,f>::operator()()
0585 {
0586 if(i == n)
0587 twist();
0588
0589 UIntType z = x[i];
0590 ++i;
0591 z ^= ((z >> u) & d);
0592 z ^= ((z << s) & b);
0593 z ^= ((z << t) & c);
0594 z ^= (z >> l);
0595 return z;
0596 }
0597
0598
0599
0600
0601
0602
0603
0604
0605
0606
0607
0608
0609 typedef mersenne_twister_engine<uint32_t,32,351,175,19,0xccab8ee7,
0610 11,0xffffffff,7,0x31b6ab00,15,0xffe50000,17,1812433253> mt11213b;
0611
0612
0613
0614
0615
0616
0617
0618
0619
0620
0621
0622
0623 typedef mersenne_twister_engine<uint32_t,32,624,397,31,0x9908b0df,
0624 11,0xffffffff,7,0x9d2c5680,15,0xefc60000,18,1812433253> mt19937;
0625
0626 #if !defined(BOOST_NO_INT64_T) && !defined(BOOST_NO_INTEGRAL_INT64_T)
0627 typedef mersenne_twister_engine<uint64_t,64,312,156,31,
0628 UINT64_C(0xb5026f5aa96619e9),29,UINT64_C(0x5555555555555555),17,
0629 UINT64_C(0x71d67fffeda60000),37,UINT64_C(0xfff7eee000000000),43,
0630 UINT64_C(6364136223846793005)> mt19937_64;
0631 #endif
0632
0633
0634
0635 template<class UIntType,
0636 int w, int n, int m, int r,
0637 UIntType a, int u, std::size_t s,
0638 UIntType b, int t,
0639 UIntType c, int l, UIntType v>
0640 class mersenne_twister :
0641 public mersenne_twister_engine<UIntType,
0642 w, n, m, r, a, u, ~(UIntType)0, s, b, t, c, l, 1812433253>
0643 {
0644 typedef mersenne_twister_engine<UIntType,
0645 w, n, m, r, a, u, ~(UIntType)0, s, b, t, c, l, 1812433253> base_type;
0646 public:
0647 mersenne_twister() {}
0648 BOOST_RANDOM_DETAIL_GENERATOR_CONSTRUCTOR(mersenne_twister, Gen, gen)
0649 { seed(gen); }
0650 BOOST_RANDOM_DETAIL_ARITHMETIC_CONSTRUCTOR(mersenne_twister, UIntType, val)
0651 { seed(val); }
0652 template<class It>
0653 mersenne_twister(It& first, It last) : base_type(first, last) {}
0654 void seed() { base_type::seed(); }
0655 BOOST_RANDOM_DETAIL_GENERATOR_SEED(mersenne_twister, Gen, gen)
0656 {
0657 detail::generator_seed_seq<Gen> seq(gen);
0658 base_type::seed(seq);
0659 }
0660 BOOST_RANDOM_DETAIL_ARITHMETIC_SEED(mersenne_twister, UIntType, val)
0661 { base_type::seed(val); }
0662 template<class It>
0663 void seed(It& first, It last) { base_type::seed(first, last); }
0664 };
0665
0666
0667
0668 }
0669
0670 using random::mt11213b;
0671 using random::mt19937;
0672 using random::mt19937_64;
0673
0674 }
0675
0676 BOOST_RANDOM_PTR_HELPER_SPEC(boost::mt11213b)
0677 BOOST_RANDOM_PTR_HELPER_SPEC(boost::mt19937)
0678 BOOST_RANDOM_PTR_HELPER_SPEC(boost::mt19937_64)
0679
0680 #include <boost/random/detail/enable_warnings.hpp>
0681
0682 #endif