summaryrefslogtreecommitdiff
path: root/boost/random
diff options
context:
space:
mode:
Diffstat (limited to 'boost/random')
-rw-r--r--boost/random/detail/const_mod.hpp18
-rw-r--r--boost/random/detail/disable_warnings.hpp1
-rw-r--r--boost/random/detail/polynomial.hpp384
-rw-r--r--boost/random/detail/seed.hpp1
-rw-r--r--boost/random/detail/seed_impl.hpp12
-rw-r--r--boost/random/generate_canonical.hpp1
-rw-r--r--boost/random/mersenne_twister.hpp121
-rw-r--r--boost/random/non_central_chi_squared_distribution.hpp221
-rw-r--r--boost/random/normal_distribution.hpp1
-rw-r--r--boost/random/uniform_int_distribution.hpp1
-rw-r--r--boost/random/uniform_on_sphere.hpp2
-rw-r--r--boost/random/uniform_real_distribution.hpp1
-rw-r--r--boost/random/uniform_smallint.hpp1
13 files changed, 747 insertions, 18 deletions
diff --git a/boost/random/detail/const_mod.hpp b/boost/random/detail/const_mod.hpp
index 07f4ea7d84..e0a43ab249 100644
--- a/boost/random/detail/const_mod.hpp
+++ b/boost/random/detail/const_mod.hpp
@@ -36,9 +36,9 @@ public:
if(((unsigned_m() - 1) & unsigned_m()) == 0)
return (unsigned_type(x)) & (unsigned_m() - 1);
else {
- IntType supress_warnings = (m == 0);
- BOOST_ASSERT(supress_warnings == 0);
- return x % (m + supress_warnings);
+ IntType suppress_warnings = (m == 0);
+ BOOST_ASSERT(suppress_warnings == 0);
+ return x % (m + suppress_warnings);
}
}
@@ -77,9 +77,9 @@ public:
else if(a == 0)
return c;
else if(m <= (traits::const_max-c)/a) { // i.e. a*m+c <= max
- IntType supress_warnings = (m == 0);
- BOOST_ASSERT(supress_warnings == 0);
- return (a*x+c) % (m + supress_warnings);
+ IntType suppress_warnings = (m == 0);
+ BOOST_ASSERT(suppress_warnings == 0);
+ return (a*x+c) % (m + suppress_warnings);
} else
return add(mult(a, x), c);
}
@@ -108,9 +108,9 @@ private:
static IntType mult_small(IntType a, IntType x)
{
- IntType supress_warnings = (m == 0);
- BOOST_ASSERT(supress_warnings == 0);
- return a*x % (m + supress_warnings);
+ IntType suppress_warnings = (m == 0);
+ BOOST_ASSERT(suppress_warnings == 0);
+ return a*x % (m + suppress_warnings);
}
static IntType mult_schrage(IntType a, IntType value)
diff --git a/boost/random/detail/disable_warnings.hpp b/boost/random/detail/disable_warnings.hpp
index d875004c44..4582dcb1a7 100644
--- a/boost/random/detail/disable_warnings.hpp
+++ b/boost/random/detail/disable_warnings.hpp
@@ -20,6 +20,7 @@
#pragma warning(disable:4512)
#pragma warning(disable:4127)
#pragma warning(disable:4724)
+#pragma warning(disable:4800) // 'int' : forcing value to bool 'true' or 'false' (performance warning)
#endif
#if defined(BOOST_GCC) && BOOST_GCC >= 40600
diff --git a/boost/random/detail/polynomial.hpp b/boost/random/detail/polynomial.hpp
new file mode 100644
index 0000000000..a8c4b269f1
--- /dev/null
+++ b/boost/random/detail/polynomial.hpp
@@ -0,0 +1,384 @@
+/* boost random/detail/polynomial.hpp header file
+ *
+ * Copyright Steven Watanabe 2014
+ * Distributed under the Boost Software License, Version 1.0. (See
+ * accompanying file LICENSE_1_0.txt or copy at
+ * http://www.boost.org/LICENSE_1_0.txt)
+ *
+ * See http://www.boost.org for most recent version including documentation.
+ *
+ * $Id$
+ */
+
+#ifndef BOOST_RANDOM_DETAIL_POLYNOMIAL_HPP
+#define BOOST_RANDOM_DETAIL_POLYNOMIAL_HPP
+
+#include <cstddef>
+#include <limits>
+#include <vector>
+#include <algorithm>
+#include <boost/assert.hpp>
+#include <boost/cstdint.hpp>
+
+namespace boost {
+namespace random {
+namespace detail {
+
+class polynomial_ops {
+public:
+ typedef unsigned long digit_t;
+
+ static void add(std::size_t size, const digit_t * lhs,
+ const digit_t * rhs, digit_t * output)
+ {
+ for(std::size_t i = 0; i < size; ++i) {
+ output[i] = lhs[i] ^ rhs[i];
+ }
+ }
+
+ static void add_shifted_inplace(std::size_t size, const digit_t * lhs,
+ digit_t * output, std::size_t shift)
+ {
+ if(shift == 0) {
+ add(size, lhs, output, output);
+ return;
+ }
+ std::size_t bits = std::numeric_limits<digit_t>::digits;
+ digit_t prev = 0;
+ for(std::size_t i = 0; i < size; ++i) {
+ digit_t tmp = lhs[i];
+ output[i] ^= (tmp << shift) | (prev >> (bits-shift));
+ prev = tmp;
+ }
+ output[size] ^= (prev >> (bits-shift));
+ }
+
+ static void multiply_simple(std::size_t size, const digit_t * lhs,
+ const digit_t * rhs, digit_t * output)
+ {
+ std::size_t bits = std::numeric_limits<digit_t>::digits;
+ for(std::size_t i = 0; i < 2*size; ++i) {
+ output[i] = 0;
+ }
+ for(std::size_t i = 0; i < size; ++i) {
+ for(std::size_t j = 0; j < bits; ++j) {
+ if((lhs[i] & (digit_t(1) << j)) != 0) {
+ add_shifted_inplace(size, rhs, output + i, j);
+ }
+ }
+ }
+ }
+
+ // memory requirements: (size - cutoff) * 4 + next_smaller
+ static void multiply_karatsuba(std::size_t size,
+ const digit_t * lhs, const digit_t * rhs,
+ digit_t * output)
+ {
+ if(size < 64) {
+ multiply_simple(size, lhs, rhs, output);
+ return;
+ }
+ // split in half
+ std::size_t cutoff = size/2;
+ multiply_karatsuba(cutoff, lhs, rhs, output);
+ multiply_karatsuba(size - cutoff, lhs + cutoff, rhs + cutoff,
+ output + cutoff*2);
+ std::vector<digit_t> local1(size - cutoff);
+ std::vector<digit_t> local2(size - cutoff);
+ // combine the digits for the inner multiply
+ add(cutoff, lhs, lhs + cutoff, &local1[0]);
+ if(size & 1) local1[cutoff] = lhs[size - 1];
+ add(cutoff, rhs + cutoff, rhs, &local2[0]);
+ if(size & 1) local2[cutoff] = rhs[size - 1];
+ std::vector<digit_t> local3((size - cutoff) * 2);
+ multiply_karatsuba(size - cutoff, &local1[0], &local2[0], &local3[0]);
+ add(cutoff * 2, output, &local3[0], &local3[0]);
+ add((size - cutoff) * 2, output + cutoff*2, &local3[0], &local3[0]);
+ // Finally, add the inner result
+ add((size - cutoff) * 2, output + cutoff, &local3[0], output + cutoff);
+ }
+
+ static void multiply_add_karatsuba(std::size_t size,
+ const digit_t * lhs, const digit_t * rhs,
+ digit_t * output)
+ {
+ std::vector<digit_t> buf(size * 2);
+ multiply_karatsuba(size, lhs, rhs, &buf[0]);
+ add(size * 2, &buf[0], output, output);
+ }
+
+ static void multiply(const digit_t * lhs, std::size_t lhs_size,
+ const digit_t * rhs, std::size_t rhs_size,
+ digit_t * output)
+ {
+ std::fill_n(output, lhs_size + rhs_size, digit_t(0));
+ multiply_add(lhs, lhs_size, rhs, rhs_size, output);
+ }
+
+ static void multiply_add(const digit_t * lhs, std::size_t lhs_size,
+ const digit_t * rhs, std::size_t rhs_size,
+ digit_t * output)
+ {
+ // split into pieces that can be passed to
+ // karatsuba multiply.
+ while(lhs_size != 0) {
+ if(lhs_size < rhs_size) {
+ std::swap(lhs, rhs);
+ std::swap(lhs_size, rhs_size);
+ }
+
+ multiply_add_karatsuba(rhs_size, lhs, rhs, output);
+
+ lhs += rhs_size;
+ lhs_size -= rhs_size;
+ output += rhs_size;
+ }
+ }
+
+ static void copy_bits(const digit_t * x, std::size_t low, std::size_t high,
+ digit_t * out)
+ {
+ const std::size_t bits = std::numeric_limits<digit_t>::digits;
+ std::size_t offset = low/bits;
+ x += offset;
+ low -= offset*bits;
+ high -= offset*bits;
+ std::size_t n = (high-low)/bits;
+ if(low == 0) {
+ for(std::size_t i = 0; i < n; ++i) {
+ out[i] = x[i];
+ }
+ } else {
+ for(std::size_t i = 0; i < n; ++i) {
+ out[i] = (x[i] >> low) | (x[i+1] << (bits-low));
+ }
+ }
+ if((high-low)%bits) {
+ digit_t low_mask = (digit_t(1) << ((high-low)%bits)) - 1;
+ digit_t result = (x[n] >> low);
+ if(low != 0 && (n+1)*bits < high) {
+ result |= (x[n+1] << (bits-low));
+ }
+ out[n] = (result & low_mask);
+ }
+ }
+
+ static void shift_left(digit_t * val, std::size_t size, std::size_t shift)
+ {
+ const std::size_t bits = std::numeric_limits<digit_t>::digits;
+ BOOST_ASSERT(shift > 0);
+ BOOST_ASSERT(shift < bits);
+ digit_t prev = 0;
+ for(std::size_t i = 0; i < size; ++i) {
+ digit_t tmp = val[i];
+ val[i] = (prev >> (bits - shift)) | (val[i] << shift);
+ prev = tmp;
+ }
+ }
+
+ static digit_t sqr(digit_t val) {
+ const std::size_t bits = std::numeric_limits<digit_t>::digits;
+ digit_t mask = (digit_t(1) << bits/2) - 1;
+ for(std::size_t i = bits; i > 1; i /= 2) {
+ val = ((val & ~mask) << i/2) | (val & mask);
+ mask = mask & (mask >> i/4);
+ mask = mask | (mask << i/2);
+ }
+ return val;
+ }
+
+ static void sqr(digit_t * val, std::size_t size)
+ {
+ const std::size_t bits = std::numeric_limits<digit_t>::digits;
+ digit_t mask = (digit_t(1) << bits/2) - 1;
+ for(std::size_t i = 0; i < size; ++i) {
+ digit_t x = val[size - i - 1];
+ val[(size - i - 1) * 2] = sqr(x & mask);
+ val[(size - i - 1) * 2 + 1] = sqr(x >> bits/2);
+ }
+ }
+
+ // optimized for the case when the modulus has few bits set.
+ struct sparse_mod {
+ sparse_mod(const digit_t * divisor, std::size_t divisor_bits)
+ {
+ const std::size_t bits = std::numeric_limits<digit_t>::digits;
+ _remainder_bits = divisor_bits - 1;
+ for(std::size_t i = 0; i < divisor_bits; ++i) {
+ if(divisor[i/bits] & (digit_t(1) << i%bits)) {
+ _bit_indices.push_back(i);
+ }
+ }
+ BOOST_ASSERT(_bit_indices.back() == divisor_bits - 1);
+ _bit_indices.pop_back();
+ if(_bit_indices.empty()) {
+ _block_bits = divisor_bits;
+ _lower_bits = 0;
+ } else {
+ _block_bits = divisor_bits - _bit_indices.back() - 1;
+ _lower_bits = _bit_indices.back() + 1;
+ }
+
+ _partial_quotient.resize((_block_bits + bits - 1)/bits);
+ }
+ void operator()(digit_t * dividend, std::size_t dividend_bits)
+ {
+ const std::size_t bits = std::numeric_limits<digit_t>::digits;
+ while(dividend_bits > _remainder_bits) {
+ std::size_t block_start = (std::max)(dividend_bits - _block_bits, _remainder_bits);
+ std::size_t block_size = (dividend_bits - block_start + bits - 1) / bits;
+ copy_bits(dividend, block_start, dividend_bits, &_partial_quotient[0]);
+ for(std::size_t i = 0; i < _bit_indices.size(); ++i) {
+ std::size_t pos = _bit_indices[i] + block_start - _remainder_bits;
+ add_shifted_inplace(block_size, &_partial_quotient[0], dividend + pos/bits, pos%bits);
+ }
+ add_shifted_inplace(block_size, &_partial_quotient[0], dividend + block_start/bits, block_start%bits);
+ dividend_bits = block_start;
+ }
+ }
+ std::vector<digit_t> _partial_quotient;
+ std::size_t _remainder_bits;
+ std::size_t _block_bits;
+ std::size_t _lower_bits;
+ std::vector<std::size_t> _bit_indices;
+ };
+
+ // base should have the same number of bits as mod
+ // base, and mod should both be able to hold a power
+ // of 2 >= mod_bits. out needs to be twice as large.
+ static void mod_pow_x(boost::uintmax_t exponent, const digit_t * mod, std::size_t mod_bits, digit_t * out)
+ {
+ const std::size_t bits = std::numeric_limits<digit_t>::digits;
+ const std::size_t n = (mod_bits + bits - 1) / bits;
+ const std::size_t highbit = mod_bits - 1;
+ if(exponent == 0) {
+ out[0] = 1;
+ std::fill_n(out + 1, n - 1, digit_t(0));
+ return;
+ }
+ boost::uintmax_t i = std::numeric_limits<boost::uintmax_t>::digits - 1;
+ while(((boost::uintmax_t(1) << i) & exponent) == 0) {
+ --i;
+ }
+ out[0] = 2;
+ std::fill_n(out + 1, n - 1, digit_t(0));
+ sparse_mod m(mod, mod_bits);
+ while(i--) {
+ sqr(out, n);
+ m(out, 2 * mod_bits - 1);
+ if((boost::uintmax_t(1) << i) & exponent) {
+ shift_left(out, n, 1);
+ if(out[highbit / bits] & (digit_t(1) << highbit%bits))
+ add(n, out, mod, out);
+ }
+ }
+ }
+};
+
+class polynomial
+{
+ typedef polynomial_ops::digit_t digit_t;
+public:
+ polynomial() : _size(0) {}
+ class reference {
+ public:
+ reference(digit_t &value, int idx)
+ : _value(value), _idx(idx) {}
+ operator bool() const { return (_value & (digit_t(1) << _idx)) != 0; }
+ reference& operator=(bool b)
+ {
+ if(b) {
+ _value |= (digit_t(1) << _idx);
+ } else {
+ _value &= ~(digit_t(1) << _idx);
+ }
+ return *this;
+ }
+ reference &operator^=(bool b)
+ {
+ _value ^= (digit_t(b) << _idx);
+ return *this;
+ }
+
+ reference &operator=(const reference &other)
+ {
+ return *this = static_cast<bool>(other);
+ }
+ private:
+ digit_t &_value;
+ int _idx;
+ };
+ reference operator[](std::size_t i)
+ {
+ static const std::size_t bits = std::numeric_limits<digit_t>::digits;
+ ensure_bit(i);
+ return reference(_storage[i/bits], i%bits);
+ }
+ bool operator[](std::size_t i) const
+ {
+ static const std::size_t bits = std::numeric_limits<digit_t>::digits;
+ if(i < size())
+ return (_storage[i/bits] & (digit_t(1) << (i%bits))) != 0;
+ else
+ return false;
+ }
+ std::size_t size() const
+ {
+ return _size;
+ }
+ void resize(std::size_t n)
+ {
+ static const std::size_t bits = std::numeric_limits<digit_t>::digits;
+ _storage.resize((n + bits - 1)/bits);
+ // clear the high order bits in case we're shrinking.
+ if(n%bits) {
+ _storage.back() &= ((digit_t(1) << (n%bits)) - 1);
+ }
+ _size = n;
+ }
+ friend polynomial operator*(const polynomial &lhs, const polynomial &rhs);
+ friend polynomial mod_pow_x(boost::uintmax_t exponent, polynomial mod);
+private:
+ std::vector<polynomial_ops::digit_t> _storage;
+ std::size_t _size;
+ void ensure_bit(std::size_t i)
+ {
+ if(i >= size()) {
+ resize(i + 1);
+ }
+ }
+ void normalize()
+ {
+ while(size() && (*this)[size() - 1] == 0)
+ resize(size() - 1);
+ }
+};
+
+inline polynomial operator*(const polynomial &lhs, const polynomial &rhs)
+{
+ polynomial result;
+ result._storage.resize(lhs._storage.size() + rhs._storage.size());
+ polynomial_ops::multiply(&lhs._storage[0], lhs._storage.size(),
+ &rhs._storage[0], rhs._storage.size(),
+ &result._storage[0]);
+ result._size = lhs._size + rhs._size;
+ return result;
+}
+
+inline polynomial mod_pow_x(boost::uintmax_t exponent, polynomial mod)
+{
+ polynomial result;
+ mod.normalize();
+ std::size_t mod_size = mod.size();
+ result._storage.resize(mod._storage.size() * 2);
+ result._size = mod.size() * 2;
+ polynomial_ops::mod_pow_x(exponent, &mod._storage[0], mod_size, &result._storage[0]);
+ result.resize(mod.size() - 1);
+ return result;
+}
+
+}
+}
+}
+
+#endif // BOOST_RANDOM_DETAIL_POLYNOMIAL_HPP
diff --git a/boost/random/detail/seed.hpp b/boost/random/detail/seed.hpp
index 55b2fa667d..557482ad4b 100644
--- a/boost/random/detail/seed.hpp
+++ b/boost/random/detail/seed.hpp
@@ -21,6 +21,7 @@
#include <boost/utility/enable_if.hpp>
#include <boost/type_traits/is_arithmetic.hpp>
+#include <boost/mpl/bool.hpp>
namespace boost {
namespace random {
diff --git a/boost/random/detail/seed_impl.hpp b/boost/random/detail/seed_impl.hpp
index f88cab299c..2451dbe521 100644
--- a/boost/random/detail/seed_impl.hpp
+++ b/boost/random/detail/seed_impl.hpp
@@ -170,9 +170,9 @@ void generate_from_int(Engine& eng, Iter begin, Iter end)
} else if(available_bits % 32 == 0) {
for(int i = 0; i < available_bits / 32; ++i) {
boost::uint_least32_t word = boost::uint_least32_t(val) & 0xFFFFFFFFu;
- int supress_warning = (bits >= 32);
- BOOST_ASSERT(supress_warning == 1);
- val >>= (32 * supress_warning);
+ int suppress_warning = (bits >= 32);
+ BOOST_ASSERT(suppress_warning == 1);
+ val >>= (32 * suppress_warning);
*begin++ = word;
if(begin == end) return;
}
@@ -191,9 +191,9 @@ void generate_from_int(Engine& eng, Iter begin, Iter end)
if(bits >= 32) {
for(; available_bits >= 32; available_bits -= 32) {
boost::uint_least32_t word = boost::uint_least32_t(val) & 0xFFFFFFFFu;
- int supress_warning = (bits >= 32);
- BOOST_ASSERT(supress_warning == 1);
- val >>= (32 * supress_warning);
+ int suppress_warning = (bits >= 32);
+ BOOST_ASSERT(suppress_warning == 1);
+ val >>= (32 * suppress_warning);
*begin++ = word;
if(begin == end) return;
}
diff --git a/boost/random/generate_canonical.hpp b/boost/random/generate_canonical.hpp
index 23ff43dfac..d8ff144d6c 100644
--- a/boost/random/generate_canonical.hpp
+++ b/boost/random/generate_canonical.hpp
@@ -19,6 +19,7 @@
#include <boost/config/no_tr1/cmath.hpp>
#include <boost/limits.hpp>
#include <boost/type_traits/is_integral.hpp>
+#include <boost/mpl/bool.hpp>
#include <boost/random/detail/signed_unsigned_tools.hpp>
#include <boost/random/detail/generator_bits.hpp>
diff --git a/boost/random/mersenne_twister.hpp b/boost/random/mersenne_twister.hpp
index 3878fee156..ce73e6825f 100644
--- a/boost/random/mersenne_twister.hpp
+++ b/boost/random/mersenne_twister.hpp
@@ -29,6 +29,9 @@
#include <boost/random/detail/seed.hpp>
#include <boost/random/detail/seed_impl.hpp>
#include <boost/random/detail/generator_seed_seq.hpp>
+#include <boost/random/detail/polynomial.hpp>
+
+#include <boost/random/detail/disable_warnings.hpp>
namespace boost {
namespace random {
@@ -199,8 +202,15 @@ public:
*/
void discard(boost::uintmax_t z)
{
- for(boost::uintmax_t j = 0; j < z; ++j) {
- (*this)();
+#ifndef BOOST_RANDOM_MERSENNE_TWISTER_DISCARD_THRESHOLD
+#define BOOST_RANDOM_MERSENNE_TWISTER_DISCARD_THRESHOLD 10000000
+#endif
+ if(z > BOOST_RANDOM_MERSENNE_TWISTER_DISCARD_THRESHOLD) {
+ discard_many(z);
+ } else {
+ for(boost::uintmax_t j = 0; j < z; ++j) {
+ (*this)();
+ }
}
}
@@ -368,6 +378,111 @@ private:
}
}
+ /**
+ * Optimized algorithm for large jumps.
+ *
+ * Hiroshi Haramoto, Makoto Matsumoto, and Pierre L'Ecuyer. 2008.
+ * A Fast Jump Ahead Algorithm for Linear Recurrences in a Polynomial
+ * Space. In Proceedings of the 5th international conference on
+ * Sequences and Their Applications (SETA '08).
+ * DOI=10.1007/978-3-540-85912-3_26
+ */
+ void discard_many(boost::uintmax_t z)
+ {
+ // Compute the minimal polynomial, phi(t)
+ // This depends only on the transition function,
+ // which is constant. The characteristic
+ // polynomial is the same as the minimal
+ // polynomial for a maximum period generator
+ // (which should be all specializations of
+ // mersenne_twister.) Even if it weren't,
+ // the characteristic polynomial is guaranteed
+ // to be a multiple of the minimal polynomial,
+ // which is good enough.
+ detail::polynomial phi = get_characteristic_polynomial();
+
+ // calculate g(t) = t^z % phi(t)
+ detail::polynomial g = mod_pow_x(z, phi);
+
+ // h(s_0, t) = \sum_{i=0}^{2k-1}o(s_i)t^{2k-i-1}
+ detail::polynomial h;
+ const std::size_t num_bits = w*n - r;
+ for(std::size_t j = 0; j < num_bits * 2; ++j) {
+ // Yes, we're advancing the generator state
+ // here, but it doesn't matter because
+ // we're going to overwrite it completely
+ // in reconstruct_state.
+ if(i >= n) twist();
+ h[2*num_bits - j - 1] = x[i++] & UIntType(1);
+ }
+ // g(t)h(s_0, t)
+ detail::polynomial gh = g * h;
+ detail::polynomial result;
+ for(std::size_t j = 0; j <= num_bits; ++j) {
+ result[j] = gh[2*num_bits - j - 1];
+ }
+ reconstruct_state(result);
+ }
+ static detail::polynomial get_characteristic_polynomial()
+ {
+ const std::size_t num_bits = w*n - r;
+ detail::polynomial helper;
+ helper[num_bits - 1] = 1;
+ mersenne_twister_engine tmp;
+ tmp.reconstruct_state(helper);
+ // Skip the first num_bits elements, since we
+ // already know what they are.
+ for(std::size_t j = 0; j < num_bits; ++j) {
+ if(tmp.i >= n) tmp.twist();
+ if(j == num_bits - 1)
+ assert((tmp.x[tmp.i] & 1) == 1);
+ else
+ assert((tmp.x[tmp.i] & 1) == 0);
+ ++tmp.i;
+ }
+ detail::polynomial phi;
+ phi[num_bits] = 1;
+ detail::polynomial next_bits = tmp.as_polynomial(num_bits);
+ for(std::size_t j = 0; j < num_bits; ++j) {
+ int val = next_bits[j] ^ phi[num_bits-j-1];
+ phi[num_bits-j-1] = val;
+ if(val) {
+ for(std::size_t k = j + 1; k < num_bits; ++k) {
+ phi[num_bits-k-1] ^= next_bits[k-j-1];
+ }
+ }
+ }
+ return phi;
+ }
+ detail::polynomial as_polynomial(std::size_t size) {
+ detail::polynomial result;
+ for(std::size_t j = 0; j < size; ++j) {
+ if(i >= n) twist();
+ result[j] = x[i++] & UIntType(1);
+ }
+ return result;
+ }
+ void reconstruct_state(const detail::polynomial& p)
+ {
+ const UIntType upper_mask = (~static_cast<UIntType>(0)) << r;
+ const UIntType lower_mask = ~upper_mask;
+ const std::size_t num_bits = w*n - r;
+ for(std::size_t j = num_bits - n + 1; j <= num_bits; ++j)
+ x[j % n] = p[j];
+
+ UIntType y0 = 0;
+ for(std::size_t j = num_bits + 1; j >= n - 1; --j) {
+ UIntType y1 = x[j % n] ^ x[(j + m) % n];
+ if(p[j - n + 1])
+ y1 = (y1 ^ a) << UIntType(1) | UIntType(1);
+ else
+ y1 = y1 << UIntType(1);
+ x[(j + 1) % n] = (y0 & upper_mask) | (y1 & lower_mask);
+ y0 = y1;
+ }
+ i = 0;
+ }
+
/// \endcond
// state representation: next output is o(x(i))
@@ -562,4 +677,6 @@ BOOST_RANDOM_PTR_HELPER_SPEC(boost::mt11213b)
BOOST_RANDOM_PTR_HELPER_SPEC(boost::mt19937)
BOOST_RANDOM_PTR_HELPER_SPEC(boost::mt19937_64)
+#include <boost/random/detail/enable_warnings.hpp>
+
#endif // BOOST_RANDOM_MERSENNE_TWISTER_HPP
diff --git a/boost/random/non_central_chi_squared_distribution.hpp b/boost/random/non_central_chi_squared_distribution.hpp
new file mode 100644
index 0000000000..28c9ff6d9a
--- /dev/null
+++ b/boost/random/non_central_chi_squared_distribution.hpp
@@ -0,0 +1,221 @@
+/* boost random/non_central_chi_squared_distribution.hpp header file
+ *
+ * Copyright Thijs van den Berg 2014
+ *
+ * Distributed under the Boost Software License, Version 1.0. (See
+ * accompanying file LICENSE_1_0.txt or copy at
+ * http://www.boost.org/LICENSE_1_0.txt)
+ *
+ * See http://www.boost.org for most recent version including documentation.
+ *
+ * $Id$
+ */
+
+#ifndef BOOST_RANDOM_NON_CENTRAL_CHI_SQUARED_DISTRIBUTION_HPP
+#define BOOST_RANDOM_NON_CENTRAL_CHI_SQUARED_DISTRIBUTION_HPP
+
+#include <boost/config/no_tr1/cmath.hpp>
+#include <iosfwd>
+#include <istream>
+#include <boost/limits.hpp>
+#include <boost/random/detail/config.hpp>
+#include <boost/random/detail/operators.hpp>
+#include <boost/random/uniform_real_distribution.hpp>
+#include <boost/random/normal_distribution.hpp>
+#include <boost/random/chi_squared_distribution.hpp>
+#include <boost/random/poisson_distribution.hpp>
+
+namespace boost {
+namespace random {
+
+/**
+ * The noncentral chi-squared distribution is a real valued distribution with
+ * two parameter, @c k and @c lambda. The distribution produces values > 0.
+ *
+ * This is the distribution of the sum of squares of k Normal distributed
+ * variates each with variance one and \f$\lambda\f$ the sum of squares of the
+ * normal means.
+ *
+ * The distribution function is
+ * \f$\displaystyle P(x) = \frac{1}{2} e^{-(x+\lambda)/2} \left( \frac{x}{\lambda} \right)^{k/4-1/2} I_{k/2-1}( \sqrt{\lambda x} )\f$.
+ * where \f$\displaystyle I_\nu(z)\f$ is a modified Bessel function of the
+ * first kind.
+ *
+ * The algorithm is taken from
+ *
+ * @blockquote
+ * "Monte Carlo Methods in Financial Engineering", Paul Glasserman,
+ * 2003, XIII, 596 p, Stochastic Modelling and Applied Probability, Vol. 53,
+ * ISBN 978-0-387-21617-1, p 124, Fig. 3.5.
+ * @endblockquote
+ */
+template <typename RealType = double>
+class non_central_chi_squared_distribution {
+public:
+ typedef RealType result_type;
+ typedef RealType input_type;
+
+ class param_type {
+ public:
+ typedef non_central_chi_squared_distribution distribution_type;
+
+ /**
+ * Constructs the parameters of a non_central_chi_squared_distribution.
+ * @c k and @c lambda are the parameter of the distribution.
+ *
+ * Requires: k > 0 && lambda > 0
+ */
+ explicit
+ param_type(RealType k_arg = RealType(1), RealType lambda_arg = RealType(1))
+ : _k(k_arg), _lambda(lambda_arg)
+ {
+ BOOST_ASSERT(k_arg > RealType(0));
+ BOOST_ASSERT(lambda_arg > RealType(0));
+ }
+
+ /** Returns the @c k parameter of the distribution */
+ RealType k() const { return _k; }
+
+ /** Returns the @c lambda parameter of the distribution */
+ RealType lambda() const { return _lambda; }
+
+ /** Writes the parameters of the distribution to a @c std::ostream. */
+ BOOST_RANDOM_DETAIL_OSTREAM_OPERATOR(os, param_type, parm)
+ {
+ os << parm._k << ' ' << parm._lambda;
+ return os;
+ }
+
+ /** Reads the parameters of the distribution from a @c std::istream. */
+ BOOST_RANDOM_DETAIL_ISTREAM_OPERATOR(is, param_type, parm)
+ {
+ is >> parm._k >> std::ws >> parm._lambda;
+ return is;
+ }
+
+ /** Returns true if the parameters have the same values. */
+ BOOST_RANDOM_DETAIL_EQUALITY_OPERATOR(param_type, lhs, rhs)
+ { return lhs._k == rhs._k && lhs._lambda == rhs._lambda; }
+
+ /** Returns true if the parameters have different values. */
+ BOOST_RANDOM_DETAIL_INEQUALITY_OPERATOR(param_type)
+
+ private:
+ RealType _k;
+ RealType _lambda;
+ };
+
+ /**
+ * Construct a @c non_central_chi_squared_distribution object. @c k and
+ * @c lambda are the parameter of the distribution.
+ *
+ * Requires: k > 0 && lambda > 0
+ */
+ explicit
+ non_central_chi_squared_distribution(RealType k_arg = RealType(1), RealType lambda_arg = RealType(1))
+ : _param(k_arg, lambda_arg)
+ {
+ BOOST_ASSERT(k_arg > RealType(0));
+ BOOST_ASSERT(lambda_arg > RealType(0));
+ }
+
+ /**
+ * Construct a @c non_central_chi_squared_distribution object from the parameter.
+ */
+ explicit
+ non_central_chi_squared_distribution(const param_type& parm)
+ : _param( parm )
+ { }
+
+ /**
+ * Returns a random variate distributed according to the
+ * non central chi squared distribution specified by @c param.
+ */
+ template<typename URNG>
+ RealType operator()(URNG& eng, const param_type& parm) const
+ { return non_central_chi_squared_distribution(parm)(eng); }
+
+ /**
+ * Returns a random variate distributed according to the
+ * non central chi squared distribution.
+ */
+ template<typename URNG>
+ RealType operator()(URNG& eng)
+ {
+ using std::sqrt;
+ if (_param.k() > 1) {
+ boost::random::normal_distribution<RealType> n_dist;
+ boost::random::chi_squared_distribution<RealType> c_dist(_param.k() - RealType(1));
+ RealType _z = n_dist(eng);
+ RealType _x = c_dist(eng);
+ RealType term1 = _z + sqrt(_param.lambda());
+ return term1*term1 + _x;
+ }
+ else {
+ boost::random::poisson_distribution<> p_dist(_param.lambda()/RealType(2));
+ boost::random::poisson_distribution<>::result_type _p = p_dist(eng);
+ boost::random::chi_squared_distribution<RealType> c_dist(_param.k() + RealType(2)*_p);
+ return c_dist(eng);
+ }
+ }
+
+ /** Returns the @c k parameter of the distribution. */
+ RealType k() const { return _param.k(); }
+
+ /** Returns the @c lambda parameter of the distribution. */
+ RealType lambda() const { return _param.lambda(); }
+
+ /** Returns the parameters of the distribution. */
+ param_type param() const { return _param; }
+
+ /** Sets parameters of the distribution. */
+ void param(const param_type& parm) { _param = parm; }
+
+ /** Resets the distribution, so that subsequent uses does not depend on values already produced by it.*/
+ void reset() {}
+
+ /** Returns the smallest value that the distribution can produce. */
+ RealType min BOOST_PREVENT_MACRO_SUBSTITUTION() const
+ { return RealType(0); }
+
+ /** Returns the largest value that the distribution can produce. */
+ RealType max BOOST_PREVENT_MACRO_SUBSTITUTION() const
+ { return (std::numeric_limits<RealType>::infinity)(); }
+
+ /** Writes the parameters of the distribution to a @c std::ostream. */
+ BOOST_RANDOM_DETAIL_OSTREAM_OPERATOR(os, non_central_chi_squared_distribution, dist)
+ {
+ os << dist.param();
+ return os;
+ }
+
+ /** reads the parameters of the distribution from a @c std::istream. */
+ BOOST_RANDOM_DETAIL_ISTREAM_OPERATOR(is, non_central_chi_squared_distribution, dist)
+ {
+ param_type parm;
+ if(is >> parm) {
+ dist.param(parm);
+ }
+ return is;
+ }
+
+ /** Returns true if two distributions have the same parameters and produce
+ the same sequence of random numbers given equal generators.*/
+ BOOST_RANDOM_DETAIL_EQUALITY_OPERATOR(non_central_chi_squared_distribution, lhs, rhs)
+ { return lhs.param() == rhs.param(); }
+
+ /** Returns true if two distributions have different parameters and/or can produce
+ different sequences of random numbers given equal generators.*/
+ BOOST_RANDOM_DETAIL_INEQUALITY_OPERATOR(non_central_chi_squared_distribution)
+
+private:
+
+ /// @cond show_private
+ param_type _param;
+ /// @endcond
+};
+
+} // namespace random
+} // namespace boost
+
+#endif
diff --git a/boost/random/normal_distribution.hpp b/boost/random/normal_distribution.hpp
index 0be1298fcc..6fffc4aef5 100644
--- a/boost/random/normal_distribution.hpp
+++ b/boost/random/normal_distribution.hpp
@@ -33,6 +33,7 @@
#include <boost/random/uniform_01.hpp>
#include <boost/random/uniform_int_distribution.hpp>
#include <boost/random/exponential_distribution.hpp>
+#include <boost/mpl/bool.hpp>
namespace boost {
namespace random {
diff --git a/boost/random/uniform_int_distribution.hpp b/boost/random/uniform_int_distribution.hpp
index a175025edc..e7ef4b8248 100644
--- a/boost/random/uniform_int_distribution.hpp
+++ b/boost/random/uniform_int_distribution.hpp
@@ -30,6 +30,7 @@
#include <boost/random/detail/signed_unsigned_tools.hpp>
#include <boost/type_traits/make_unsigned.hpp>
#include <boost/type_traits/is_integral.hpp>
+#include <boost/mpl/bool.hpp>
namespace boost {
namespace random {
diff --git a/boost/random/uniform_on_sphere.hpp b/boost/random/uniform_on_sphere.hpp
index fdd4e950bb..0f744551a2 100644
--- a/boost/random/uniform_on_sphere.hpp
+++ b/boost/random/uniform_on_sphere.hpp
@@ -194,7 +194,7 @@ public:
{
uniform_01<RealType> uniform;
RealType sqsum;
- RealType x, y, z;
+ RealType x, y;
do {
x = uniform(eng) * 2 - 1;
y = uniform(eng) * 2 - 1;
diff --git a/boost/random/uniform_real_distribution.hpp b/boost/random/uniform_real_distribution.hpp
index 437ee69330..c53d85835b 100644
--- a/boost/random/uniform_real_distribution.hpp
+++ b/boost/random/uniform_real_distribution.hpp
@@ -24,6 +24,7 @@
#include <boost/random/detail/operators.hpp>
#include <boost/random/detail/signed_unsigned_tools.hpp>
#include <boost/type_traits/is_integral.hpp>
+#include <boost/mpl/bool.hpp>
namespace boost {
namespace random {
diff --git a/boost/random/uniform_smallint.hpp b/boost/random/uniform_smallint.hpp
index 4f5f437517..96f2fdc2c1 100644
--- a/boost/random/uniform_smallint.hpp
+++ b/boost/random/uniform_smallint.hpp
@@ -28,6 +28,7 @@
#include <boost/random/detail/signed_unsigned_tools.hpp>
#include <boost/random/uniform_01.hpp>
#include <boost/detail/workaround.hpp>
+#include <boost/mpl/bool.hpp>
namespace boost {
namespace random {
generated by cgit v1.2.3 (git 2.25.1) at 2024-06-17 09:49:58 +0000