Index: boost/random/uniform_int.hpp =================================================================== --- boost/random/uniform_int.hpp (revision 53231) +++ boost/random/uniform_int.hpp (working copy) @@ -23,13 +23,9 @@ #include #include #include -#include #include #include #include -#ifdef BOOST_NO_LIMITS_COMPILE_TIME_CONSTANTS -#include -#endif namespace boost { @@ -119,46 +115,136 @@ for(;;) { // concatenate several invocations of the base RNG // take extra care to avoid overflows + + // limit == floor((range+1)/(brange+1)) + // Therefore limit*(brange+1) <= range+1 range_type limit; if(range == (std::numeric_limits::max)()) { limit = range/(range_type(brange)+1); - if(range % range_type(brange)+1 == range_type(brange)) + if(range % (range_type(brange)+1) == range_type(brange)) ++limit; } else { limit = (range+1)/(range_type(brange)+1); } + // We consider "result" as expressed to base (brange+1): // For every power of (brange+1), we determine a random factor range_type result = range_type(0); range_type mult = range_type(1); + + // loop invariants: + // result < mult + // mult <= range while(mult <= limit) { + // Postcondition: result <= range, thus no overflow + // + // limit*(brange+1)<=range+1 def. of limit (1) + // eng()-bmin<=brange eng() post. (2) + // and mult<=limit. loop condition (3) + // Therefore mult*(eng()-bmin+1)<=range+1 by (1),(2),(3) (4) + // Therefore mult*(eng()-bmin)+mult<=range+1 rearranging (4) (5) + // result()(eng(), bmin) * mult; + + // equivalent to (mult * (brange+1)) == range+1, but avoids overflow. + if(mult * range_type(brange) == range - mult + 1) { + // The destination range is an integer power of + // the generator's range. + return(result); + } + + // Postcondition: mult <= range + // + // limit*(brange+1)<=range+1 def. of limit (1) + // mult<=limit loop condition (2) + // Therefore mult*(brange+1)<=range+1 by (1), (2) (3) + // mult*(brange+1)!=range+1 preceding if (4) + // Therefore mult*(brange+1) limit loop condition (1) + // Suppose range/mult >= brange+1 Assumption (2) + // range >= mult*(brange+1) by (2) (3) + // range+1 > mult*(brange+1) by (3) (4) + // range+1 > (limit+1)*(brange+1) by (1), (4) (5) + // (range+1)/(brange+1) > limit+1 by (5) (6) + // limit < floor((range+1)/(brange+1)) by (6) (7) + // limit==floor((range+1)/(brange+1)) def. of limit (8) + // not (2) reductio (9) + // + // loop postcondition: (range/mult)*mult+(mult-1) >= range + // + // (range/mult)*mult + range%mult == range identity (1) + // range%mult < mult def. of % (2) + // (range/mult)*mult+mult > range by (1), (2) (3) + // (range/mult)*mult+(mult-1) >= range by (3) (4) + // + // Note that the maximum value of result at this point is (mult-1), + // so after this final step, we generate numbers that can be + // at least as large as range. We have to really careful to avoid + // overflow in this final addition and in the rejection. Anything + // that overflows is larger than range and can thus be rejected. + // range/mult < brange+1 -> no endless loop - result += uniform_int(0, range/mult)(eng) * mult; - if(result <= range) - return random::detail::add()(result, min_value); + range_type result_increment = uniform_int(0, range/mult)(eng); + if((std::numeric_limits::max)() / mult < result_increment) { + // The multiplcation would overflow. Reject immediately. + continue; + } + result_increment *= mult; + // unsigned integers are guaranteed to wrap on overflow. + result += result_increment; + if(result < result_increment) { + // The addition overflowed. Reject. + continue; + } + if(result > range) { + // Too big. Reject. + continue; + } + return random::detail::add()(result, min_value); } } else { // brange > range - if(brange / range > 4 /* quantization_cutoff */ ) { - // the new range is vastly smaller than the source range, - // so quantization effects are not relevant - return boost::uniform_smallint(min_value, max_value)(eng); + base_unsigned bucket_size; + // it's safe to add 1 to range, as long as we cast it first, + // because we know that it is less than brange. However, + // we do need to be careful not to cause overflow by adding 1 + // to brange. + if(brange == (std::numeric_limits::max)()) { + bucket_size = brange / (static_cast(range)+1); + if(brange % (static_cast(range)+1) == static_cast(range)) { + ++bucket_size; + } } else { - // use rejection method to handle cases like 0..5 -> 0..4 - for(;;) { - base_unsigned result = - random::detail::subtract()(eng(), bmin); - // result and range are non-negative, and result is possibly larger - // than range, so the cast is safe - if(result <= static_cast(range)) - return random::detail::add()(result, min_value); - } + bucket_size = (brange+1) / (static_cast(range)+1); } + for(;;) { + base_unsigned result = + random::detail::subtract()(eng(), bmin); + result /= bucket_size; + // result and range are non-negative, and result is possibly larger + // than range, so the cast is safe + if(result <= static_cast(range)) + return random::detail::add()(result, min_value); + } } } Index: libs/random/random_test.cpp =================================================================== --- libs/random/random_test.cpp (revision 53231) +++ libs/random/random_test.cpp (working copy) @@ -64,7 +64,8 @@ for(int k = 0; k < range; k++) sum += bucket[k]; double avg = static_cast(sum)/range; - double threshold = 2*avg/std::sqrt(static_cast(iter)); + double p = 1 / static_cast(range); + double threshold = 2*std::sqrt(static_cast(iter)*p*(1-p)); for(int i = 0; i < range; i++) { if(std::fabs(bucket[i] - avg) > threshold) { // 95% confidence interval @@ -100,11 +101,37 @@ // small range => larger range level_two uint05(uint12, int_gen(-3, 2)); check_uniform_int(uint05, 100000); + + // small range => larger range + level_two uint099(uint12, int_gen(0, 99)); + check_uniform_int(uint099, 100000); // larger => small range, rejection case typedef boost::variate_generator level_three; level_three uint1_4(uint05, int_gen(1, 4)); check_uniform_int(uint1_4, 100000); + + typedef boost::uniform_int int8_gen; + typedef boost::variate_generator gen8_t; + + gen8_t gen8_03(gen, int8_gen(0, 3)); + + // use the full range of the type, where the destination + // range is a power of the source range + typedef boost::variate_generator uniform_uint8; + uniform_uint8 uint8_0255(gen8_03, int8_gen(0, 255)); + check_uniform_int(uint8_0255, 100000); + + // use the full range, but a generator whose range is not + // a root of the destination range. + gen8_t gen8_02(gen, int8_gen(0, 2)); + uniform_uint8 uint8_0255_2(gen8_02, int8_gen(0, 255)); + check_uniform_int(uint8_0255_2, 100000); + + // expand the range to a larger type. + typedef boost::variate_generator uniform_uint_from8; + uniform_uint_from8 uint0300(gen8_03, int_gen(0, 300)); + check_uniform_int(uint0300, 100000); } #if defined(BOOST_MSVC) && _MSC_VER < 1300 @@ -140,10 +167,13 @@ class ruetti_gen { public: + ruetti_gen() : state((max)() - 1) {} typedef boost::uint64_t result_type; result_type min BOOST_PREVENT_MACRO_SUBSTITUTION () const { return 0; } result_type max BOOST_PREVENT_MACRO_SUBSTITUTION () const { return std::numeric_limits::max BOOST_PREVENT_MACRO_SUBSTITUTION (); } - result_type operator()() { return (max)()-1; } + result_type operator()() { return state--; } +private: + result_type state; }; void test_overflow_range()