lagged_fibonacci.hpp 18 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525526527528529530531532533534535536537
  1. /* boost random/lagged_fibonacci.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. * 2013-10-14 fixed some warnings with Wshadow (mgaunard)
  14. * 2001-02-18 moved to individual header files
  15. */
  16. #ifndef BOOST_RANDOM_LAGGED_FIBONACCI_HPP
  17. #define BOOST_RANDOM_LAGGED_FIBONACCI_HPP
  18. #include <istream>
  19. #include <iosfwd>
  20. #include <algorithm> // std::max
  21. #include <iterator>
  22. #include <boost/config/no_tr1/cmath.hpp> // std::pow
  23. #include <boost/config.hpp>
  24. #include <boost/limits.hpp>
  25. #include <boost/cstdint.hpp>
  26. #include <boost/integer/integer_mask.hpp>
  27. #include <boost/random/linear_congruential.hpp>
  28. #include <boost/random/uniform_01.hpp>
  29. #include <boost/random/detail/config.hpp>
  30. #include <boost/random/detail/seed.hpp>
  31. #include <boost/random/detail/operators.hpp>
  32. #include <boost/random/detail/generator_seed_seq.hpp>
  33. namespace boost {
  34. namespace random {
  35. /**
  36. * Instantiations of class template \lagged_fibonacci_engine model a
  37. * \pseudo_random_number_generator. It uses a lagged Fibonacci
  38. * algorithm with two lags @c p and @c q:
  39. * x(i) = x(i-p) + x(i-q) (mod 2<sup>w</sup>) with p > q.
  40. */
  41. template<class UIntType, int w, unsigned int p, unsigned int q>
  42. class lagged_fibonacci_engine
  43. {
  44. public:
  45. typedef UIntType result_type;
  46. BOOST_STATIC_CONSTANT(bool, has_fixed_range = false);
  47. BOOST_STATIC_CONSTANT(int, word_size = w);
  48. BOOST_STATIC_CONSTANT(unsigned int, long_lag = p);
  49. BOOST_STATIC_CONSTANT(unsigned int, short_lag = q);
  50. BOOST_STATIC_CONSTANT(UIntType, default_seed = 331u);
  51. /** Returns the smallest value that the generator can produce. */
  52. static result_type min BOOST_PREVENT_MACRO_SUBSTITUTION () { return 0; }
  53. /** Returns the largest value that the generator can produce. */
  54. static result_type max BOOST_PREVENT_MACRO_SUBSTITUTION ()
  55. { return low_bits_mask_t<w>::sig_bits; }
  56. /** Creates a new @c lagged_fibonacci_engine and calls @c seed(). */
  57. lagged_fibonacci_engine() { seed(); }
  58. /** Creates a new @c lagged_fibonacci_engine and calls @c seed(value). */
  59. BOOST_RANDOM_DETAIL_ARITHMETIC_CONSTRUCTOR(lagged_fibonacci_engine,
  60. UIntType, value)
  61. { seed(value); }
  62. /** Creates a new @c lagged_fibonacci_engine and calls @c seed(seq). */
  63. BOOST_RANDOM_DETAIL_SEED_SEQ_CONSTRUCTOR(lagged_fibonacci_engine,
  64. SeedSeq, seq)
  65. { seed(seq); }
  66. /**
  67. * Creates a new @c lagged_fibonacci_engine and calls @c seed(first, last).
  68. */
  69. template<class It> lagged_fibonacci_engine(It& first, It last)
  70. { seed(first, last); }
  71. // compiler-generated copy ctor and assignment operator are fine
  72. /** Calls @c seed(default_seed). */
  73. void seed() { seed(default_seed); }
  74. /**
  75. * Sets the state of the generator to values produced by
  76. * a \minstd_rand0 generator.
  77. */
  78. BOOST_RANDOM_DETAIL_ARITHMETIC_SEED(lagged_fibonacci_engine,
  79. UIntType, value)
  80. {
  81. minstd_rand0 intgen(static_cast<boost::uint32_t>(value));
  82. detail::generator_seed_seq<minstd_rand0> gen(intgen);
  83. seed(gen);
  84. }
  85. /**
  86. * Sets the state of the generator using values produced by seq.
  87. */
  88. BOOST_RANDOM_DETAIL_SEED_SEQ_SEED(lagged_fibonacci_engine, SeedSeq, seq)
  89. {
  90. detail::seed_array_int<w>(seq, x);
  91. i = long_lag;
  92. }
  93. /**
  94. * Sets the state of the generator to values from the iterator
  95. * range [first, last). If there are not enough elements in the
  96. * range [first, last) throws @c std::invalid_argument.
  97. */
  98. template<class It>
  99. void seed(It& first, It last)
  100. {
  101. detail::fill_array_int<w>(first, last, x);
  102. i = long_lag;
  103. }
  104. /** Returns the next value of the generator. */
  105. result_type operator()()
  106. {
  107. if(i >= long_lag)
  108. fill();
  109. return x[i++];
  110. }
  111. /** Fills a range with random values */
  112. template<class Iter>
  113. void generate(Iter first, Iter last)
  114. { detail::generate_from_int(*this, first, last); }
  115. /** Advances the state of the generator by @c z. */
  116. void discard(boost::uintmax_t z)
  117. {
  118. for(boost::uintmax_t j = 0; j < z; ++j) {
  119. (*this)();
  120. }
  121. }
  122. /**
  123. * Writes the textual representation of the generator to a @c std::ostream.
  124. */
  125. BOOST_RANDOM_DETAIL_OSTREAM_OPERATOR(os, lagged_fibonacci_engine, f)
  126. {
  127. os << f.i;
  128. for(unsigned int j = 0; j < f.long_lag; ++j)
  129. os << ' ' << f.x[j];
  130. return os;
  131. }
  132. /**
  133. * Reads the textual representation of the generator from a @c std::istream.
  134. */
  135. BOOST_RANDOM_DETAIL_ISTREAM_OPERATOR(is, lagged_fibonacci_engine, f)
  136. {
  137. is >> f.i >> std::ws;
  138. for(unsigned int j = 0; j < f.long_lag; ++j)
  139. is >> f.x[j] >> std::ws;
  140. return is;
  141. }
  142. /**
  143. * Returns true if the two generators will produce identical
  144. * sequences of outputs.
  145. */
  146. BOOST_RANDOM_DETAIL_EQUALITY_OPERATOR(lagged_fibonacci_engine, x_, y_)
  147. { return x_.i == y_.i && std::equal(x_.x, x_.x+long_lag, y_.x); }
  148. /**
  149. * Returns true if the two generators will produce different
  150. * sequences of outputs.
  151. */
  152. BOOST_RANDOM_DETAIL_INEQUALITY_OPERATOR(lagged_fibonacci_engine)
  153. private:
  154. /// \cond show_private
  155. void fill();
  156. /// \endcond
  157. unsigned int i;
  158. UIntType x[long_lag];
  159. };
  160. #ifndef BOOST_NO_INCLASS_MEMBER_INITIALIZATION
  161. // A definition is required even for integral static constants
  162. template<class UIntType, int w, unsigned int p, unsigned int q>
  163. const bool lagged_fibonacci_engine<UIntType, w, p, q>::has_fixed_range;
  164. template<class UIntType, int w, unsigned int p, unsigned int q>
  165. const unsigned int lagged_fibonacci_engine<UIntType, w, p, q>::long_lag;
  166. template<class UIntType, int w, unsigned int p, unsigned int q>
  167. const unsigned int lagged_fibonacci_engine<UIntType, w, p, q>::short_lag;
  168. template<class UIntType, int w, unsigned int p, unsigned int q>
  169. const UIntType lagged_fibonacci_engine<UIntType, w, p, q>::default_seed;
  170. #endif
  171. /// \cond show_private
  172. template<class UIntType, int w, unsigned int p, unsigned int q>
  173. void lagged_fibonacci_engine<UIntType, w, p, q>::fill()
  174. {
  175. // two loops to avoid costly modulo operations
  176. { // extra scope for MSVC brokenness w.r.t. for scope
  177. for(unsigned int j = 0; j < short_lag; ++j)
  178. x[j] = (x[j] + x[j+(long_lag-short_lag)]) & low_bits_mask_t<w>::sig_bits;
  179. }
  180. for(unsigned int j = short_lag; j < long_lag; ++j)
  181. x[j] = (x[j] + x[j-short_lag]) & low_bits_mask_t<w>::sig_bits;
  182. i = 0;
  183. }
  184. /// \endcond
  185. /// \cond show_deprecated
  186. // provided for backwards compatibility
  187. template<class UIntType, int w, unsigned int p, unsigned int q, UIntType v = 0>
  188. class lagged_fibonacci : public lagged_fibonacci_engine<UIntType, w, p, q>
  189. {
  190. typedef lagged_fibonacci_engine<UIntType, w, p, q> base_type;
  191. public:
  192. lagged_fibonacci() {}
  193. BOOST_RANDOM_DETAIL_ARITHMETIC_CONSTRUCTOR(lagged_fibonacci, UIntType, val)
  194. { this->seed(val); }
  195. BOOST_RANDOM_DETAIL_SEED_SEQ_CONSTRUCTOR(lagged_fibonacci, SeedSeq, seq)
  196. { this->seed(seq); }
  197. template<class It>
  198. lagged_fibonacci(It& first, It last) : base_type(first, last) {}
  199. };
  200. /// \endcond
  201. // lagged Fibonacci generator for the range [0..1)
  202. // contributed by Matthias Troyer
  203. // for p=55, q=24 originally by G. J. Mitchell and D. P. Moore 1958
  204. /**
  205. * Instantiations of class template @c lagged_fibonacci_01 model a
  206. * \pseudo_random_number_generator. It uses a lagged Fibonacci
  207. * algorithm with two lags @c p and @c q, evaluated in floating-point
  208. * arithmetic: x(i) = x(i-p) + x(i-q) (mod 1) with p > q. See
  209. *
  210. * @blockquote
  211. * "Uniform random number generators for supercomputers", Richard Brent,
  212. * Proc. of Fifth Australian Supercomputer Conference, Melbourne,
  213. * Dec. 1992, pp. 704-706.
  214. * @endblockquote
  215. *
  216. * @xmlnote
  217. * The quality of the generator crucially depends on the choice
  218. * of the parameters. User code should employ one of the sensibly
  219. * parameterized generators such as \lagged_fibonacci607 instead.
  220. * @endxmlnote
  221. *
  222. * The generator requires considerable amounts of memory for the storage
  223. * of its state array. For example, \lagged_fibonacci607 requires about
  224. * 4856 bytes and \lagged_fibonacci44497 requires about 350 KBytes.
  225. */
  226. template<class RealType, int w, unsigned int p, unsigned int q>
  227. class lagged_fibonacci_01_engine
  228. {
  229. public:
  230. typedef RealType result_type;
  231. BOOST_STATIC_CONSTANT(bool, has_fixed_range = false);
  232. BOOST_STATIC_CONSTANT(int, word_size = w);
  233. BOOST_STATIC_CONSTANT(unsigned int, long_lag = p);
  234. BOOST_STATIC_CONSTANT(unsigned int, short_lag = q);
  235. BOOST_STATIC_CONSTANT(boost::uint32_t, default_seed = 331u);
  236. /** Constructs a @c lagged_fibonacci_01 generator and calls @c seed(). */
  237. lagged_fibonacci_01_engine() { seed(); }
  238. /** Constructs a @c lagged_fibonacci_01 generator and calls @c seed(value). */
  239. BOOST_RANDOM_DETAIL_ARITHMETIC_CONSTRUCTOR(lagged_fibonacci_01_engine, uint32_t, value)
  240. { seed(value); }
  241. /** Constructs a @c lagged_fibonacci_01 generator and calls @c seed(gen). */
  242. BOOST_RANDOM_DETAIL_SEED_SEQ_CONSTRUCTOR(lagged_fibonacci_01_engine, SeedSeq, seq)
  243. { seed(seq); }
  244. template<class It> lagged_fibonacci_01_engine(It& first, It last)
  245. { seed(first, last); }
  246. // compiler-generated copy ctor and assignment operator are fine
  247. /** Calls seed(default_seed). */
  248. void seed() { seed(default_seed); }
  249. /**
  250. * Constructs a \minstd_rand0 generator with the constructor parameter
  251. * value and calls seed with it. Distinct seeds in the range
  252. * [1, 2147483647) will produce generators with different states. Other
  253. * seeds will be equivalent to some seed within this range. See
  254. * \linear_congruential_engine for details.
  255. */
  256. BOOST_RANDOM_DETAIL_ARITHMETIC_SEED(lagged_fibonacci_01_engine, boost::uint32_t, value)
  257. {
  258. minstd_rand0 intgen(value);
  259. detail::generator_seed_seq<minstd_rand0> gen(intgen);
  260. seed(gen);
  261. }
  262. /**
  263. * Seeds this @c lagged_fibonacci_01_engine using values produced by
  264. * @c seq.generate.
  265. */
  266. BOOST_RANDOM_DETAIL_SEED_SEQ_SEED(lagged_fibonacci_01_engine, SeedSeq, seq)
  267. {
  268. detail::seed_array_real<w>(seq, x);
  269. i = long_lag;
  270. }
  271. /**
  272. * Seeds this @c lagged_fibonacci_01_engine using values from the
  273. * iterator range [first, last). If there are not enough elements
  274. * in the range, throws @c std::invalid_argument.
  275. */
  276. template<class It>
  277. void seed(It& first, It last)
  278. {
  279. detail::fill_array_real<w>(first, last, x);
  280. i = long_lag;
  281. }
  282. /** Returns the smallest value that the generator can produce. */
  283. static result_type min BOOST_PREVENT_MACRO_SUBSTITUTION () { return result_type(0); }
  284. /** Returns the upper bound of the generators outputs. */
  285. static result_type max BOOST_PREVENT_MACRO_SUBSTITUTION () { return result_type(1); }
  286. /** Returns the next value of the generator. */
  287. result_type operator()()
  288. {
  289. if(i >= long_lag)
  290. fill();
  291. return x[i++];
  292. }
  293. /** Fills a range with random values */
  294. template<class Iter>
  295. void generate(Iter first, Iter last)
  296. { return detail::generate_from_real(*this, first, last); }
  297. /** Advances the state of the generator by @c z. */
  298. void discard(boost::uintmax_t z)
  299. {
  300. for(boost::uintmax_t j = 0; j < z; ++j) {
  301. (*this)();
  302. }
  303. }
  304. /**
  305. * Writes the textual representation of the generator to a @c std::ostream.
  306. */
  307. BOOST_RANDOM_DETAIL_OSTREAM_OPERATOR(os, lagged_fibonacci_01_engine, f)
  308. {
  309. // allow for Koenig lookup
  310. using std::pow;
  311. os << f.i;
  312. std::ios_base::fmtflags oldflags = os.flags(os.dec | os.fixed | os.left);
  313. for(unsigned int j = 0; j < f.long_lag; ++j)
  314. os << ' ' << f.x[j] * f.modulus();
  315. os.flags(oldflags);
  316. return os;
  317. }
  318. /**
  319. * Reads the textual representation of the generator from a @c std::istream.
  320. */
  321. BOOST_RANDOM_DETAIL_ISTREAM_OPERATOR(is, lagged_fibonacci_01_engine, f)
  322. {
  323. is >> f.i;
  324. for(unsigned int j = 0; j < f.long_lag; ++j) {
  325. typename lagged_fibonacci_01_engine::result_type value;
  326. is >> std::ws >> value;
  327. f.x[j] = value / f.modulus();
  328. }
  329. return is;
  330. }
  331. /**
  332. * Returns true if the two generators will produce identical
  333. * sequences of outputs.
  334. */
  335. BOOST_RANDOM_DETAIL_EQUALITY_OPERATOR(lagged_fibonacci_01_engine, x_, y_)
  336. { return x_.i == y_.i && std::equal(x_.x, x_.x+long_lag, y_.x); }
  337. /**
  338. * Returns true if the two generators will produce different
  339. * sequences of outputs.
  340. */
  341. BOOST_RANDOM_DETAIL_INEQUALITY_OPERATOR(lagged_fibonacci_01_engine)
  342. private:
  343. /// \cond show_private
  344. void fill();
  345. static RealType modulus()
  346. {
  347. using std::pow;
  348. return pow(RealType(2), word_size);
  349. }
  350. /// \endcond
  351. unsigned int i;
  352. RealType x[long_lag];
  353. };
  354. #ifndef BOOST_NO_INCLASS_MEMBER_INITIALIZATION
  355. // A definition is required even for integral static constants
  356. template<class RealType, int w, unsigned int p, unsigned int q>
  357. const bool lagged_fibonacci_01_engine<RealType, w, p, q>::has_fixed_range;
  358. template<class RealType, int w, unsigned int p, unsigned int q>
  359. const unsigned int lagged_fibonacci_01_engine<RealType, w, p, q>::long_lag;
  360. template<class RealType, int w, unsigned int p, unsigned int q>
  361. const unsigned int lagged_fibonacci_01_engine<RealType, w, p, q>::short_lag;
  362. template<class RealType, int w, unsigned int p, unsigned int q>
  363. const int lagged_fibonacci_01_engine<RealType,w,p,q>::word_size;
  364. template<class RealType, int w, unsigned int p, unsigned int q>
  365. const boost::uint32_t lagged_fibonacci_01_engine<RealType,w,p,q>::default_seed;
  366. #endif
  367. /// \cond show_private
  368. template<class RealType, int w, unsigned int p, unsigned int q>
  369. void lagged_fibonacci_01_engine<RealType, w, p, q>::fill()
  370. {
  371. // two loops to avoid costly modulo operations
  372. { // extra scope for MSVC brokenness w.r.t. for scope
  373. for(unsigned int j = 0; j < short_lag; ++j) {
  374. RealType t = x[j] + x[j+(long_lag-short_lag)];
  375. if(t >= RealType(1))
  376. t -= RealType(1);
  377. x[j] = t;
  378. }
  379. }
  380. for(unsigned int j = short_lag; j < long_lag; ++j) {
  381. RealType t = x[j] + x[j-short_lag];
  382. if(t >= RealType(1))
  383. t -= RealType(1);
  384. x[j] = t;
  385. }
  386. i = 0;
  387. }
  388. /// \endcond
  389. /// \cond show_deprecated
  390. // provided for backwards compatibility
  391. template<class RealType, int w, unsigned int p, unsigned int q>
  392. class lagged_fibonacci_01 : public lagged_fibonacci_01_engine<RealType, w, p, q>
  393. {
  394. typedef lagged_fibonacci_01_engine<RealType, w, p, q> base_type;
  395. public:
  396. lagged_fibonacci_01() {}
  397. BOOST_RANDOM_DETAIL_ARITHMETIC_CONSTRUCTOR(lagged_fibonacci_01, boost::uint32_t, val)
  398. { this->seed(val); }
  399. BOOST_RANDOM_DETAIL_SEED_SEQ_CONSTRUCTOR(lagged_fibonacci_01, SeedSeq, seq)
  400. { this->seed(seq); }
  401. template<class It>
  402. lagged_fibonacci_01(It& first, It last) : base_type(first, last) {}
  403. };
  404. /// \endcond
  405. namespace detail {
  406. template<class Engine>
  407. struct generator_bits;
  408. template<class RealType, int w, unsigned int p, unsigned int q>
  409. struct generator_bits<lagged_fibonacci_01_engine<RealType, w, p, q> >
  410. {
  411. static std::size_t value() { return w; }
  412. };
  413. template<class RealType, int w, unsigned int p, unsigned int q>
  414. struct generator_bits<lagged_fibonacci_01<RealType, w, p, q> >
  415. {
  416. static std::size_t value() { return w; }
  417. };
  418. }
  419. #ifdef BOOST_RANDOM_DOXYGEN
  420. namespace detail {
  421. /**
  422. * The specializations lagged_fibonacci607 ... lagged_fibonacci44497
  423. * use well tested lags.
  424. *
  425. * See
  426. *
  427. * @blockquote
  428. * "On the Periods of Generalized Fibonacci Recurrences", Richard P. Brent
  429. * Computer Sciences Laboratory Australian National University, December 1992
  430. * @endblockquote
  431. *
  432. * The lags used here can be found in
  433. *
  434. * @blockquote
  435. * "Uniform random number generators for supercomputers", Richard Brent,
  436. * Proc. of Fifth Australian Supercomputer Conference, Melbourne,
  437. * Dec. 1992, pp. 704-706.
  438. * @endblockquote
  439. */
  440. struct lagged_fibonacci_doc {};
  441. }
  442. #endif
  443. /** @copydoc boost::random::detail::lagged_fibonacci_doc */
  444. typedef lagged_fibonacci_01_engine<double, 48, 607, 273> lagged_fibonacci607;
  445. /** @copydoc boost::random::detail::lagged_fibonacci_doc */
  446. typedef lagged_fibonacci_01_engine<double, 48, 1279, 418> lagged_fibonacci1279;
  447. /** @copydoc boost::random::detail::lagged_fibonacci_doc */
  448. typedef lagged_fibonacci_01_engine<double, 48, 2281, 1252> lagged_fibonacci2281;
  449. /** @copydoc boost::random::detail::lagged_fibonacci_doc */
  450. typedef lagged_fibonacci_01_engine<double, 48, 3217, 576> lagged_fibonacci3217;
  451. /** @copydoc boost::random::detail::lagged_fibonacci_doc */
  452. typedef lagged_fibonacci_01_engine<double, 48, 4423, 2098> lagged_fibonacci4423;
  453. /** @copydoc boost::random::detail::lagged_fibonacci_doc */
  454. typedef lagged_fibonacci_01_engine<double, 48, 9689, 5502> lagged_fibonacci9689;
  455. /** @copydoc boost::random::detail::lagged_fibonacci_doc */
  456. typedef lagged_fibonacci_01_engine<double, 48, 19937, 9842> lagged_fibonacci19937;
  457. /** @copydoc boost::random::detail::lagged_fibonacci_doc */
  458. typedef lagged_fibonacci_01_engine<double, 48, 23209, 13470> lagged_fibonacci23209;
  459. /** @copydoc boost::random::detail::lagged_fibonacci_doc */
  460. typedef lagged_fibonacci_01_engine<double, 48, 44497, 21034> lagged_fibonacci44497;
  461. } // namespace random
  462. using random::lagged_fibonacci607;
  463. using random::lagged_fibonacci1279;
  464. using random::lagged_fibonacci2281;
  465. using random::lagged_fibonacci3217;
  466. using random::lagged_fibonacci4423;
  467. using random::lagged_fibonacci9689;
  468. using random::lagged_fibonacci19937;
  469. using random::lagged_fibonacci23209;
  470. using random::lagged_fibonacci44497;
  471. } // namespace boost
  472. #endif // BOOST_RANDOM_LAGGED_FIBONACCI_HPP