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