random/detail/seed_impl.hpp

0001 /* boost random/detail/seed.hpp header file
0002  *
0003  * Copyright Steven Watanabe 2009
0004  * Distributed under the Boost Software License, Version 1.0. (See
0005  * accompanying file LICENSE_1_0.txt or copy at
0006  * http://www.boost.org/LICENSE_1_0.txt)
0007  *
0008  * See http://www.boost.org for most recent version including documentation.
0009  *
0010  * $Id$
0011  */
0012
0013 #ifndef BOOST_RANDOM_DETAIL_SEED_IMPL_HPP
0014 #define BOOST_RANDOM_DETAIL_SEED_IMPL_HPP
0015
0016 #include <stdexcept>
0017 #include <boost/cstdint.hpp>
0018 #include <boost/throw_exception.hpp>
0019 #include <boost/config/no_tr1/cmath.hpp>
0020 #include <boost/integer/integer_mask.hpp>
0021 #include <boost/integer/static_log2.hpp>
0022 #include <boost/random/traits.hpp>
0023 #include <boost/random/detail/const_mod.hpp>
0024 #include <boost/random/detail/integer_log2.hpp>
0025 #include <boost/random/detail/signed_unsigned_tools.hpp>
0026 #include <boost/random/detail/generator_bits.hpp>
0027 #include <boost/type_traits/conditional.hpp>
0028 #include <boost/type_traits/integral_constant.hpp>
0029
0030 #include <boost/random/detail/disable_warnings.hpp>
0031
0032 namespace boost {
0033 namespace random {
0034 namespace detail {
0035
0036 // finds the seed type of an engine, given its
0037 // result_type.  If the result_type is integral
0038 // the seed type is the same.  If the result_type
0039 // is floating point, the seed type is uint32_t
0040 template<class T>
0041 struct seed_type
0042 {
0043     typedef typename boost::conditional<boost::is_integral<T>::value,
0044         T,
0045         boost::uint32_t
0046     >::type type;
0047 };
0048
0049 template<int N>
0050 struct const_pow_impl
0051 {
0052     template<class T>
0053     static T call(T arg, int n, T result)
0054     {
0055         return const_pow_impl<N / 2>::call(T(arg * arg), n / 2,
0056             n%2 == 0? result : T(result * arg));
0057     }
0058 };
0059
0060 template<>
0061 struct const_pow_impl<0>
0062 {
0063     template<class T>
0064     static T call(T, int, T result)
0065     {
0066         return result;
0067     }
0068 };
0069
0070 // requires N is an upper bound on n
0071 template<int N, class T>
0072 inline T const_pow(T arg, int n) { return const_pow_impl<N>::call(arg, n, T(1)); }
0073
0074 template<class T>
0075 inline T pow2(int n)
0076 {
0077     typedef unsigned int_type;
0078     const int max_bits = std::numeric_limits<int_type>::digits;
0079     T multiplier = T(int_type(1) << (max_bits - 1)) * 2;
0080     return (int_type(1) << (n % max_bits)) *
0081         const_pow<std::numeric_limits<T>::digits / max_bits>(multiplier, n / max_bits);
0082 }
0083
0084 template<class Engine, class Iter>
0085 void generate_from_real(Engine& eng, Iter begin, Iter end)
0086 {
0087     using std::fmod;
0088     typedef typename Engine::result_type RealType;
0089     const int Bits = detail::generator_bits<Engine>::value();
0090     int remaining_bits = 0;
0091     boost::uint_least32_t saved_bits = 0;
0092     RealType multiplier = pow2<RealType>( Bits);
0093     RealType mult32 = RealType(4294967296.0); // 2^32
0094     while(true) {
0095         RealType val = eng() * multiplier;
0096         int available_bits = Bits;
0097         // Make sure the compiler can optimize this out
0098         // if it isn't possible.
0099         if(Bits < 32 && available_bits < 32 - remaining_bits) {
0100             saved_bits |= boost::uint_least32_t(val) << remaining_bits;
0101             remaining_bits += Bits;
0102         } else {
0103             // If Bits < 32, then remaining_bits != 0, since
0104             // if remaining_bits == 0, available_bits < 32 - 0,
0105             // and we won't get here to begin with.
0106             if(Bits < 32 || remaining_bits != 0) {
0107                 boost::uint_least32_t divisor =
0108                     (boost::uint_least32_t(1) << (32 - remaining_bits));
0109                 boost::uint_least32_t extra_bits = boost::uint_least32_t(fmod(val, mult32)) & (divisor - 1);
0110                 val = val / divisor;
0111                 *begin++ = saved_bits | (extra_bits << remaining_bits);
0112                 if(begin == end) return;
0113                 available_bits -= 32 - remaining_bits;
0114                 remaining_bits = 0;
0115             }
0116             // If Bits < 32 we should never enter this loop
0117             if(Bits >= 32) {
0118                 for(; available_bits >= 32; available_bits -= 32) {
0119                     boost::uint_least32_t word = boost::uint_least32_t(fmod(val, mult32));
0120                     val /= mult32;
0121                     *begin++ = word;
0122                     if(begin == end) return;
0123                 }
0124             }
0125             remaining_bits = available_bits;
0126             saved_bits = static_cast<boost::uint_least32_t>(val);
0127         }
0128     }
0129 }
0130
0131 template<class Engine, class Iter>
0132 void generate_from_int(Engine& eng, Iter begin, Iter end)
0133 {
0134     typedef typename Engine::result_type IntType;
0135     typedef typename boost::random::traits::make_unsigned<IntType>::type unsigned_type;
0136     int remaining_bits = 0;
0137     boost::uint_least32_t saved_bits = 0;
0138     unsigned_type range = boost::random::detail::subtract<IntType>()((eng.max)(), (eng.min)());
0139
0140     int bits =
0141         (range == (std::numeric_limits<unsigned_type>::max)()) ?
0142             std::numeric_limits<unsigned_type>::digits :
0143             detail::integer_log2(range + 1);
0144
0145     {
0146         int discarded_bits = detail::integer_log2(bits);
0147         unsigned_type excess = (range + 1) >> (bits - discarded_bits);
0148         if(excess != 0) {
0149             int extra_bits = detail::integer_log2((excess - 1) ^ excess);
0150             bits = bits - discarded_bits + extra_bits;
0151         }
0152     }
0153
0154     unsigned_type mask = (static_cast<unsigned_type>(2) << (bits - 1)) - 1;
0155     unsigned_type limit = ((range + 1) & ~mask) - 1;
0156
0157     while(true) {
0158         unsigned_type val;
0159         do {
0160             val = boost::random::detail::subtract<IntType>()(eng(), (eng.min)());
0161         } while(limit != range && val > limit);
0162         val &= mask;
0163         int available_bits = bits;
0164         if(available_bits == 32) {
0165             *begin++ = static_cast<boost::uint_least32_t>(val) & 0xFFFFFFFFu;
0166             if(begin == end) return;
0167         } else if(available_bits % 32 == 0) {
0168             for(int i = 0; i < available_bits / 32; ++i) {
0169                 boost::uint_least32_t word = boost::uint_least32_t(val) & 0xFFFFFFFFu;
0170                 int suppress_warning = (bits >= 32);
0171                 BOOST_ASSERT(suppress_warning == 1);
0172                 val >>= (32 * suppress_warning);
0173                 *begin++ = word;
0174                 if(begin == end) return;
0175             }
0176         } else if(bits < 32 && available_bits < 32 - remaining_bits) {
0177             saved_bits |= boost::uint_least32_t(val) << remaining_bits;
0178             remaining_bits += bits;
0179         } else {
0180             if(bits < 32 || remaining_bits != 0) {
0181                 boost::uint_least32_t extra_bits = boost::uint_least32_t(val) & ((boost::uint_least32_t(1) << (32 - remaining_bits)) - 1);
0182                 val >>= 32 - remaining_bits;
0183                 *begin++ = saved_bits | (extra_bits << remaining_bits);
0184                 if(begin == end) return;
0185                 available_bits -= 32 - remaining_bits;
0186                 remaining_bits = 0;
0187             }
0188             if(bits >= 32) {
0189                 for(; available_bits >= 32; available_bits -= 32) {
0190                     boost::uint_least32_t word = boost::uint_least32_t(val) & 0xFFFFFFFFu;
0191                     int suppress_warning = (bits >= 32);
0192                     BOOST_ASSERT(suppress_warning == 1);
0193                     val >>= (32 * suppress_warning);
0194                     *begin++ = word;
0195                     if(begin == end) return;
0196                 }
0197             }
0198             remaining_bits = available_bits;
0199             saved_bits = static_cast<boost::uint_least32_t>(val);
0200         }
0201     }
0202 }
0203
0204 template<class Engine, class Iter>
0205 void generate_impl(Engine& eng, Iter first, Iter last, boost::true_type)
0206 {
0207     return detail::generate_from_int(eng, first, last);
0208 }
0209
0210 template<class Engine, class Iter>
0211 void generate_impl(Engine& eng, Iter first, Iter last, boost::false_type)
0212 {
0213     return detail::generate_from_real(eng, first, last);
0214 }
0215
0216 template<class Engine, class Iter>
0217 void generate(Engine& eng, Iter first, Iter last)
0218 {
0219     return detail::generate_impl(eng, first, last, boost::random::traits::is_integral<typename Engine::result_type>());
0220 }
0221
0222
0223
0224 template<class IntType, IntType m, class SeedSeq>
0225 IntType seed_one_int(SeedSeq& seq)
0226 {
0227     static const int log = ::boost::conditional<(m == 0),
0228         ::boost::integral_constant<int, (::std::numeric_limits<IntType>::digits)>,
0229         ::boost::static_log2<m> >::type::value;
0230     static const int k =
0231         (log + ((~(static_cast<IntType>(2) << (log - 1)) & m)? 32 : 31)) / 32;
0232     ::boost::uint_least32_t array[log / 32 + 4];
0233     seq.generate(&array[0], &array[0] + k + 3);
0234     IntType s = 0;
0235     for(int j = 0; j < k; ++j) {
0236         IntType digit = const_mod<IntType, m>::apply(IntType(array[j+3]));
0237         IntType mult = IntType(1) << 32*j;
0238         s = const_mod<IntType, m>::mult_add(mult, digit, s);
0239     }
0240     return s;
0241 }
0242
0243 template<class IntType, IntType m, class Iter>
0244 IntType get_one_int(Iter& first, Iter last)
0245 {
0246     static const int log = ::boost::conditional<(m == 0),
0247         ::boost::integral_constant<int, (::std::numeric_limits<IntType>::digits)>,
0248         ::boost::static_log2<m> >::type::value;
0249     static const int k =
0250         (log + ((~(static_cast<IntType>(2) << (log - 1)) & m)? 32 : 31)) / 32;
0251     IntType s = 0;
0252     for(int j = 0; j < k; ++j) {
0253         if(first == last) {
0254             boost::throw_exception(::std::invalid_argument("Not enough elements in call to seed."));
0255         }
0256         IntType digit = const_mod<IntType, m>::apply(IntType(*first++));
0257         IntType mult = IntType(1) << 32*j;
0258         s = const_mod<IntType, m>::mult_add(mult, digit, s);
0259     }
0260     return s;
0261 }
0262
0263 // TODO: work in-place whenever possible
0264 template<int w, std::size_t n, class SeedSeq, class UIntType>
0265 void seed_array_int_impl(SeedSeq& seq, UIntType (&x)[n])
0266 {
0267     boost::uint_least32_t storage[((w+31)/32) * n];
0268     seq.generate(&storage[0], &storage[0] + ((w+31)/32) * n);
0269     for(std::size_t j = 0; j < n; j++) {
0270         UIntType val = 0;
0271         for(std::size_t k = 0; k < (w+31)/32; ++k) {
0272             val += static_cast<UIntType>(storage[(w+31)/32*j + k]) << 32*k;
0273         }
0274         x[j] = val & ::boost::low_bits_mask_t<w>::sig_bits;
0275     }
0276 }
0277
0278 template<int w, std::size_t n, class SeedSeq, class IntType>
0279 inline void seed_array_int_impl(SeedSeq& seq, IntType (&x)[n], boost::true_type)
0280 {
0281     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.");
0282     typedef typename boost::make_unsigned<IntType>::type unsigned_array[n];
0283     seed_array_int_impl<w>(seq, reinterpret_cast<unsigned_array&>(x));
0284 }
0285
0286 template<int w, std::size_t n, class SeedSeq, class IntType>
0287 inline void seed_array_int_impl(SeedSeq& seq, IntType (&x)[n], boost::false_type)
0288 {
0289     seed_array_int_impl<w>(seq, x);
0290 }
0291
0292 template<int w, std::size_t n, class SeedSeq, class IntType>
0293 inline void seed_array_int(SeedSeq& seq, IntType (&x)[n])
0294 {
0295     seed_array_int_impl<w>(seq, x, boost::random::traits::is_signed<IntType>());
0296 }
0297
0298 template<int w, std::size_t n, class Iter, class UIntType>
0299 void fill_array_int_impl(Iter& first, Iter last, UIntType (&x)[n])
0300 {
0301     for(std::size_t j = 0; j < n; j++) {
0302         UIntType val = 0;
0303         for(std::size_t k = 0; k < (w+31)/32; ++k) {
0304             if(first == last) {
0305                 boost::throw_exception(std::invalid_argument("Not enough elements in call to seed."));
0306             }
0307             val += static_cast<UIntType>(*first++) << 32*k;
0308         }
0309         x[j] = val & ::boost::low_bits_mask_t<w>::sig_bits;
0310     }
0311 }
0312
0313 template<int w, std::size_t n, class Iter, class IntType>
0314 inline void fill_array_int_impl(Iter& first, Iter last, IntType (&x)[n], boost::true_type)
0315 {
0316     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.");
0317     typedef typename boost::make_unsigned<IntType>::type unsigned_array[n];
0318     fill_array_int_impl<w>(first, last, reinterpret_cast<unsigned_array&>(x));
0319 }
0320
0321 template<int w, std::size_t n, class Iter, class IntType>
0322 inline void fill_array_int_impl(Iter& first, Iter last, IntType (&x)[n], boost::false_type)
0323 {
0324     fill_array_int_impl<w>(first, last, x);
0325 }
0326
0327 template<int w, std::size_t n, class Iter, class IntType>
0328 inline void fill_array_int(Iter& first, Iter last, IntType (&x)[n])
0329 {
0330     fill_array_int_impl<w>(first, last, x, boost::random::traits::is_signed<IntType>());
0331 }
0332
0333 template<int w, std::size_t n, class RealType>
0334 void seed_array_real_impl(const boost::uint_least32_t* storage, RealType (&x)[n])
0335 {
0336     boost::uint_least32_t mask = ~((~boost::uint_least32_t(0)) << (w%32));
0337     RealType two32 = 4294967296.0;
0338     const RealType divisor = RealType(1)/detail::pow2<RealType>(w);
0339     unsigned int j;
0340     for(j = 0; j < n; ++j) {
0341         RealType val = RealType(0);
0342         RealType mult = divisor;
0343         for(int k = 0; k < w/32; ++k) {
0344             val += *storage++ * mult;
0345             mult *= two32;
0346         }
0347         if(mask != 0) {
0348             val += (*storage++ & mask) * mult;
0349         }
0350         BOOST_ASSERT(val >= 0);
0351         BOOST_ASSERT(val < 1);
0352         x[j] = val;
0353     }
0354 }
0355
0356 template<int w, std::size_t n, class SeedSeq, class RealType>
0357 void seed_array_real(SeedSeq& seq, RealType (&x)[n])
0358 {
0359     using std::pow;
0360     boost::uint_least32_t storage[((w+31)/32) * n];
0361     seq.generate(&storage[0], &storage[0] + ((w+31)/32) * n);
0362     seed_array_real_impl<w>(storage, x);
0363 }
0364
0365 template<int w, std::size_t n, class Iter, class RealType>
0366 void fill_array_real(Iter& first, Iter last, RealType (&x)[n])
0367 {
0368     boost::uint_least32_t mask = ~((~boost::uint_least32_t(0)) << (w%32));
0369     RealType two32 = 4294967296.0;
0370     const RealType divisor = RealType(1)/detail::pow2<RealType>(w);
0371     unsigned int j;
0372     for(j = 0; j < n; ++j) {
0373         RealType val = RealType(0);
0374         RealType mult = divisor;
0375         for(int k = 0; k < w/32; ++k, ++first) {
0376             if(first == last) boost::throw_exception(std::invalid_argument("Not enough elements in call to seed."));
0377             val += *first * mult;
0378             mult *= two32;
0379         }
0380         if(mask != 0) {
0381             if(first == last) boost::throw_exception(std::invalid_argument("Not enough elements in call to seed."));
0382             val += (*first & mask) * mult;
0383             ++first;
0384         }
0385         BOOST_ASSERT(val >= 0);
0386         BOOST_ASSERT(val < 1);
0387         x[j] = val;
0388     }
0389 }
0390
0391 }
0392 }
0393 }
0394
0395 #include <boost/random/detail/enable_warnings.hpp>
0396
0397 #endif