uniform_int_distribution.hpp 16 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419
  1. /* boost random/uniform_int_distribution.hpp header file
  2. *
  3. * Copyright Jens Maurer 2000-2001
  4. * Copyright Steven Watanabe 2011
  5. * Distributed under the Boost Software License, Version 1.0. (See
  6. * accompanying file LICENSE_1_0.txt or copy at
  7. * http://www.boost.org/LICENSE_1_0.txt)
  8. *
  9. * See http://www.boost.org for most recent version including documentation.
  10. *
  11. * $Id$
  12. *
  13. * Revision history
  14. * 2001-04-08 added min<max assertion (N. Becker)
  15. * 2001-02-18 moved to individual header files
  16. */
  17. #ifndef BOOST_RANDOM_UNIFORM_INT_DISTRIBUTION_HPP
  18. #define BOOST_RANDOM_UNIFORM_INT_DISTRIBUTION_HPP
  19. #include <iosfwd>
  20. #include <ios>
  21. #include <istream>
  22. #include <boost/config.hpp>
  23. #include <boost/limits.hpp>
  24. #include <boost/assert.hpp>
  25. #include <boost/random/detail/config.hpp>
  26. #include <boost/random/detail/operators.hpp>
  27. #include <boost/random/detail/uniform_int_float.hpp>
  28. #include <boost/random/detail/signed_unsigned_tools.hpp>
  29. #include <boost/random/traits.hpp>
  30. #include <boost/type_traits/integral_constant.hpp>
  31. #ifdef BOOST_NO_CXX11_EXPLICIT_CONVERSION_OPERATORS
  32. #include <boost/type_traits/conditional.hpp>
  33. #endif
  34. namespace boost {
  35. namespace random {
  36. namespace detail {
  37. #ifdef BOOST_MSVC
  38. #pragma warning(push)
  39. // disable division by zero warning, since we can't
  40. // actually divide by zero.
  41. #pragma warning(disable:4723)
  42. #endif
  43. template<class Engine, class T>
  44. T generate_uniform_int(
  45. Engine& eng, T min_value, T max_value,
  46. boost::true_type /** is_integral<Engine::result_type> */)
  47. {
  48. typedef T result_type;
  49. typedef typename boost::random::traits::make_unsigned_or_unbounded<T>::type range_type;
  50. typedef typename Engine::result_type base_result;
  51. // ranges are always unsigned or unbounded
  52. typedef typename boost::random::traits::make_unsigned_or_unbounded<base_result>::type base_unsigned;
  53. const range_type range = random::detail::subtract<result_type>()(max_value, min_value);
  54. const base_result bmin = (eng.min)();
  55. const base_unsigned brange =
  56. random::detail::subtract<base_result>()((eng.max)(), (eng.min)());
  57. if(range == 0) {
  58. return min_value;
  59. } else if(brange == range) {
  60. // this will probably never happen in real life
  61. // basically nothing to do; just take care we don't overflow / underflow
  62. base_unsigned v = random::detail::subtract<base_result>()(eng(), bmin);
  63. return random::detail::add<base_unsigned, result_type>()(v, min_value);
  64. } else if(brange < range) {
  65. // use rejection method to handle things like 0..3 --> 0..4
  66. for(;;) {
  67. // concatenate several invocations of the base RNG
  68. // take extra care to avoid overflows
  69. // limit == floor((range+1)/(brange+1))
  70. // Therefore limit*(brange+1) <= range+1
  71. range_type limit;
  72. if(range == (std::numeric_limits<range_type>::max)()) {
  73. limit = range/(range_type(brange)+1);
  74. if(range % (range_type(brange)+1) == range_type(brange))
  75. ++limit;
  76. } else {
  77. limit = (range+1)/(range_type(brange)+1);
  78. }
  79. // We consider "result" as expressed to base (brange+1):
  80. // For every power of (brange+1), we determine a random factor
  81. range_type result = range_type(0);
  82. range_type mult = range_type(1);
  83. // loop invariants:
  84. // result < mult
  85. // mult <= range
  86. while(mult <= limit) {
  87. // Postcondition: result <= range, thus no overflow
  88. //
  89. // limit*(brange+1)<=range+1 def. of limit (1)
  90. // eng()-bmin<=brange eng() post. (2)
  91. // and mult<=limit. loop condition (3)
  92. // Therefore mult*(eng()-bmin+1)<=range+1 by (1),(2),(3) (4)
  93. // Therefore mult*(eng()-bmin)+mult<=range+1 rearranging (4) (5)
  94. // result<mult loop invariant (6)
  95. // Therefore result+mult*(eng()-bmin)<range+1 by (5), (6) (7)
  96. //
  97. // Postcondition: result < mult*(brange+1)
  98. //
  99. // result<mult loop invariant (1)
  100. // eng()-bmin<=brange eng() post. (2)
  101. // Therefore result+mult*(eng()-bmin) <
  102. // mult+mult*(eng()-bmin) by (1) (3)
  103. // Therefore result+(eng()-bmin)*mult <
  104. // mult+mult*brange by (2), (3) (4)
  105. // Therefore result+(eng()-bmin)*mult <
  106. // mult*(brange+1) by (4)
  107. result += static_cast<range_type>(static_cast<range_type>(random::detail::subtract<base_result>()(eng(), bmin)) * mult);
  108. // equivalent to (mult * (brange+1)) == range+1, but avoids overflow.
  109. if(mult * range_type(brange) == range - mult + 1) {
  110. // The destination range is an integer power of
  111. // the generator's range.
  112. return(result);
  113. }
  114. // Postcondition: mult <= range
  115. //
  116. // limit*(brange+1)<=range+1 def. of limit (1)
  117. // mult<=limit loop condition (2)
  118. // Therefore mult*(brange+1)<=range+1 by (1), (2) (3)
  119. // mult*(brange+1)!=range+1 preceding if (4)
  120. // Therefore mult*(brange+1)<range+1 by (3), (4) (5)
  121. //
  122. // Postcondition: result < mult
  123. //
  124. // See the second postcondition on the change to result.
  125. mult *= range_type(brange)+range_type(1);
  126. }
  127. // loop postcondition: range/mult < brange+1
  128. //
  129. // mult > limit loop condition (1)
  130. // Suppose range/mult >= brange+1 Assumption (2)
  131. // range >= mult*(brange+1) by (2) (3)
  132. // range+1 > mult*(brange+1) by (3) (4)
  133. // range+1 > (limit+1)*(brange+1) by (1), (4) (5)
  134. // (range+1)/(brange+1) > limit+1 by (5) (6)
  135. // limit < floor((range+1)/(brange+1)) by (6) (7)
  136. // limit==floor((range+1)/(brange+1)) def. of limit (8)
  137. // not (2) reductio (9)
  138. //
  139. // loop postcondition: (range/mult)*mult+(mult-1) >= range
  140. //
  141. // (range/mult)*mult + range%mult == range identity (1)
  142. // range%mult < mult def. of % (2)
  143. // (range/mult)*mult+mult > range by (1), (2) (3)
  144. // (range/mult)*mult+(mult-1) >= range by (3) (4)
  145. //
  146. // Note that the maximum value of result at this point is (mult-1),
  147. // so after this final step, we generate numbers that can be
  148. // at least as large as range. We have to really careful to avoid
  149. // overflow in this final addition and in the rejection. Anything
  150. // that overflows is larger than range and can thus be rejected.
  151. // range/mult < brange+1 -> no endless loop
  152. range_type result_increment =
  153. generate_uniform_int(
  154. eng,
  155. static_cast<range_type>(0),
  156. static_cast<range_type>(range/mult),
  157. boost::true_type());
  158. if(std::numeric_limits<range_type>::is_bounded && ((std::numeric_limits<range_type>::max)() / mult < result_increment)) {
  159. // The multiplcation would overflow. Reject immediately.
  160. continue;
  161. }
  162. result_increment *= mult;
  163. // unsigned integers are guaranteed to wrap on overflow.
  164. result += result_increment;
  165. if(result < result_increment) {
  166. // The addition overflowed. Reject.
  167. continue;
  168. }
  169. if(result > range) {
  170. // Too big. Reject.
  171. continue;
  172. }
  173. return random::detail::add<range_type, result_type>()(result, min_value);
  174. }
  175. } else { // brange > range
  176. #ifdef BOOST_NO_CXX11_EXPLICIT_CONVERSION_OPERATORS
  177. typedef typename conditional<
  178. std::numeric_limits<range_type>::is_specialized && std::numeric_limits<base_unsigned>::is_specialized
  179. && (std::numeric_limits<range_type>::digits >= std::numeric_limits<base_unsigned>::digits),
  180. range_type, base_unsigned>::type mixed_range_type;
  181. #else
  182. typedef base_unsigned mixed_range_type;
  183. #endif
  184. mixed_range_type bucket_size;
  185. // it's safe to add 1 to range, as long as we cast it first,
  186. // because we know that it is less than brange. However,
  187. // we do need to be careful not to cause overflow by adding 1
  188. // to brange. We use mixed_range_type throughout for mixed
  189. // arithmetic between base_unsigned and range_type - in the case
  190. // that range_type has more bits than base_unsigned it is always
  191. // safe to use range_type for this albeit it may be more effient
  192. // to use base_unsigned. The latter is a narrowing conversion though
  193. // which may be disallowed if range_type is a multiprecision type
  194. // and there are no explicit converison operators.
  195. if(brange == (std::numeric_limits<base_unsigned>::max)()) {
  196. bucket_size = static_cast<mixed_range_type>(brange) / (static_cast<mixed_range_type>(range)+1);
  197. if(static_cast<mixed_range_type>(brange) % (static_cast<mixed_range_type>(range)+1) == static_cast<mixed_range_type>(range)) {
  198. ++bucket_size;
  199. }
  200. } else {
  201. bucket_size = static_cast<mixed_range_type>(brange + 1) / (static_cast<mixed_range_type>(range)+1);
  202. }
  203. for(;;) {
  204. mixed_range_type result =
  205. random::detail::subtract<base_result>()(eng(), bmin);
  206. result /= bucket_size;
  207. // result and range are non-negative, and result is possibly larger
  208. // than range, so the cast is safe
  209. if(result <= static_cast<mixed_range_type>(range))
  210. return random::detail::add<mixed_range_type, result_type>()(result, min_value);
  211. }
  212. }
  213. }
  214. #ifdef BOOST_MSVC
  215. #pragma warning(pop)
  216. #endif
  217. template<class Engine, class T>
  218. inline T generate_uniform_int(
  219. Engine& eng, T min_value, T max_value,
  220. boost::false_type /** is_integral<Engine::result_type> */)
  221. {
  222. uniform_int_float<Engine> wrapper(eng);
  223. return generate_uniform_int(wrapper, min_value, max_value, boost::true_type());
  224. }
  225. template<class Engine, class T>
  226. inline T generate_uniform_int(Engine& eng, T min_value, T max_value)
  227. {
  228. typedef typename Engine::result_type base_result;
  229. return generate_uniform_int(eng, min_value, max_value,
  230. boost::random::traits::is_integral<base_result>());
  231. }
  232. }
  233. /**
  234. * The class template uniform_int_distribution models a \random_distribution.
  235. * On each invocation, it returns a random integer value uniformly
  236. * distributed in the set of integers {min, min+1, min+2, ..., max}.
  237. *
  238. * The template parameter IntType shall denote an integer-like value type.
  239. */
  240. template<class IntType = int>
  241. class uniform_int_distribution
  242. {
  243. public:
  244. typedef IntType input_type;
  245. typedef IntType result_type;
  246. class param_type
  247. {
  248. public:
  249. typedef uniform_int_distribution distribution_type;
  250. /**
  251. * Constructs the parameters of a uniform_int_distribution.
  252. *
  253. * Requires min <= max
  254. */
  255. explicit param_type(
  256. IntType min_arg = 0,
  257. IntType max_arg = (std::numeric_limits<IntType>::max)())
  258. : _min(min_arg), _max(max_arg)
  259. {
  260. BOOST_ASSERT(_min <= _max);
  261. }
  262. /** Returns the minimum value of the distribution. */
  263. IntType a() const { return _min; }
  264. /** Returns the maximum value of the distribution. */
  265. IntType b() const { return _max; }
  266. /** Writes the parameters to a @c std::ostream. */
  267. BOOST_RANDOM_DETAIL_OSTREAM_OPERATOR(os, param_type, parm)
  268. {
  269. os << parm._min << " " << parm._max;
  270. return os;
  271. }
  272. /** Reads the parameters from a @c std::istream. */
  273. BOOST_RANDOM_DETAIL_ISTREAM_OPERATOR(is, param_type, parm)
  274. {
  275. IntType min_in, max_in;
  276. if(is >> min_in >> std::ws >> max_in) {
  277. if(min_in <= max_in) {
  278. parm._min = min_in;
  279. parm._max = max_in;
  280. } else {
  281. is.setstate(std::ios_base::failbit);
  282. }
  283. }
  284. return is;
  285. }
  286. /** Returns true if the two sets of parameters are equal. */
  287. BOOST_RANDOM_DETAIL_EQUALITY_OPERATOR(param_type, lhs, rhs)
  288. { return lhs._min == rhs._min && lhs._max == rhs._max; }
  289. /** Returns true if the two sets of parameters are different. */
  290. BOOST_RANDOM_DETAIL_INEQUALITY_OPERATOR(param_type)
  291. private:
  292. IntType _min;
  293. IntType _max;
  294. };
  295. /**
  296. * Constructs a uniform_int_distribution. @c min and @c max are
  297. * the parameters of the distribution.
  298. *
  299. * Requires: min <= max
  300. */
  301. explicit uniform_int_distribution(
  302. IntType min_arg = 0,
  303. IntType max_arg = (std::numeric_limits<IntType>::max)())
  304. : _min(min_arg), _max(max_arg)
  305. {
  306. BOOST_ASSERT(min_arg <= max_arg);
  307. }
  308. /** Constructs a uniform_int_distribution from its parameters. */
  309. explicit uniform_int_distribution(const param_type& parm)
  310. : _min(parm.a()), _max(parm.b()) {}
  311. /** Returns the minimum value of the distribution */
  312. IntType min BOOST_PREVENT_MACRO_SUBSTITUTION () const { return _min; }
  313. /** Returns the maximum value of the distribution */
  314. IntType max BOOST_PREVENT_MACRO_SUBSTITUTION () const { return _max; }
  315. /** Returns the minimum value of the distribution */
  316. IntType a() const { return _min; }
  317. /** Returns the maximum value of the distribution */
  318. IntType b() const { return _max; }
  319. /** Returns the parameters of the distribution. */
  320. param_type param() const { return param_type(_min, _max); }
  321. /** Sets the parameters of the distribution. */
  322. void param(const param_type& parm)
  323. {
  324. _min = parm.a();
  325. _max = parm.b();
  326. }
  327. /**
  328. * Effects: Subsequent uses of the distribution do not depend
  329. * on values produced by any engine prior to invoking reset.
  330. */
  331. void reset() { }
  332. /** Returns an integer uniformly distributed in the range [min, max]. */
  333. template<class Engine>
  334. result_type operator()(Engine& eng) const
  335. { return detail::generate_uniform_int(eng, _min, _max); }
  336. /**
  337. * Returns an integer uniformly distributed in the range
  338. * [param.a(), param.b()].
  339. */
  340. template<class Engine>
  341. result_type operator()(Engine& eng, const param_type& parm) const
  342. { return detail::generate_uniform_int(eng, parm.a(), parm.b()); }
  343. /** Writes the distribution to a @c std::ostream. */
  344. BOOST_RANDOM_DETAIL_OSTREAM_OPERATOR(os, uniform_int_distribution, ud)
  345. {
  346. os << ud.param();
  347. return os;
  348. }
  349. /** Reads the distribution from a @c std::istream. */
  350. BOOST_RANDOM_DETAIL_ISTREAM_OPERATOR(is, uniform_int_distribution, ud)
  351. {
  352. param_type parm;
  353. if(is >> parm) {
  354. ud.param(parm);
  355. }
  356. return is;
  357. }
  358. /**
  359. * Returns true if the two distributions will produce identical sequences
  360. * of values given equal generators.
  361. */
  362. BOOST_RANDOM_DETAIL_EQUALITY_OPERATOR(uniform_int_distribution, lhs, rhs)
  363. { return lhs._min == rhs._min && lhs._max == rhs._max; }
  364. /**
  365. * Returns true if the two distributions may produce different sequences
  366. * of values given equal generators.
  367. */
  368. BOOST_RANDOM_DETAIL_INEQUALITY_OPERATOR(uniform_int_distribution)
  369. private:
  370. IntType _min;
  371. IntType _max;
  372. };
  373. } // namespace random
  374. } // namespace boost
  375. #endif // BOOST_RANDOM_UNIFORM_INT_HPP