boost/random/detail/seed_impl.hpp

   1 /* boost random/detail/seed.hpp header file
   2  *
   3  * Copyright Steven Watanabe 2009
   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: seed_impl.hpp 72951 2011-07-07 04:57:37Z steven_watanabe $
  11  */
  12
  13 #ifndef BOOST_RANDOM_DETAIL_SEED_IMPL_HPP
  14 #define BOOST_RANDOM_DETAIL_SEED_IMPL_HPP
  15
  16 #include <stdexcept>
  17 #include <boost/cstdint.hpp>
  18 #include <boost/config/no_tr1/cmath.hpp>
  19 #include <boost/integer/integer_mask.hpp>
  20 #include <boost/integer/static_log2.hpp>
  21 #include <boost/type_traits/is_signed.hpp>
  22 #include <boost/type_traits/is_integral.hpp>
  23 #include <boost/type_traits/make_unsigned.hpp>
  24 #include <boost/mpl/bool.hpp>
  25 #include <boost/mpl/if.hpp>
  26 #include <boost/mpl/int.hpp>
  27 #include <boost/random/detail/const_mod.hpp>
  28 #include <boost/random/detail/integer_log2.hpp>
  29 #include <boost/random/detail/signed_unsigned_tools.hpp>
  30 #include <boost/random/detail/generator_bits.hpp>
  31
  32 #include <boost/random/detail/disable_warnings.hpp>
  33
  34 namespace boost {
  35 namespace random {
  36 namespace detail {
  37
  38 // finds the seed type of an engine, given its
  39 // result_type.  If the result_type is integral
  40 // the seed type is the same.  If the result_type
  41 // is floating point, the seed type is uint32_t
  42 template<class T>
  43 struct seed_type
  44 {
  45     typedef typename boost::mpl::if_<boost::is_integral<T>,
  46         T,
  47         boost::uint32_t
  48     >::type type;
  49 };
  50
  51 template<int N>
  52 struct const_pow_impl
  53 {
  54     template<class T>
  55     static T call(T arg, int n, T result)
  56     {
  57         return const_pow_impl<N / 2>::call(arg * arg, n / 2,
  58             n%2 == 0? result : result * arg);
  59     }
  60 };
  61
  62 template<>
  63 struct const_pow_impl<0>
  64 {
  65     template<class T>
  66     static T call(T, int, T result)
  67     {
  68         return result;
  69     }
  70 };
  71
  72 // requires N is an upper bound on n
  73 template<int N, class T>
  74 inline T const_pow(T arg, int n) { return const_pow_impl<N>::call(arg, n, T(1)); }
  75
  76 template<class T>
  77 inline T pow2(int n)
  78 {
  79     typedef unsigned int_type;
  80     const int max_bits = std::numeric_limits<int_type>::digits;
  81     T multiplier = T(int_type(1) << (max_bits - 1)) * 2;
  82     return (int_type(1) << (n % max_bits)) *
  83         const_pow<std::numeric_limits<T>::digits / max_bits>(multiplier, n / max_bits);
  84 }
  85
  86 template<class Engine, class Iter>
  87 void generate_from_real(Engine& eng, Iter begin, Iter end)
  88 {
  89     using std::fmod;
  90     typedef typename Engine::result_type RealType;
  91     const int Bits = detail::generator_bits<Engine>::value();
  92     int remaining_bits = 0;
  93     boost::uint_least32_t saved_bits = 0;
  94     RealType multiplier = pow2<RealType>( Bits);
  95     RealType mult32 = RealType(4294967296.0); // 2^32
  96     while(true) {
  97         RealType val = eng() * multiplier;
  98         int available_bits = Bits;
  99         // Make sure the compiler can optimize this out
 100         // if it isn't possible.
 101         if(Bits < 32 && available_bits < 32 - remaining_bits) {
 102             saved_bits |= boost::uint_least32_t(val) << remaining_bits;
 103             remaining_bits += Bits;
 104         } else {
 105             // If Bits < 32, then remaining_bits != 0, since
 106             // if remaining_bits == 0, available_bits < 32 - 0,
 107             // and we won't get here to begin with.
 108             if(Bits < 32 || remaining_bits != 0) {
 109                 boost::uint_least32_t divisor =
 110                     (boost::uint_least32_t(1) << (32 - remaining_bits));
 111                 boost::uint_least32_t extra_bits = boost::uint_least32_t(fmod(val, mult32)) & (divisor - 1);
 112                 val = val / divisor;
 113                 *begin++ = saved_bits | (extra_bits << remaining_bits);
 114                 if(begin == end) return;
 115                 available_bits -= 32 - remaining_bits;
 116                 remaining_bits = 0;
 117             }
 118             // If Bits < 32 we should never enter this loop
 119             if(Bits >= 32) {
 120                 for(; available_bits >= 32; available_bits -= 32) {
 121                     boost::uint_least32_t word = boost::uint_least32_t(fmod(val, mult32));
 122                     val /= mult32;
 123                     *begin++ = word;
 124                     if(begin == end) return;
 125                 }
 126             }
 127             remaining_bits = available_bits;
 128             saved_bits = static_cast<boost::uint_least32_t>(val);
 129         }
 130     }
 131 }
 132
 133 template<class Engine, class Iter>
 134 void generate_from_int(Engine& eng, Iter begin, Iter end)
 135 {
 136     typedef typename Engine::result_type IntType;
 137     typedef typename boost::make_unsigned<IntType>::type unsigned_type;
 138     int remaining_bits = 0;
 139     boost::uint_least32_t saved_bits = 0;
 140     unsigned_type range = boost::random::detail::subtract<IntType>()((eng.max)(), (eng.min)());
 141
 142     int bits =
 143         (range == (std::numeric_limits<unsigned_type>::max)()) ?
 144             std::numeric_limits<unsigned_type>::digits :
 145             detail::integer_log2(range + 1);
 146
 147     {
 148         int discarded_bits = detail::integer_log2(bits);
 149         unsigned_type excess = (range + 1) >> (bits - discarded_bits);
 150         if(excess != 0) {
 151             int extra_bits = detail::integer_log2((excess - 1) ^ excess);
 152             bits = bits - discarded_bits + extra_bits;
 153         }
 154     }
 155
 156     unsigned_type mask = (static_cast<unsigned_type>(2) << (bits - 1)) - 1;
 157     unsigned_type limit = ((range + 1) & ~mask) - 1;
 158
 159     while(true) {
 160         unsigned_type val;
 161         do {
 162             val = boost::random::detail::subtract<IntType>()(eng(), (eng.min)());
 163         } while(limit != range && val > limit);
 164         val &= mask;
 165         int available_bits = bits;
 166         if(available_bits == 32) {
 167             *begin++ = static_cast<boost::uint_least32_t>(val) & 0xFFFFFFFFu;
 168             if(begin == end) return;
 169         } else if(available_bits % 32 == 0) {
 170             for(int i = 0; i < available_bits / 32; ++i) {
 171                 boost::uint_least32_t word = boost::uint_least32_t(val) & 0xFFFFFFFFu;
 172                 int supress_warning = (bits >= 32);
 173                 BOOST_ASSERT(supress_warning == 1);
 174                 val >>= (32 * supress_warning);
 175                 *begin++ = word;
 176                 if(begin == end) return;
 177             }
 178         } else if(bits < 32 && available_bits < 32 - remaining_bits) {
 179             saved_bits |= boost::uint_least32_t(val) << remaining_bits;
 180             remaining_bits += bits;
 181         } else {
 182             if(bits < 32 || remaining_bits != 0) {
 183                 boost::uint_least32_t extra_bits = boost::uint_least32_t(val) & ((boost::uint_least32_t(1) << (32 - remaining_bits)) - 1);
 184                 val >>= 32 - remaining_bits;
 185                 *begin++ = saved_bits | (extra_bits << remaining_bits);
 186                 if(begin == end) return;
 187                 available_bits -= 32 - remaining_bits;
 188                 remaining_bits = 0;
 189             }
 190             if(bits >= 32) {
 191                 for(; available_bits >= 32; available_bits -= 32) {
 192                     boost::uint_least32_t word = boost::uint_least32_t(val) & 0xFFFFFFFFu;
 193                     int supress_warning = (bits >= 32);
 194                     BOOST_ASSERT(supress_warning == 1);
 195                     val >>= (32 * supress_warning);
 196                     *begin++ = word;
 197                     if(begin == end) return;
 198                 }
 199             }
 200             remaining_bits = available_bits;
 201             saved_bits = static_cast<boost::uint_least32_t>(val);
 202         }
 203     }
 204 }
 205
 206 template<class Engine, class Iter>
 207 void generate_impl(Engine& eng, Iter first, Iter last, boost::mpl::true_)
 208 {
 209     return detail::generate_from_int(eng, first, last);
 210 }
 211
 212 template<class Engine, class Iter>
 213 void generate_impl(Engine& eng, Iter first, Iter last, boost::mpl::false_)
 214 {
 215     return detail::generate_from_real(eng, first, last);
 216 }
 217
 218 template<class Engine, class Iter>
 219 void generate(Engine& eng, Iter first, Iter last)
 220 {
 221     return detail::generate_impl(eng, first, last, boost::is_integral<typename Engine::result_type>());
 222 }
 223
 224
 225
 226 template<class IntType, IntType m, class SeedSeq>
 227 IntType seed_one_int(SeedSeq& seq)
 228 {
 229     static const int log = ::boost::mpl::if_c<(m == 0),
 230         ::boost::mpl::int_<(::std::numeric_limits<IntType>::digits)>,
 231         ::boost::static_log2<m> >::type::value;
 232     static const int k =
 233         (log + ((~(static_cast<IntType>(2) << (log - 1)) & m)? 32 : 31)) / 32;
 234     ::boost::uint_least32_t array[log / 32 + 4];
 235     seq.generate(&array[0], &array[0] + k + 3);
 236     IntType s = 0;
 237     for(int j = 0; j < k; ++j) {
 238         IntType digit = const_mod<IntType, m>::apply(IntType(array[j+3]));
 239         IntType mult = IntType(1) << 32*j;
 240         s = const_mod<IntType, m>::mult_add(mult, digit, s);
 241     }
 242     return s;
 243 }
 244
 245 template<class IntType, IntType m, class Iter>
 246 IntType get_one_int(Iter& first, Iter last)
 247 {
 248     static const int log = ::boost::mpl::if_c<(m == 0),
 249         ::boost::mpl::int_<(::std::numeric_limits<IntType>::digits)>,
 250         ::boost::static_log2<m> >::type::value;
 251     static const int k =
 252         (log + ((~(static_cast<IntType>(2) << (log - 1)) & m)? 32 : 31)) / 32;
 253     IntType s = 0;
 254     for(int j = 0; j < k; ++j) {
 255         if(first == last) {
 256             throw ::std::invalid_argument("Not enough elements in call to seed.");
 257         }
 258         IntType digit = const_mod<IntType, m>::apply(IntType(*first++));
 259         IntType mult = IntType(1) << 32*j;
 260         s = const_mod<IntType, m>::mult_add(mult, digit, s);
 261     }
 262     return s;
 263 }
 264
 265 // TODO: work in-place whenever possible
 266 template<int w, std::size_t n, class SeedSeq, class UIntType>
 267 void seed_array_int_impl(SeedSeq& seq, UIntType (&x)[n])
 268 {
 269     boost::uint_least32_t storage[((w+31)/32) * n];
 270     seq.generate(&storage[0], &storage[0] + ((w+31)/32) * n);
 271     for(std::size_t j = 0; j < n; j++) {
 272         UIntType val = 0;
 273         for(std::size_t k = 0; k < (w+31)/32; ++k) {
 274             val += static_cast<UIntType>(storage[(w+31)/32*j + k]) << 32*k;
 275         }
 276         x[j] = val & ::boost::low_bits_mask_t<w>::sig_bits;
 277     }
 278 }
 279
 280 template<int w, std::size_t n, class SeedSeq, class IntType>
 281 inline void seed_array_int_impl(SeedSeq& seq, IntType (&x)[n], boost::mpl::true_)
 282 {
 283     typedef typename boost::make_unsigned<IntType>::type unsigned_array[n];
 284     seed_array_int_impl<w>(seq, reinterpret_cast<unsigned_array&>(x));
 285 }
 286
 287 template<int w, std::size_t n, class SeedSeq, class IntType>
 288 inline void seed_array_int_impl(SeedSeq& seq, IntType (&x)[n], boost::mpl::false_)
 289 {
 290     seed_array_int_impl<w>(seq, x);
 291 }
 292
 293 template<int w, std::size_t n, class SeedSeq, class IntType>
 294 inline void seed_array_int(SeedSeq& seq, IntType (&x)[n])
 295 {
 296     seed_array_int_impl<w>(seq, x, boost::is_signed<IntType>());
 297 }
 298
 299 template<int w, std::size_t n, class Iter, class UIntType>
 300 void fill_array_int_impl(Iter& first, Iter last, UIntType (&x)[n])
 301 {
 302     for(std::size_t j = 0; j < n; j++) {
 303         UIntType val = 0;
 304         for(std::size_t k = 0; k < (w+31)/32; ++k) {
 305             if(first == last) {
 306                 throw std::invalid_argument("Not enough elements in call to seed.");
 307             }
 308             val += static_cast<UIntType>(*first++) << 32*k;
 309         }
 310         x[j] = val & ::boost::low_bits_mask_t<w>::sig_bits;
 311     }
 312 }
 313
 314 template<int w, std::size_t n, class Iter, class IntType>
 315 inline void fill_array_int_impl(Iter& first, Iter last, IntType (&x)[n], boost::mpl::true_)
 316 {
 317     typedef typename boost::make_unsigned<IntType>::type unsigned_array[n];
 318     fill_array_int_impl<w>(first, last, reinterpret_cast<unsigned_array&>(x));
 319 }
 320
 321 template<int w, std::size_t n, class Iter, class IntType>
 322 inline void fill_array_int_impl(Iter& first, Iter last, IntType (&x)[n], boost::mpl::false_)
 323 {
 324     fill_array_int_impl<w>(first, last, x);
 325 }
 326
 327 template<int w, std::size_t n, class Iter, class IntType>
 328 inline void fill_array_int(Iter& first, Iter last, IntType (&x)[n])
 329 {
 330     fill_array_int_impl<w>(first, last, x, boost::is_signed<IntType>());
 331 }
 332
 333 template<int w, std::size_t n, class RealType>
 334 void seed_array_real_impl(const boost::uint_least32_t* storage, RealType (&x)[n])
 335 {
 336     boost::uint_least32_t mask = ~((~boost::uint_least32_t(0)) << (w%32));
 337     RealType two32 = 4294967296.0;
 338     const RealType divisor = RealType(1)/detail::pow2<RealType>(w);
 339     unsigned int j;
 340     for(j = 0; j < n; ++j) {
 341         RealType val = RealType(0);
 342         RealType mult = divisor;
 343         for(int k = 0; k < w/32; ++k) {
 344             val += *storage++ * mult;
 345             mult *= two32;
 346         }
 347         if(mask != 0) {
 348             val += (*storage++ & mask) * mult;
 349         }
 350         BOOST_ASSERT(val >= 0);
 351         BOOST_ASSERT(val < 1);
 352         x[j] = val;
 353     }
 354 }
 355
 356 template<int w, std::size_t n, class SeedSeq, class RealType>
 357 void seed_array_real(SeedSeq& seq, RealType (&x)[n])
 358 {
 359     using std::pow;
 360     boost::uint_least32_t storage[((w+31)/32) * n];
 361     seq.generate(&storage[0], &storage[0] + ((w+31)/32) * n);
 362     seed_array_real_impl<w>(storage, x);
 363 }
 364
 365 template<int w, std::size_t n, class Iter, class RealType>
 366 void fill_array_real(Iter& first, Iter last, RealType (&x)[n])
 367 {
 368     boost::uint_least32_t mask = ~((~boost::uint_least32_t(0)) << (w%32));
 369     RealType two32 = 4294967296.0;
 370     const RealType divisor = RealType(1)/detail::pow2<RealType>(w);
 371     unsigned int j;
 372     for(j = 0; j < n; ++j) {
 373         RealType val = RealType(0);
 374         RealType mult = divisor;
 375         for(int k = 0; k < w/32; ++k, ++first) {
 376             if(first == last) throw std::invalid_argument("Not enough elements in call to seed.");
 377             val += *first * mult;
 378             mult *= two32;
 379         }
 380         if(mask != 0) {
 381             if(first == last) throw std::invalid_argument("Not enough elements in call to seed.");
 382             val += (*first & mask) * mult;
 383             ++first;
 384         }
 385         BOOST_ASSERT(val >= 0);
 386         BOOST_ASSERT(val < 1);
 387         x[j] = val;
 388     }
 389 }
 390
 391 }
 392 }
 393 }
 394
 395 #include <boost/random/detail/enable_warnings.hpp>
 396
 397 #endif