1 /* boost random/linear_congruential.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$ 11 * 12 * Revision history 13 * 2001-02-18 moved to individual header files 14 */ 15 16 #ifndef BOOST_RANDOM_LINEAR_CONGRUENTIAL_HPP 17 #define BOOST_RANDOM_LINEAR_CONGRUENTIAL_HPP 18 19 #include <iostream> 20 #include <stdexcept> 21 #include <boost/assert.hpp> 22 #include <boost/config.hpp> 23 #include <boost/cstdint.hpp> 24 #include <boost/limits.hpp> 25 #include <boost/static_assert.hpp> 26 #include <boost/integer/static_log2.hpp> 27 #include <boost/mpl/if.hpp> 28 #include <boost/type_traits/is_arithmetic.hpp> 29 #include <boost/random/detail/config.hpp> 30 #include <boost/random/detail/const_mod.hpp> 31 #include <boost/random/detail/seed.hpp> 32 #include <boost/random/detail/seed_impl.hpp> 33 #include <boost/detail/workaround.hpp> 34 35 #include <boost/random/detail/disable_warnings.hpp> 36 37 namespace boost { 38 namespace random { 39 40 /** 41 * Instantiations of class template linear_congruential_engine model a 42 * \pseudo_random_number_generator. Linear congruential pseudo-random 43 * number generators are described in: 44 * 45 * @blockquote 46 * "Mathematical methods in large-scale computing units", D. H. Lehmer, 47 * Proc. 2nd Symposium on Large-Scale Digital Calculating Machines, 48 * Harvard University Press, 1951, pp. 141-146 49 * @endblockquote 50 * 51 * Let x(n) denote the sequence of numbers returned by some pseudo-random 52 * number generator. Then for the linear congruential generator, 53 * x(n+1) := (a * x(n) + c) mod m. Parameters for the generator are 54 * x(0), a, c, m. The template parameter IntType shall denote an integral 55 * type. It must be large enough to hold values a, c, and m. The template 56 * parameters a and c must be smaller than m. 57 * 58 * Note: The quality of the generator crucially depends on the choice of 59 * the parameters. User code should use one of the sensibly parameterized 60 * generators such as minstd_rand instead. 61 */ 62 template<class IntType, IntType a, IntType c, IntType m> 63 class linear_congruential_engine 64 { 65 public: 66 typedef IntType result_type; 67 68 // Required for old Boost.Random concept 69 BOOST_STATIC_CONSTANT(bool, has_fixed_range = false); 70 71 BOOST_STATIC_CONSTANT(IntType, multiplier = a); 72 BOOST_STATIC_CONSTANT(IntType, increment = c); 73 BOOST_STATIC_CONSTANT(IntType, modulus = m); 74 BOOST_STATIC_CONSTANT(IntType, default_seed = 1); 75 76 BOOST_STATIC_ASSERT(std::numeric_limits<IntType>::is_integer); 77 BOOST_STATIC_ASSERT(m == 0 || a < m); 78 BOOST_STATIC_ASSERT(m == 0 || c < m); 79 80 /** 81 * Constructs a @c linear_congruential_engine, using the default seed 82 */ linear_congruential_engine()83 linear_congruential_engine() { seed(); } 84 85 /** 86 * Constructs a @c linear_congruential_engine, seeding it with @c x0. 87 */ BOOST_RANDOM_DETAIL_ARITHMETIC_CONSTRUCTOR(linear_congruential_engine,IntType,x0)88 BOOST_RANDOM_DETAIL_ARITHMETIC_CONSTRUCTOR(linear_congruential_engine, 89 IntType, x0) 90 { seed(x0); } 91 92 /** 93 * Constructs a @c linear_congruential_engine, seeding it with values 94 * produced by a call to @c seq.generate(). 95 */ BOOST_RANDOM_DETAIL_SEED_SEQ_CONSTRUCTOR(linear_congruential_engine,SeedSeq,seq)96 BOOST_RANDOM_DETAIL_SEED_SEQ_CONSTRUCTOR(linear_congruential_engine, 97 SeedSeq, seq) 98 { seed(seq); } 99 100 /** 101 * Constructs a @c linear_congruential_engine and seeds it 102 * with values taken from the itrator range [first, last) 103 * and adjusts first to point to the element after the last one 104 * used. If there are not enough elements, throws @c std::invalid_argument. 105 * 106 * first and last must be input iterators. 107 */ 108 template<class It> linear_congruential_engine(It & first,It last)109 linear_congruential_engine(It& first, It last) 110 { 111 seed(first, last); 112 } 113 114 // compiler-generated copy constructor and assignment operator are fine 115 116 /** 117 * Calls seed(default_seed) 118 */ seed()119 void seed() { seed(default_seed); } 120 121 /** 122 * If c mod m is zero and x0 mod m is zero, changes the current value of 123 * the generator to 1. Otherwise, changes it to x0 mod m. If c is zero, 124 * distinct seeds in the range [1,m) will leave the generator in distinct 125 * states. If c is not zero, the range is [0,m). 126 */ BOOST_RANDOM_DETAIL_ARITHMETIC_SEED(linear_congruential_engine,IntType,x0_)127 BOOST_RANDOM_DETAIL_ARITHMETIC_SEED(linear_congruential_engine, IntType, x0_) 128 { 129 // Work around a msvc 12/14 optimizer bug, which causes 130 // the line _x = 1 to run unconditionally sometimes. 131 // Creating a local copy seems to make it work. 132 IntType x0 = x0_; 133 // wrap _x if it doesn't fit in the destination 134 if(modulus == 0) { 135 _x = x0; 136 } else { 137 _x = x0 % modulus; 138 } 139 // handle negative seeds 140 if(_x <= 0 && _x != 0) { 141 _x += modulus; 142 } 143 // adjust to the correct range 144 if(increment == 0 && _x == 0) { 145 _x = 1; 146 } 147 BOOST_ASSERT(_x >= (min)()); 148 BOOST_ASSERT(_x <= (max)()); 149 } 150 151 /** 152 * Seeds a @c linear_congruential_engine using values from a SeedSeq. 153 */ BOOST_RANDOM_DETAIL_SEED_SEQ_SEED(linear_congruential_engine,SeedSeq,seq)154 BOOST_RANDOM_DETAIL_SEED_SEQ_SEED(linear_congruential_engine, SeedSeq, seq) 155 { seed(detail::seed_one_int<IntType, m>(seq)); } 156 157 /** 158 * seeds a @c linear_congruential_engine with values taken 159 * from the itrator range [first, last) and adjusts @c first to 160 * point to the element after the last one used. If there are 161 * not enough elements, throws @c std::invalid_argument. 162 * 163 * @c first and @c last must be input iterators. 164 */ 165 template<class It> seed(It & first,It last)166 void seed(It& first, It last) 167 { seed(detail::get_one_int<IntType, m>(first, last)); } 168 169 /** 170 * Returns the smallest value that the @c linear_congruential_engine 171 * can produce. 172 */ BOOST_PREVENT_MACRO_SUBSTITUTION()173 static result_type min BOOST_PREVENT_MACRO_SUBSTITUTION () 174 { return c == 0 ? 1 : 0; } 175 /** 176 * Returns the largest value that the @c linear_congruential_engine 177 * can produce. 178 */ BOOST_PREVENT_MACRO_SUBSTITUTION()179 static result_type max BOOST_PREVENT_MACRO_SUBSTITUTION () 180 { return modulus-1; } 181 182 /** Returns the next value of the @c linear_congruential_engine. */ operator ()()183 IntType operator()() 184 { 185 _x = const_mod<IntType, m>::mult_add(a, _x, c); 186 return _x; 187 } 188 189 /** Fills a range with random values */ 190 template<class Iter> generate(Iter first,Iter last)191 void generate(Iter first, Iter last) 192 { detail::generate_from_int(*this, first, last); } 193 194 /** Advances the state of the generator by @c z. */ discard(boost::uintmax_t z)195 void discard(boost::uintmax_t z) 196 { 197 typedef const_mod<IntType, m> mod_type; 198 IntType b_inv = mod_type::invert(a-1); 199 IntType b_gcd = mod_type::mult(a-1, b_inv); 200 if(b_gcd == 1) { 201 IntType a_z = mod_type::pow(a, z); 202 _x = mod_type::mult_add(a_z, _x, 203 mod_type::mult(mod_type::mult(c, b_inv), a_z - 1)); 204 } else { 205 // compute (a^z - 1)*c % (b_gcd * m) / (b / b_gcd) * inv(b / b_gcd) 206 // we're storing the intermediate result / b_gcd 207 IntType a_zm1_over_gcd = 0; 208 IntType a_km1_over_gcd = (a - 1) / b_gcd; 209 boost::uintmax_t exponent = z; 210 while(exponent != 0) { 211 if(exponent % 2 == 1) { 212 a_zm1_over_gcd = 213 mod_type::mult_add( 214 b_gcd, 215 mod_type::mult(a_zm1_over_gcd, a_km1_over_gcd), 216 mod_type::add(a_zm1_over_gcd, a_km1_over_gcd)); 217 } 218 a_km1_over_gcd = mod_type::mult_add( 219 b_gcd, 220 mod_type::mult(a_km1_over_gcd, a_km1_over_gcd), 221 mod_type::add(a_km1_over_gcd, a_km1_over_gcd)); 222 exponent /= 2; 223 } 224 225 IntType a_z = mod_type::mult_add(b_gcd, a_zm1_over_gcd, 1); 226 IntType num = mod_type::mult(c, a_zm1_over_gcd); 227 b_inv = mod_type::invert((a-1)/b_gcd); 228 _x = mod_type::mult_add(a_z, _x, mod_type::mult(b_inv, num)); 229 } 230 } 231 operator ==(const linear_congruential_engine & x,const linear_congruential_engine & y)232 friend bool operator==(const linear_congruential_engine& x, 233 const linear_congruential_engine& y) 234 { return x._x == y._x; } operator !=(const linear_congruential_engine & x,const linear_congruential_engine & y)235 friend bool operator!=(const linear_congruential_engine& x, 236 const linear_congruential_engine& y) 237 { return !(x == y); } 238 239 #if !defined(BOOST_RANDOM_NO_STREAM_OPERATORS) 240 /** Writes a @c linear_congruential_engine to a @c std::ostream. */ 241 template<class CharT, class Traits> 242 friend std::basic_ostream<CharT,Traits>& operator <<(std::basic_ostream<CharT,Traits> & os,const linear_congruential_engine & lcg)243 operator<<(std::basic_ostream<CharT,Traits>& os, 244 const linear_congruential_engine& lcg) 245 { 246 return os << lcg._x; 247 } 248 249 /** Reads a @c linear_congruential_engine from a @c std::istream. */ 250 template<class CharT, class Traits> 251 friend std::basic_istream<CharT,Traits>& operator >>(std::basic_istream<CharT,Traits> & is,linear_congruential_engine & lcg)252 operator>>(std::basic_istream<CharT,Traits>& is, 253 linear_congruential_engine& lcg) 254 { 255 lcg.read(is); 256 return is; 257 } 258 #endif 259 260 private: 261 262 /// \cond show_private 263 264 template<class CharT, class Traits> read(std::basic_istream<CharT,Traits> & is)265 void read(std::basic_istream<CharT, Traits>& is) { 266 IntType x; 267 if(is >> x) { 268 if(x >= (min)() && x <= (max)()) { 269 _x = x; 270 } else { 271 is.setstate(std::ios_base::failbit); 272 } 273 } 274 } 275 276 /// \endcond 277 278 IntType _x; 279 }; 280 281 #ifndef BOOST_NO_INCLASS_MEMBER_INITIALIZATION 282 // A definition is required even for integral static constants 283 template<class IntType, IntType a, IntType c, IntType m> 284 const bool linear_congruential_engine<IntType, a, c, m>::has_fixed_range; 285 template<class IntType, IntType a, IntType c, IntType m> 286 const IntType linear_congruential_engine<IntType,a,c,m>::multiplier; 287 template<class IntType, IntType a, IntType c, IntType m> 288 const IntType linear_congruential_engine<IntType,a,c,m>::increment; 289 template<class IntType, IntType a, IntType c, IntType m> 290 const IntType linear_congruential_engine<IntType,a,c,m>::modulus; 291 template<class IntType, IntType a, IntType c, IntType m> 292 const IntType linear_congruential_engine<IntType,a,c,m>::default_seed; 293 #endif 294 295 /// \cond show_deprecated 296 297 // provided for backwards compatibility 298 template<class IntType, IntType a, IntType c, IntType m, IntType val = 0> 299 class linear_congruential : public linear_congruential_engine<IntType, a, c, m> 300 { 301 typedef linear_congruential_engine<IntType, a, c, m> base_type; 302 public: linear_congruential(IntType x0=1)303 linear_congruential(IntType x0 = 1) : base_type(x0) {} 304 template<class It> linear_congruential(It & first,It last)305 linear_congruential(It& first, It last) : base_type(first, last) {} 306 }; 307 308 /// \endcond 309 310 /** 311 * The specialization \minstd_rand0 was originally suggested in 312 * 313 * @blockquote 314 * A pseudo-random number generator for the System/360, P.A. Lewis, 315 * A.S. Goodman, J.M. Miller, IBM Systems Journal, Vol. 8, No. 2, 316 * 1969, pp. 136-146 317 * @endblockquote 318 * 319 * It is examined more closely together with \minstd_rand in 320 * 321 * @blockquote 322 * "Random Number Generators: Good ones are hard to find", 323 * Stephen K. Park and Keith W. Miller, Communications of 324 * the ACM, Vol. 31, No. 10, October 1988, pp. 1192-1201 325 * @endblockquote 326 */ 327 typedef linear_congruential_engine<uint32_t, 16807, 0, 2147483647> minstd_rand0; 328 329 /** The specialization \minstd_rand was suggested in 330 * 331 * @blockquote 332 * "Random Number Generators: Good ones are hard to find", 333 * Stephen K. Park and Keith W. Miller, Communications of 334 * the ACM, Vol. 31, No. 10, October 1988, pp. 1192-1201 335 * @endblockquote 336 */ 337 typedef linear_congruential_engine<uint32_t, 48271, 0, 2147483647> minstd_rand; 338 339 340 #if !defined(BOOST_NO_INT64_T) && !defined(BOOST_NO_INTEGRAL_INT64_T) 341 /** 342 * Class @c rand48 models a \pseudo_random_number_generator. It uses 343 * the linear congruential algorithm with the parameters a = 0x5DEECE66D, 344 * c = 0xB, m = 2**48. It delivers identical results to the @c lrand48() 345 * function available on some systems (assuming lcong48 has not been called). 346 * 347 * It is only available on systems where @c uint64_t is provided as an 348 * integral type, so that for example static in-class constants and/or 349 * enum definitions with large @c uint64_t numbers work. 350 */ 351 class rand48 352 { 353 public: 354 typedef boost::uint32_t result_type; 355 356 BOOST_STATIC_CONSTANT(bool, has_fixed_range = false); 357 /** 358 * Returns the smallest value that the generator can produce 359 */ BOOST_PREVENT_MACRO_SUBSTITUTION()360 static uint32_t min BOOST_PREVENT_MACRO_SUBSTITUTION () { return 0; } 361 /** 362 * Returns the largest value that the generator can produce 363 */ BOOST_PREVENT_MACRO_SUBSTITUTION()364 static uint32_t max BOOST_PREVENT_MACRO_SUBSTITUTION () 365 { return 0x7FFFFFFF; } 366 367 /** Seeds the generator with the default seed. */ rand48()368 rand48() : lcf(cnv(static_cast<uint32_t>(1))) {} 369 /** 370 * Constructs a \rand48 generator with x(0) := (x0 << 16) | 0x330e. 371 */ BOOST_RANDOM_DETAIL_ARITHMETIC_CONSTRUCTOR(rand48,result_type,x0)372 BOOST_RANDOM_DETAIL_ARITHMETIC_CONSTRUCTOR(rand48, result_type, x0) 373 { seed(x0); } 374 /** 375 * Seeds the generator with values produced by @c seq.generate(). 376 */ BOOST_RANDOM_DETAIL_SEED_SEQ_CONSTRUCTOR(rand48,SeedSeq,seq)377 BOOST_RANDOM_DETAIL_SEED_SEQ_CONSTRUCTOR(rand48, SeedSeq, seq) 378 { seed(seq); } 379 /** 380 * Seeds the generator using values from an iterator range, 381 * and updates first to point one past the last value consumed. 382 */ rand48(It & first,It last)383 template<class It> rand48(It& first, It last) : lcf(first, last) { } 384 385 // compiler-generated copy ctor and assignment operator are fine 386 387 /** Seeds the generator with the default seed. */ seed()388 void seed() { seed(static_cast<uint32_t>(1)); } 389 /** 390 * Changes the current value x(n) of the generator to (x0 << 16) | 0x330e. 391 */ BOOST_RANDOM_DETAIL_ARITHMETIC_SEED(rand48,result_type,x0)392 BOOST_RANDOM_DETAIL_ARITHMETIC_SEED(rand48, result_type, x0) 393 { lcf.seed(cnv(x0)); } 394 /** 395 * Seeds the generator using values from an iterator range, 396 * and updates first to point one past the last value consumed. 397 */ seed(It & first,It last)398 template<class It> void seed(It& first, It last) { lcf.seed(first,last); } 399 /** 400 * Seeds the generator with values produced by @c seq.generate(). 401 */ BOOST_RANDOM_DETAIL_SEED_SEQ_SEED(rand48,SeedSeq,seq)402 BOOST_RANDOM_DETAIL_SEED_SEQ_SEED(rand48, SeedSeq, seq) 403 { lcf.seed(seq); } 404 405 /** Returns the next value of the generator. */ operator ()()406 uint32_t operator()() { return static_cast<uint32_t>(lcf() >> 17); } 407 408 /** Advances the state of the generator by @c z. */ discard(boost::uintmax_t z)409 void discard(boost::uintmax_t z) { lcf.discard(z); } 410 411 /** Fills a range with random values */ 412 template<class Iter> generate(Iter first,Iter last)413 void generate(Iter first, Iter last) 414 { 415 for(; first != last; ++first) { 416 *first = (*this)(); 417 } 418 } 419 420 #ifndef BOOST_RANDOM_NO_STREAM_OPERATORS 421 /** Writes a @c rand48 to a @c std::ostream. */ 422 template<class CharT,class Traits> 423 friend std::basic_ostream<CharT,Traits>& operator <<(std::basic_ostream<CharT,Traits> & os,const rand48 & r)424 operator<<(std::basic_ostream<CharT,Traits>& os, const rand48& r) 425 { os << r.lcf; return os; } 426 427 /** Reads a @c rand48 from a @c std::istream. */ 428 template<class CharT,class Traits> 429 friend std::basic_istream<CharT,Traits>& operator >>(std::basic_istream<CharT,Traits> & is,rand48 & r)430 operator>>(std::basic_istream<CharT,Traits>& is, rand48& r) 431 { is >> r.lcf; return is; } 432 #endif 433 434 /** 435 * Returns true if the two generators will produce identical 436 * sequences of values. 437 */ operator ==(const rand48 & x,const rand48 & y)438 friend bool operator==(const rand48& x, const rand48& y) 439 { return x.lcf == y.lcf; } 440 /** 441 * Returns true if the two generators will produce different 442 * sequences of values. 443 */ operator !=(const rand48 & x,const rand48 & y)444 friend bool operator!=(const rand48& x, const rand48& y) 445 { return !(x == y); } 446 private: 447 /// \cond show_private 448 typedef random::linear_congruential_engine<uint64_t, 449 // xxxxULL is not portable 450 uint64_t(0xDEECE66DUL) | (uint64_t(0x5) << 32), 451 0xB, uint64_t(1)<<48> lcf_t; 452 lcf_t lcf; 453 cnv(boost::uint32_t x)454 static boost::uint64_t cnv(boost::uint32_t x) 455 { return (static_cast<uint64_t>(x) << 16) | 0x330e; } 456 /// \endcond 457 }; 458 #endif /* !BOOST_NO_INT64_T && !BOOST_NO_INTEGRAL_INT64_T */ 459 460 } // namespace random 461 462 using random::minstd_rand0; 463 using random::minstd_rand; 464 using random::rand48; 465 466 } // namespace boost 467 468 #include <boost/random/detail/enable_warnings.hpp> 469 470 #endif // BOOST_RANDOM_LINEAR_CONGRUENTIAL_HPP 471