File indexing completed on 2025-01-18 09:51:11
0001
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
0016 #ifndef BOOST_RANDOM_LINEAR_CONGRUENTIAL_HPP
0017 #define BOOST_RANDOM_LINEAR_CONGRUENTIAL_HPP
0018
0019 #include <iostream>
0020 #include <stdexcept>
0021 #include <boost/assert.hpp>
0022 #include <boost/config.hpp>
0023 #include <boost/cstdint.hpp>
0024 #include <boost/limits.hpp>
0025 #include <boost/static_assert.hpp>
0026 #include <boost/type_traits/is_arithmetic.hpp>
0027 #include <boost/random/detail/config.hpp>
0028 #include <boost/random/detail/const_mod.hpp>
0029 #include <boost/random/detail/seed.hpp>
0030 #include <boost/random/detail/seed_impl.hpp>
0031 #include <boost/detail/workaround.hpp>
0032
0033 #include <boost/random/detail/disable_warnings.hpp>
0034
0035 namespace boost {
0036 namespace random {
0037
0038
0039
0040
0041
0042
0043
0044
0045
0046
0047
0048
0049
0050
0051
0052
0053
0054
0055
0056
0057
0058
0059
0060 template<class IntType, IntType a, IntType c, IntType m>
0061 class linear_congruential_engine
0062 {
0063 public:
0064 typedef IntType result_type;
0065
0066
0067 BOOST_STATIC_CONSTANT(bool, has_fixed_range = false);
0068
0069 BOOST_STATIC_CONSTANT(IntType, multiplier = a);
0070 BOOST_STATIC_CONSTANT(IntType, increment = c);
0071 BOOST_STATIC_CONSTANT(IntType, modulus = m);
0072 BOOST_STATIC_CONSTANT(IntType, default_seed = 1);
0073
0074 BOOST_STATIC_ASSERT(std::numeric_limits<IntType>::is_integer);
0075 BOOST_STATIC_ASSERT(m == 0 || a < m);
0076 BOOST_STATIC_ASSERT(m == 0 || c < m);
0077
0078
0079
0080
0081 linear_congruential_engine() { seed(); }
0082
0083
0084
0085
0086 BOOST_RANDOM_DETAIL_ARITHMETIC_CONSTRUCTOR(linear_congruential_engine,
0087 IntType, x0)
0088 { seed(x0); }
0089
0090
0091
0092
0093
0094 BOOST_RANDOM_DETAIL_SEED_SEQ_CONSTRUCTOR(linear_congruential_engine,
0095 SeedSeq, seq)
0096 { seed(seq); }
0097
0098
0099
0100
0101
0102
0103
0104
0105
0106 template<class It>
0107 linear_congruential_engine(It& first, It last)
0108 {
0109 seed(first, last);
0110 }
0111
0112
0113
0114
0115
0116
0117 void seed() { seed(default_seed); }
0118
0119
0120
0121
0122
0123
0124
0125 BOOST_RANDOM_DETAIL_ARITHMETIC_SEED(linear_congruential_engine, IntType, x0_)
0126 {
0127
0128
0129
0130 IntType x0 = x0_;
0131
0132 if(modulus == 0) {
0133 _x = x0;
0134 } else {
0135 _x = x0 % modulus;
0136 }
0137
0138 if(_x < 0) {
0139 _x += modulus;
0140 }
0141
0142 if(increment == 0 && _x == 0) {
0143 _x = 1;
0144 }
0145 BOOST_ASSERT(_x >= (min)());
0146 BOOST_ASSERT(_x <= (max)());
0147 }
0148
0149
0150
0151
0152 BOOST_RANDOM_DETAIL_SEED_SEQ_SEED(linear_congruential_engine, SeedSeq, seq)
0153 { seed(detail::seed_one_int<IntType, m>(seq)); }
0154
0155
0156
0157
0158
0159
0160
0161
0162
0163 template<class It>
0164 void seed(It& first, It last)
0165 { seed(detail::get_one_int<IntType, m>(first, last)); }
0166
0167
0168
0169
0170
0171 static BOOST_CONSTEXPR result_type min BOOST_PREVENT_MACRO_SUBSTITUTION ()
0172 { return c == 0 ? 1 : 0; }
0173
0174
0175
0176
0177 static BOOST_CONSTEXPR result_type max BOOST_PREVENT_MACRO_SUBSTITUTION ()
0178 { return modulus-1; }
0179
0180
0181 IntType operator()()
0182 {
0183 _x = const_mod<IntType, m>::mult_add(a, _x, c);
0184 return _x;
0185 }
0186
0187
0188 template<class Iter>
0189 void generate(Iter first, Iter last)
0190 { detail::generate_from_int(*this, first, last); }
0191
0192
0193 void discard(boost::uintmax_t z)
0194 {
0195 typedef const_mod<IntType, m> mod_type;
0196 IntType b_inv = mod_type::invert(a-1);
0197 IntType b_gcd = mod_type::mult(a-1, b_inv);
0198 if(b_gcd == 1) {
0199 IntType a_z = mod_type::pow(a, z);
0200 _x = mod_type::mult_add(a_z, _x,
0201 mod_type::mult(mod_type::mult(c, b_inv), a_z - 1));
0202 } else {
0203
0204
0205 IntType a_zm1_over_gcd = 0;
0206 IntType a_km1_over_gcd = (a - 1) / b_gcd;
0207 boost::uintmax_t exponent = z;
0208 while(exponent != 0) {
0209 if(exponent % 2 == 1) {
0210 a_zm1_over_gcd =
0211 mod_type::mult_add(
0212 b_gcd,
0213 mod_type::mult(a_zm1_over_gcd, a_km1_over_gcd),
0214 mod_type::add(a_zm1_over_gcd, a_km1_over_gcd));
0215 }
0216 a_km1_over_gcd = mod_type::mult_add(
0217 b_gcd,
0218 mod_type::mult(a_km1_over_gcd, a_km1_over_gcd),
0219 mod_type::add(a_km1_over_gcd, a_km1_over_gcd));
0220 exponent /= 2;
0221 }
0222
0223 IntType a_z = mod_type::mult_add(b_gcd, a_zm1_over_gcd, 1);
0224 IntType num = mod_type::mult(c, a_zm1_over_gcd);
0225 b_inv = mod_type::invert((a-1)/b_gcd);
0226 _x = mod_type::mult_add(a_z, _x, mod_type::mult(b_inv, num));
0227 }
0228 }
0229
0230 friend bool operator==(const linear_congruential_engine& x,
0231 const linear_congruential_engine& y)
0232 { return x._x == y._x; }
0233 friend bool operator!=(const linear_congruential_engine& x,
0234 const linear_congruential_engine& y)
0235 { return !(x == y); }
0236
0237 #if !defined(BOOST_RANDOM_NO_STREAM_OPERATORS)
0238
0239 template<class CharT, class Traits>
0240 friend std::basic_ostream<CharT,Traits>&
0241 operator<<(std::basic_ostream<CharT,Traits>& os,
0242 const linear_congruential_engine& lcg)
0243 {
0244 return os << lcg._x;
0245 }
0246
0247
0248 template<class CharT, class Traits>
0249 friend std::basic_istream<CharT,Traits>&
0250 operator>>(std::basic_istream<CharT,Traits>& is,
0251 linear_congruential_engine& lcg)
0252 {
0253 lcg.read(is);
0254 return is;
0255 }
0256 #endif
0257
0258 private:
0259
0260
0261
0262 template<class CharT, class Traits>
0263 void read(std::basic_istream<CharT, Traits>& is) {
0264 IntType x;
0265 if(is >> x) {
0266 if(x >= (min)() && x <= (max)()) {
0267 _x = x;
0268 } else {
0269 is.setstate(std::ios_base::failbit);
0270 }
0271 }
0272 }
0273
0274
0275
0276 IntType _x;
0277 };
0278
0279 #ifndef BOOST_NO_INCLASS_MEMBER_INITIALIZATION
0280
0281 template<class IntType, IntType a, IntType c, IntType m>
0282 const bool linear_congruential_engine<IntType, a, c, m>::has_fixed_range;
0283 template<class IntType, IntType a, IntType c, IntType m>
0284 const IntType linear_congruential_engine<IntType,a,c,m>::multiplier;
0285 template<class IntType, IntType a, IntType c, IntType m>
0286 const IntType linear_congruential_engine<IntType,a,c,m>::increment;
0287 template<class IntType, IntType a, IntType c, IntType m>
0288 const IntType linear_congruential_engine<IntType,a,c,m>::modulus;
0289 template<class IntType, IntType a, IntType c, IntType m>
0290 const IntType linear_congruential_engine<IntType,a,c,m>::default_seed;
0291 #endif
0292
0293
0294
0295
0296 template<class IntType, IntType a, IntType c, IntType m, IntType val = 0>
0297 class linear_congruential : public linear_congruential_engine<IntType, a, c, m>
0298 {
0299 typedef linear_congruential_engine<IntType, a, c, m> base_type;
0300 public:
0301 linear_congruential(IntType x0 = 1) : base_type(x0) {}
0302 template<class It>
0303 linear_congruential(It& first, It last) : base_type(first, last) {}
0304 };
0305
0306
0307
0308
0309
0310
0311
0312
0313
0314
0315
0316
0317
0318
0319
0320
0321
0322
0323
0324
0325 typedef linear_congruential_engine<uint32_t, 16807, 0, 2147483647> minstd_rand0;
0326
0327
0328
0329
0330
0331
0332
0333
0334
0335 typedef linear_congruential_engine<uint32_t, 48271, 0, 2147483647> minstd_rand;
0336
0337
0338 #if !defined(BOOST_NO_INT64_T) && !defined(BOOST_NO_INTEGRAL_INT64_T)
0339
0340
0341
0342
0343
0344
0345
0346
0347
0348
0349 class rand48
0350 {
0351 public:
0352 typedef boost::uint32_t result_type;
0353
0354 BOOST_STATIC_CONSTANT(bool, has_fixed_range = false);
0355
0356
0357
0358 static BOOST_CONSTEXPR uint32_t min BOOST_PREVENT_MACRO_SUBSTITUTION () { return 0; }
0359
0360
0361
0362 static BOOST_CONSTEXPR uint32_t max BOOST_PREVENT_MACRO_SUBSTITUTION ()
0363 { return 0x7FFFFFFF; }
0364
0365
0366 rand48() : lcf(cnv(static_cast<uint32_t>(1))) {}
0367
0368
0369
0370 BOOST_RANDOM_DETAIL_ARITHMETIC_CONSTRUCTOR(rand48, result_type, x0)
0371 { seed(x0); }
0372
0373
0374
0375 BOOST_RANDOM_DETAIL_SEED_SEQ_CONSTRUCTOR(rand48, SeedSeq, seq)
0376 { seed(seq); }
0377
0378
0379
0380
0381 template<class It> rand48(It& first, It last) : lcf(first, last) { }
0382
0383
0384
0385
0386 void seed() { seed(static_cast<uint32_t>(1)); }
0387
0388
0389
0390 BOOST_RANDOM_DETAIL_ARITHMETIC_SEED(rand48, result_type, x0)
0391 { lcf.seed(cnv(x0)); }
0392
0393
0394
0395
0396 template<class It> void seed(It& first, It last) { lcf.seed(first,last); }
0397
0398
0399
0400 BOOST_RANDOM_DETAIL_SEED_SEQ_SEED(rand48, SeedSeq, seq)
0401 { lcf.seed(seq); }
0402
0403
0404 uint32_t operator()() { return static_cast<uint32_t>(lcf() >> 17); }
0405
0406
0407 void discard(boost::uintmax_t z) { lcf.discard(z); }
0408
0409
0410 template<class Iter>
0411 void generate(Iter first, Iter last)
0412 {
0413 for(; first != last; ++first) {
0414 *first = (*this)();
0415 }
0416 }
0417
0418 #ifndef BOOST_RANDOM_NO_STREAM_OPERATORS
0419
0420 template<class CharT,class Traits>
0421 friend std::basic_ostream<CharT,Traits>&
0422 operator<<(std::basic_ostream<CharT,Traits>& os, const rand48& r)
0423 { os << r.lcf; return os; }
0424
0425
0426 template<class CharT,class Traits>
0427 friend std::basic_istream<CharT,Traits>&
0428 operator>>(std::basic_istream<CharT,Traits>& is, rand48& r)
0429 { is >> r.lcf; return is; }
0430 #endif
0431
0432
0433
0434
0435
0436 friend bool operator==(const rand48& x, const rand48& y)
0437 { return x.lcf == y.lcf; }
0438
0439
0440
0441
0442 friend bool operator!=(const rand48& x, const rand48& y)
0443 { return !(x == y); }
0444 private:
0445
0446 typedef random::linear_congruential_engine<uint64_t,
0447
0448 uint64_t(0xDEECE66DUL) | (uint64_t(0x5) << 32),
0449 0xB, uint64_t(1)<<48> lcf_t;
0450 lcf_t lcf;
0451
0452 static boost::uint64_t cnv(boost::uint32_t x)
0453 { return (static_cast<uint64_t>(x) << 16) | 0x330e; }
0454
0455 };
0456 #endif
0457
0458 }
0459
0460 using random::minstd_rand0;
0461 using random::minstd_rand;
0462 using random::rand48;
0463
0464 }
0465
0466 #include <boost/random/detail/enable_warnings.hpp>
0467
0468 #endif