ceph/src/boost/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$
  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/throw_exception.hpp>
  19 #include <boost/config/no_tr1/cmath.hpp>
  20 #include <boost/integer/integer_mask.hpp>
  21 #include <boost/integer/static_log2.hpp>
  22 #include <boost/random/traits.hpp>
  23 #include <boost/random/detail/const_mod.hpp>
  24 #include <boost/random/detail/integer_log2.hpp>
  25 #include <boost/random/detail/signed_unsigned_tools.hpp>
  26 #include <boost/random/detail/generator_bits.hpp>
  27 #include <boost/type_traits/conditional.hpp>
  28 #include <boost/type_traits/integral_constant.hpp>
  29
  30 #include <boost/random/detail/disable_warnings.hpp>
  31
  32 namespace boost {
  33 namespace random {
  34 namespace detail {
  35
  36 // finds the seed type of an engine, given its
  37 // result_type.  If the result_type is integral
  38 // the seed type is the same.  If the result_type
  39 // is floating point, the seed type is uint32_t
  40 template<class T>
  41 struct seed_type
  42 {
  43     typedef typename boost::conditional<boost::is_integral<T>::value,
  44         T,
  45         boost::uint32_t
  46     >::type type;
  47 };
  48
  49 template<int N>
  50 struct const_pow_impl
  51 {
  52     template<class T>
  53     static T call(T arg, int n, T result)
  54     {
  55         return const_pow_impl<N / 2>::call(T(arg * arg), n / 2,
  56             n%2 == 0? result : T(result * arg));
  57     }
  58 };
  59
  60 template<>
  61 struct const_pow_impl<0>
  62 {
  63     template<class T>
  64     static T call(T, int, T result)
  65     {
  66         return result;
  67     }
  68 };
  69
  70 // requires N is an upper bound on n
  71 template<int N, class T>
  72 inline T const_pow(T arg, int n) { return const_pow_impl<N>::call(arg, n, T(1)); }
  73
  74 template<class T>
  75 inline T pow2(int n)
  76 {
  77     typedef unsigned int_type;
  78     const int max_bits = std::numeric_limits<int_type>::digits;
  79     T multiplier = T(int_type(1) << (max_bits - 1)) * 2;
  80     return (int_type(1) << (n % max_bits)) *
  81         const_pow<std::numeric_limits<T>::digits / max_bits>(multiplier, n / max_bits);
  82 }
  83
  84 template<class Engine, class Iter>
  85 void generate_from_real(Engine& eng, Iter begin, Iter end)
  86 {
  87     using std::fmod;
  88     typedef typename Engine::result_type RealType;
  89     const int Bits = detail::generator_bits<Engine>::value();
  90     int remaining_bits = 0;
  91     boost::uint_least32_t saved_bits = 0;
  92     RealType multiplier = pow2<RealType>( Bits);
  93     RealType mult32 = RealType(4294967296.0); // 2^32
  94     while(true) {
  95         RealType val = eng() * multiplier;
  96         int available_bits = Bits;
  97         // Make sure the compiler can optimize this out
  98         // if it isn't possible.
  99         if(Bits < 32 && available_bits < 32 - remaining_bits) {
 100             saved_bits |= boost::uint_least32_t(val) << remaining_bits;
 101             remaining_bits += Bits;
 102         } else {
 103             // If Bits < 32, then remaining_bits != 0, since
 104             // if remaining_bits == 0, available_bits < 32 - 0,
 105             // and we won't get here to begin with.
 106             if(Bits < 32 || remaining_bits != 0) {
 107                 boost::uint_least32_t divisor =
 108                     (boost::uint_least32_t(1) << (32 - remaining_bits));
 109                 boost::uint_least32_t extra_bits = boost::uint_least32_t(fmod(val, mult32)) & (divisor - 1);
 110                 val = val / divisor;
 111                 *begin++ = saved_bits | (extra_bits << remaining_bits);
 112                 if(begin == end) return;
 113                 available_bits -= 32 - remaining_bits;
 114                 remaining_bits = 0;
 115             }
 116             // If Bits < 32 we should never enter this loop
 117             if(Bits >= 32) {
 118                 for(; available_bits >= 32; available_bits -= 32) {
 119                     boost::uint_least32_t word = boost::uint_least32_t(fmod(val, mult32));
 120                     val /= mult32;
 121                     *begin++ = word;
 122                     if(begin == end) return;
 123                 }
 124             }
 125             remaining_bits = available_bits;
 126             saved_bits = static_cast<boost::uint_least32_t>(val);
 127         }
 128     }
 129 }
 130
 131 template<class Engine, class Iter>
 132 void generate_from_int(Engine& eng, Iter begin, Iter end)
 133 {
 134     typedef typename Engine::result_type IntType;
 135     typedef typename boost::random::traits::make_unsigned<IntType>::type unsigned_type;
 136     int remaining_bits = 0;
 137     boost::uint_least32_t saved_bits = 0;
 138     unsigned_type range = boost::random::detail::subtract<IntType>()((eng.max)(), (eng.min)());
 139
 140     int bits =
 141         (range == (std::numeric_limits<unsigned_type>::max)()) ?
 142             std::numeric_limits<unsigned_type>::digits :
 143             detail::integer_log2(range + 1);
 144
 145     {
 146         int discarded_bits = detail::integer_log2(bits);
 147         unsigned_type excess = (range + 1) >> (bits - discarded_bits);
 148         if(excess != 0) {
 149             int extra_bits = detail::integer_log2((excess - 1) ^ excess);
 150             bits = bits - discarded_bits + extra_bits;
 151         }
 152     }
 153
 154     unsigned_type mask = (static_cast<unsigned_type>(2) << (bits - 1)) - 1;
 155     unsigned_type limit = ((range + 1) & ~mask) - 1;
 156
 157     while(true) {
 158         unsigned_type val;
 159         do {
 160             val = boost::random::detail::subtract<IntType>()(eng(), (eng.min)());
 161         } while(limit != range && val > limit);
 162         val &= mask;
 163         int available_bits = bits;
 164         if(available_bits == 32) {
 165             *begin++ = static_cast<boost::uint_least32_t>(val) & 0xFFFFFFFFu;
 166             if(begin == end) return;
 167         } else if(available_bits % 32 == 0) {
 168             for(int i = 0; i < available_bits / 32; ++i) {
 169                 boost::uint_least32_t word = boost::uint_least32_t(val) & 0xFFFFFFFFu;
 170                 int suppress_warning = (bits >= 32);
 171                 BOOST_ASSERT(suppress_warning == 1);
 172                 val >>= (32 * suppress_warning);
 173                 *begin++ = word;
 174                 if(begin == end) return;
 175             }
 176         } else if(bits < 32 && available_bits < 32 - remaining_bits) {
 177             saved_bits |= boost::uint_least32_t(val) << remaining_bits;
 178             remaining_bits += bits;
 179         } else {
 180             if(bits < 32 || remaining_bits != 0) {
 181                 boost::uint_least32_t extra_bits = boost::uint_least32_t(val) & ((boost::uint_least32_t(1) << (32 - remaining_bits)) - 1);
 182                 val >>= 32 - remaining_bits;
 183                 *begin++ = saved_bits | (extra_bits << remaining_bits);
 184                 if(begin == end) return;
 185                 available_bits -= 32 - remaining_bits;
 186                 remaining_bits = 0;
 187             }
 188             if(bits >= 32) {
 189                 for(; available_bits >= 32; available_bits -= 32) {
 190                     boost::uint_least32_t word = boost::uint_least32_t(val) & 0xFFFFFFFFu;
 191                     int suppress_warning = (bits >= 32);
 192                     BOOST_ASSERT(suppress_warning == 1);
 193                     val >>= (32 * suppress_warning);
 194                     *begin++ = word;
 195                     if(begin == end) return;
 196                 }
 197             }
 198             remaining_bits = available_bits;
 199             saved_bits = static_cast<boost::uint_least32_t>(val);
 200         }
 201     }
 202 }
 203
 204 template<class Engine, class Iter>
 205 void generate_impl(Engine& eng, Iter first, Iter last, boost::true_type)
 206 {
 207     return detail::generate_from_int(eng, first, last);
 208 }
 209
 210 template<class Engine, class Iter>
 211 void generate_impl(Engine& eng, Iter first, Iter last, boost::false_type)
 212 {
 213     return detail::generate_from_real(eng, first, last);
 214 }
 215
 216 template<class Engine, class Iter>
 217 void generate(Engine& eng, Iter first, Iter last)
 218 {
 219     return detail::generate_impl(eng, first, last, boost::random::traits::is_integral<typename Engine::result_type>());
 220 }
 221
 222
 223
 224 template<class IntType, IntType m, class SeedSeq>
 225 IntType seed_one_int(SeedSeq& seq)
 226 {
 227     static const int log = ::boost::conditional<(m == 0),
 228         ::boost::integral_constant<int, (::std::numeric_limits<IntType>::digits)>,
 229         ::boost::static_log2<m> >::type::value;
 230     static const int k =
 231         (log + ((~(static_cast<IntType>(2) << (log - 1)) & m)? 32 : 31)) / 32;
 232     ::boost::uint_least32_t array[log / 32 + 4];
 233     seq.generate(&array[0], &array[0] + k + 3);
 234     IntType s = 0;
 235     for(int j = 0; j < k; ++j) {
 236         IntType digit = const_mod<IntType, m>::apply(IntType(array[j+3]));
 237         IntType mult = IntType(1) << 32*j;
 238         s = const_mod<IntType, m>::mult_add(mult, digit, s);
 239     }
 240     return s;
 241 }
 242
 243 template<class IntType, IntType m, class Iter>
 244 IntType get_one_int(Iter& first, Iter last)
 245 {
 246     static const int log = ::boost::conditional<(m == 0),
 247         ::boost::integral_constant<int, (::std::numeric_limits<IntType>::digits)>,
 248         ::boost::static_log2<m> >::type::value;
 249     static const int k =
 250         (log + ((~(static_cast<IntType>(2) << (log - 1)) & m)? 32 : 31)) / 32;
 251     IntType s = 0;
 252     for(int j = 0; j < k; ++j) {
 253         if(first == last) {
 254             boost::throw_exception(::std::invalid_argument("Not enough elements in call to seed."));
 255         }
 256         IntType digit = const_mod<IntType, m>::apply(IntType(*first++));
 257         IntType mult = IntType(1) << 32*j;
 258         s = const_mod<IntType, m>::mult_add(mult, digit, s);
 259     }
 260     return s;
 261 }
 262
 263 // TODO: work in-place whenever possible
 264 template<int w, std::size_t n, class SeedSeq, class UIntType>
 265 void seed_array_int_impl(SeedSeq& seq, UIntType (&x)[n])
 266 {
 267     boost::uint_least32_t storage[((w+31)/32) * n];
 268     seq.generate(&storage[0], &storage[0] + ((w+31)/32) * n);
 269     for(std::size_t j = 0; j < n; j++) {
 270         UIntType val = 0;
 271         for(std::size_t k = 0; k < (w+31)/32; ++k) {
 272             val += static_cast<UIntType>(storage[(w+31)/32*j + k]) << 32*k;
 273         }
 274         x[j] = val & ::boost::low_bits_mask_t<w>::sig_bits;
 275     }
 276 }
 277
 278 template<int w, std::size_t n, class SeedSeq, class IntType>
 279 inline void seed_array_int_impl(SeedSeq& seq, IntType (&x)[n], boost::true_type)
 280 {
 281     BOOST_STATIC_ASSERT_MSG(boost::is_integral<IntType>::value, "Sorry but this routine has not been ported to non built-in integers as it relies on a reinterpret_cast.");
 282     typedef typename boost::make_unsigned<IntType>::type unsigned_array[n];
 283     seed_array_int_impl<w>(seq, reinterpret_cast<unsigned_array&>(x));
 284 }
 285
 286 template<int w, std::size_t n, class SeedSeq, class IntType>
 287 inline void seed_array_int_impl(SeedSeq& seq, IntType (&x)[n], boost::false_type)
 288 {
 289     seed_array_int_impl<w>(seq, x);
 290 }
 291
 292 template<int w, std::size_t n, class SeedSeq, class IntType>
 293 inline void seed_array_int(SeedSeq& seq, IntType (&x)[n])
 294 {
 295     seed_array_int_impl<w>(seq, x, boost::random::traits::is_signed<IntType>());
 296 }
 297
 298 template<int w, std::size_t n, class Iter, class UIntType>
 299 void fill_array_int_impl(Iter& first, Iter last, UIntType (&x)[n])
 300 {
 301     for(std::size_t j = 0; j < n; j++) {
 302         UIntType val = 0;
 303         for(std::size_t k = 0; k < (w+31)/32; ++k) {
 304             if(first == last) {
 305                 boost::throw_exception(std::invalid_argument("Not enough elements in call to seed."));
 306             }
 307             val += static_cast<UIntType>(*first++) << 32*k;
 308         }
 309         x[j] = val & ::boost::low_bits_mask_t<w>::sig_bits;
 310     }
 311 }
 312
 313 template<int w, std::size_t n, class Iter, class IntType>
 314 inline void fill_array_int_impl(Iter& first, Iter last, IntType (&x)[n], boost::true_type)
 315 {
 316     BOOST_STATIC_ASSERT_MSG(boost::is_integral<IntType>::value, "Sorry but this routine has not been ported to non built-in integers as it relies on a reinterpret_cast.");
 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::false_type)
 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::random::traits::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) boost::throw_exception(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) boost::throw_exception(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