diff --git a/Boost/boost/random.hpp b/Boost/boost/random.hpp new file mode 100644 index 000000000..4939e7d9b --- /dev/null +++ b/Boost/boost/random.hpp @@ -0,0 +1,84 @@ +/* boost random.hpp header file + * + * Copyright Jens Maurer 2000-2001 + * 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/libs/random for documentation. + * + * $Id: random.hpp 60755 2010-03-22 00:45:06Z steven_watanabe $ + * + * Revision history + * 2000-02-18 portability fixes (thanks to Beman Dawes) + * 2000-02-21 shuffle_output, inversive_congruential_schrage, + * generator_iterator, uniform_smallint + * 2000-02-23 generic modulus arithmetic helper, removed *_schrage classes, + * implemented Streamable and EqualityComparable concepts for + * generators, added Bernoulli distribution and Box-Muller + * transform + * 2000-03-01 cauchy, lognormal, triangle distributions; fixed + * uniform_smallint; renamed gaussian to normal distribution + * 2000-03-05 implemented iterator syntax for distribution functions + * 2000-04-21 removed some optimizations for better BCC/MSVC compatibility + * 2000-05-10 adapted to BCC and MSVC + * 2000-06-13 incorporated review results + * 2000-07-06 moved basic templates from namespace detail to random + * 2000-09-23 warning removals and int64 fixes (Ed Brey) + * 2000-09-24 added lagged_fibonacci generator (Matthias Troyer) + * 2001-02-18 moved to individual header files + */ + +#ifndef BOOST_RANDOM_HPP +#define BOOST_RANDOM_HPP + +// generators +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +namespace boost { + /** + * The specialization taus88 was suggested in + * + * @blockquote + * "Maximally Equidistributed Combined Tausworthe Generators", + * Pierre L'Ecuyer, Mathematics of Computation, Volume 65, + * Number 213, January 1996, Pages 203-213 + * @endblockquote + */ + typedef random::xor_combine, 0, + random::linear_feedback_shift, 0, 0>, 0, + random::linear_feedback_shift, 0, 0> taus88; +} // namespace boost + +// misc +#include + +// distributions +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#endif // BOOST_RANDOM_HPP diff --git a/Boost/boost/random/additive_combine.hpp b/Boost/boost/random/additive_combine.hpp new file mode 100644 index 000000000..3e8b0c312 --- /dev/null +++ b/Boost/boost/random/additive_combine.hpp @@ -0,0 +1,234 @@ +/* boost random/additive_combine.hpp header file + * + * Copyright Jens Maurer 2000-2001 + * 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: additive_combine.hpp 60755 2010-03-22 00:45:06Z steven_watanabe $ + * + * Revision history + * 2001-02-18 moved to individual header files + */ + +#ifndef BOOST_RANDOM_ADDITIVE_COMBINE_HPP +#define BOOST_RANDOM_ADDITIVE_COMBINE_HPP + +#include +#include // for std::min and std::max +#include +#include +#include +#include + +namespace boost { +namespace random { + +/** + * An instantiation of class template \additive_combine model a + * \pseudo_random_number_generator. It combines two multiplicative + * \linear_congruential number generators, i.e. those with @c c = 0. + * It is described in + * + * @blockquote + * "Efficient and Portable Combined Random Number Generators", Pierre L'Ecuyer, + * Communications of the ACM, Vol. 31, No. 6, June 1988, pp. 742-749, 774 + * @endblockquote + * + * The template parameters MLCG1 and MLCG2 shall denote two different + * \linear_congruential number generators, each with c = 0. Each invocation + * returns a random number X(n) := (MLCG1(n) - MLCG2(n)) mod (m1 - 1), where + * m1 denotes the modulus of MLCG1. + * + * The template parameter @c val is the validation value checked by validation. + */ +template +class additive_combine +{ +public: + typedef MLCG1 first_base; + typedef MLCG2 second_base; + typedef typename MLCG1::result_type result_type; +#ifndef BOOST_NO_INCLASS_MEMBER_INITIALIZATION + static const bool has_fixed_range = true; + static const result_type min_value = 1; + static const result_type max_value = MLCG1::max_value-1; +#else + enum { has_fixed_range = false }; +#endif + /** + * Returns: The smallest value that the generator can produce + */ + result_type min BOOST_PREVENT_MACRO_SUBSTITUTION () const { return 1; } + /** + * Returns: The largest value that the generator can produce + */ + result_type max BOOST_PREVENT_MACRO_SUBSTITUTION () const { return (_mlcg1.max)()-1; } + + /** + * Constructs an \additive_combine generator using the + * default constructors of the two base generators. + */ + additive_combine() : _mlcg1(), _mlcg2() { } + /** + * Constructs an \additive_combine generator, using aseed as + * the constructor argument for both base generators. + */ + explicit additive_combine(result_type aseed) + : _mlcg1(aseed), _mlcg2(aseed) { } + /** + * Constructs an \additive_combine generator, using + * @c seed1 and @c seed2 as the constructor argument to + * the first and second base generators, respectively. + */ + additive_combine(typename MLCG1::result_type seed1, + typename MLCG2::result_type seed2) + : _mlcg1(seed1), _mlcg2(seed2) { } + /** + * Contructs an \additive_combine generator with + * values from the range defined by the input iterators first + * and last. first will be modified to point to the element + * after the last one used. + * + * Throws: @c std::invalid_argument if the input range is too small. + * + * Exception Safety: Basic + */ + template additive_combine(It& first, It last) + : _mlcg1(first, last), _mlcg2(first, last) { } + + /** + * Seeds an \additive_combine generator using the default + * seeds of the two base generators. + */ + void seed() + { + _mlcg1.seed(); + _mlcg2.seed(); + } + + /** + * Seeds an \additive_combine generator, using @c aseed as the + * seed for both base generators. + */ + void seed(result_type aseed) + { + _mlcg1.seed(aseed); + _mlcg2.seed(aseed); + } + + /** + * Seeds an \additive_combine generator, using @c seed1 and @c seed2 as + * the seeds to the first and second base generators, respectively. + */ + void seed(typename MLCG1::result_type seed1, + typename MLCG2::result_type seed2) + { + _mlcg1.seed(seed1); + _mlcg2.seed(seed2); + } + + /** + * Seeds an \additive_combine generator with + * values from the range defined by the input iterators first + * and last. first will be modified to point to the element + * after the last one used. + * + * Throws: @c std::invalid_argument if the input range is too small. + * + * Exception Safety: Basic + */ + template void seed(It& first, It last) + { + _mlcg1.seed(first, last); + _mlcg2.seed(first, last); + } + + /** + * Returns: the next value of the generator + */ + result_type operator()() { + result_type z = _mlcg1() - _mlcg2(); + if(z < 1) + z += MLCG1::modulus-1; + return z; + } + + static bool validation(result_type x) { return val == x; } + +#ifndef BOOST_NO_OPERATORS_IN_NAMESPACE + +#ifndef BOOST_RANDOM_NO_STREAM_OPERATORS + /** + * Writes the state of an \additive_combine generator to a @c + * std::ostream. The textual representation of an \additive_combine + * generator is the textual representation of the first base + * generator followed by the textual representation of the + * second base generator. + */ + template + friend std::basic_ostream& + operator<<(std::basic_ostream& os, const additive_combine& r) + { os << r._mlcg1 << " " << r._mlcg2; return os; } + + /** + * Reads the state of an \additive_combine generator from a + * @c std::istream. + */ + template + friend std::basic_istream& + operator>>(std::basic_istream& is, additive_combine& r) + { is >> r._mlcg1 >> std::ws >> r._mlcg2; return is; } +#endif + + /** + * Returns: true iff the two \additive_combine generators will + * produce the same sequence of values. + */ + friend bool operator==(const additive_combine& x, const additive_combine& y) + { return x._mlcg1 == y._mlcg1 && x._mlcg2 == y._mlcg2; } + /** + * Returns: true iff the two \additive_combine generators will + * produce different sequences of values. + */ + friend bool operator!=(const additive_combine& x, const additive_combine& y) + { return !(x == y); } +#else + // Use a member function; Streamable concept not supported. + bool operator==(const additive_combine& rhs) const + { return _mlcg1 == rhs._mlcg1 && _mlcg2 == rhs._mlcg2; } + bool operator!=(const additive_combine& rhs) const + { return !(*this == rhs); } +#endif + +private: + MLCG1 _mlcg1; + MLCG2 _mlcg2; +}; + +} // namespace random + +/** + * The specialization \ecuyer1988 was suggested in + * + * @blockquote + * "Efficient and Portable Combined Random Number Generators", Pierre L'Ecuyer, + * Communications of the ACM, Vol. 31, No. 6, June 1988, pp. 742-749, 774 + * @endblockquote + */ +typedef random::additive_combine< + random::linear_congruential, + random::linear_congruential, + 2060321752> ecuyer1988; + +} // namespace boost + +#endif // BOOST_RANDOM_ADDITIVE_COMBINE_HPP diff --git a/Boost/boost/random/bernoulli_distribution.hpp b/Boost/boost/random/bernoulli_distribution.hpp new file mode 100644 index 000000000..8b4276314 --- /dev/null +++ b/Boost/boost/random/bernoulli_distribution.hpp @@ -0,0 +1,109 @@ +/* boost random/bernoulli_distribution.hpp header file + * + * Copyright Jens Maurer 2000-2001 + * 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: bernoulli_distribution.hpp 60755 2010-03-22 00:45:06Z steven_watanabe $ + * + * Revision history + * 2001-02-18 moved to individual header files + */ + +#ifndef BOOST_RANDOM_BERNOULLI_DISTRIBUTION_HPP +#define BOOST_RANDOM_BERNOULLI_DISTRIBUTION_HPP + +#include +#include +#include + +namespace boost { + +/** + * Instantiations of class template \bernoulli_distribution model a + * \random_distribution. Such a random distribution produces bool values + * distributed with probabilities P(true) = p and P(false) = 1-p. p is + * the parameter of the distribution. + */ +template +class bernoulli_distribution +{ +public: + // In principle, this could work with both integer and floating-point + // types. Generating floating-point random numbers in the first + // place is probably more expensive, so use integer as input. + typedef int input_type; + typedef bool result_type; + + /** + * Constructs a \bernoulli_distribution object. + * p is the parameter of the distribution. + * + * Requires: 0 <= p <= 1 + */ + explicit bernoulli_distribution(const RealType& p_arg = RealType(0.5)) + : _p(p_arg) + { + assert(_p >= 0); + assert(_p <= 1); + } + + // compiler-generated copy ctor and assignment operator are fine + + /** + * Returns: The "p" parameter of the distribution. + */ + RealType p() const { return _p; } + /** + * Effects: Subsequent uses of the distribution do not depend + * on values produced by any engine prior to invoking reset. + */ + void reset() { } + + /** + * Returns: a random variate distributed according to the + * \bernoulli_distribution. + */ + template + result_type operator()(Engine& eng) + { + if(_p == RealType(0)) + return false; + else + return RealType(eng() - (eng.min)()) <= _p * RealType((eng.max)()-(eng.min)()); + } + +#ifndef BOOST_RANDOM_NO_STREAM_OPERATORS + /** + * Writes the parameters of the distribution to a @c std::ostream. + */ + template + friend std::basic_ostream& + operator<<(std::basic_ostream& os, const bernoulli_distribution& bd) + { + os << bd._p; + return os; + } + + /** + * Reads the parameters of the distribution from a @c std::istream. + */ + template + friend std::basic_istream& + operator>>(std::basic_istream& is, bernoulli_distribution& bd) + { + is >> std::ws >> bd._p; + return is; + } +#endif + +private: + RealType _p; +}; + +} // namespace boost + +#endif // BOOST_RANDOM_BERNOULLI_DISTRIBUTION_HPP diff --git a/Boost/boost/random/binomial_distribution.hpp b/Boost/boost/random/binomial_distribution.hpp new file mode 100644 index 000000000..1ca110d08 --- /dev/null +++ b/Boost/boost/random/binomial_distribution.hpp @@ -0,0 +1,111 @@ +/* boost random/binomial_distribution.hpp header file + * + * Copyright Jens Maurer 2002 + * 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: binomial_distribution.hpp 60755 2010-03-22 00:45:06Z steven_watanabe $ + * + */ + +#ifndef BOOST_RANDOM_BINOMIAL_DISTRIBUTION_HPP +#define BOOST_RANDOM_BINOMIAL_DISTRIBUTION_HPP + +#include +#include +#include +#include + +namespace boost { + +/** + * The binomial distribution is an integer valued distribution with + * two parameters, @c t and @c p. The values of the distribution + * are within the range [0,t]. + * + * The probability that the distribution produces a value k is + * \f${t \choose k}p^k(1-p)^{t-k}\f$. + */ +template +class binomial_distribution +{ +public: + typedef typename bernoulli_distribution::input_type input_type; + typedef IntType result_type; + + /** + * Construct an @c binomial_distribution object. @c t and @c p + * are the parameters of the distribution. + * + * Requires: t >=0 && 0 <= p <= 1 + */ + explicit binomial_distribution(IntType t = 1, + const RealType& p = RealType(0.5)) + : _bernoulli(p), _t(t) + { + assert(_t >= 0); + assert(RealType(0) <= p && p <= RealType(1)); + } + + // compiler-generated copy ctor and assignment operator are fine + + /** Returns: the @c t parameter of the distribution */ + IntType t() const { return _t; } + /** Returns: the @c p parameter of the distribution */ + RealType p() const { return _bernoulli.p(); } + /** + * Effects: Subsequent uses of the distribution do not depend + * on values produced by any engine prior to invoking reset. + */ + void reset() { } + + /** + * Returns: a random variate distributed according to the + * binomial distribution. + */ + template + result_type operator()(Engine& eng) + { + // TODO: This is O(_t), but it should be O(log(_t)) for large _t + result_type n = 0; + for(IntType i = 0; i < _t; ++i) + if(_bernoulli(eng)) + ++n; + return n; + } + +#ifndef BOOST_RANDOM_NO_STREAM_OPERATORS + /** + * Writes the parameters of the distribution to a @c std::ostream. + */ + template + friend std::basic_ostream& + operator<<(std::basic_ostream& os, const binomial_distribution& bd) + { + os << bd._bernoulli << " " << bd._t; + return os; + } + + /** + * Reads the parameters of the distribution from a @c std::istream. + */ + template + friend std::basic_istream& + operator>>(std::basic_istream& is, binomial_distribution& bd) + { + is >> std::ws >> bd._bernoulli >> std::ws >> bd._t; + return is; + } +#endif + +private: + bernoulli_distribution _bernoulli; + IntType _t; +}; + +} // namespace boost + +#endif // BOOST_RANDOM_BINOMIAL_DISTRIBUTION_HPP diff --git a/Boost/boost/random/cauchy_distribution.hpp b/Boost/boost/random/cauchy_distribution.hpp new file mode 100644 index 000000000..8555104d1 --- /dev/null +++ b/Boost/boost/random/cauchy_distribution.hpp @@ -0,0 +1,121 @@ +/* boost random/cauchy_distribution.hpp header file + * + * Copyright Jens Maurer 2000-2001 + * 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: cauchy_distribution.hpp 60755 2010-03-22 00:45:06Z steven_watanabe $ + * + * Revision history + * 2001-02-18 moved to individual header files + */ + +#ifndef BOOST_RANDOM_CAUCHY_DISTRIBUTION_HPP +#define BOOST_RANDOM_CAUCHY_DISTRIBUTION_HPP + +#include +#include +#include +#include +#include + +namespace boost { + +#if defined(__GNUC__) && (__GNUC__ < 3) +// Special gcc workaround: gcc 2.95.x ignores using-declarations +// in template classes (confirmed by gcc author Martin v. Loewis) + using std::tan; +#endif + +// Cauchy distribution: + +/** + * The cauchy distribution is a continuous distribution with two + * parameters, sigma and median. + * + * It has \f$p(x) = \frac{\sigma}{\pi(\sigma^2 + (x-m)^2)}\f$ + */ +template +class cauchy_distribution +{ +public: + typedef RealType input_type; + typedef RealType result_type; + +#ifndef BOOST_NO_LIMITS_COMPILE_TIME_CONSTANTS + BOOST_STATIC_ASSERT(!std::numeric_limits::is_integer); +#endif + + /** + * Constructs a \cauchy_distribution with the paramters @c median + * and @c sigma. + */ + explicit cauchy_distribution(result_type median_arg = result_type(0), + result_type sigma_arg = result_type(1)) + : _median(median_arg), _sigma(sigma_arg) { } + + // compiler-generated copy ctor and assignment operator are fine + + /** + * Returns: the "median" parameter of the distribution + */ + result_type median() const { return _median; } + /** + * Returns: the "sigma" parameter of the distribution + */ + result_type sigma() const { return _sigma; } + /** + * Effects: Subsequent uses of the distribution do not depend + * on values produced by any engine prior to invoking reset. + */ + void reset() { } + + /** + * Returns: A random variate distributed according to the + * cauchy distribution. + */ + template + result_type operator()(Engine& eng) + { + // Can we have a boost::mathconst please? + const result_type pi = result_type(3.14159265358979323846); +#ifndef BOOST_NO_STDC_NAMESPACE + using std::tan; +#endif + return _median + _sigma * tan(pi*(eng()-result_type(0.5))); + } + +#ifndef BOOST_RANDOM_NO_STREAM_OPERATORS + /** + * Writes the parameters of the distribution to a @c std::ostream. + */ + template + friend std::basic_ostream& + operator<<(std::basic_ostream& os, const cauchy_distribution& cd) + { + os << cd._median << " " << cd._sigma; + return os; + } + + /** + * Reads the parameters of the distribution from a @c std::istream. + */ + template + friend std::basic_istream& + operator>>(std::basic_istream& is, cauchy_distribution& cd) + { + is >> std::ws >> cd._median >> std::ws >> cd._sigma; + return is; + } +#endif + +private: + result_type _median, _sigma; +}; + +} // namespace boost + +#endif // BOOST_RANDOM_CAUCHY_DISTRIBUTION_HPP diff --git a/Boost/boost/random/detail/auto_link.hpp b/Boost/boost/random/detail/auto_link.hpp new file mode 100644 index 000000000..aec0af630 --- /dev/null +++ b/Boost/boost/random/detail/auto_link.hpp @@ -0,0 +1,42 @@ +/* boost random auto_link.hpp header file + * + * Copyright Steven Watanabe 2010 + * 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) + * + * $Id: auto_link.hpp 60755 2010-03-22 00:45:06Z steven_watanabe $ + */ + +#ifndef BOOST_RANDOM_DETAIL_AUTO_LINK_HPP +#define BOOST_RANDOM_DETAIL_AUTO_LINK_HPP + +#include + +#ifdef BOOST_HAS_DECLSPEC + #if defined(BOOST_ALL_DYN_LINK) || defined(BOOST_RANDOM_DYN_LINK) + #if defined(BOOST_RANDOM_SOURCE) + #define BOOST_RANDOM_DECL __declspec(dllexport) + #else + #define BOOST_RANDOM_DECL __declspec(dllimport) + #endif + #endif +#endif + +#ifndef BOOST_RANDOM_DECL + #define BOOST_RANDOM_DECL +#endif + +#if !defined(BOOST_RANDOM_NO_LIB) && !defined(BOOST_ALL_NO_LIB) && !defined(BOOST_RANDOM_SOURCE) + +#define BOOST_LIB_NAME boost_random + +#if defined(BOOST_RANDOM_DYN_LINK) || defined(BOOST_ALL_DYN_LINK) + #define BOOST_DYN_LINK +#endif + +#include + +#endif + +#endif diff --git a/Boost/boost/random/detail/config.hpp b/Boost/boost/random/detail/config.hpp new file mode 100644 index 000000000..bb834534b --- /dev/null +++ b/Boost/boost/random/detail/config.hpp @@ -0,0 +1,18 @@ +/* boost random/detail/config.hpp header file + * + * Copyright Steven Watanabe 2009 + * 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: config.hpp 52492 2009-04-19 14:55:57Z steven_watanabe $ + */ + +#include + +#if (defined(BOOST_NO_OPERATORS_IN_NAMESPACE) || defined(BOOST_NO_MEMBER_TEMPLATE_FRIENDS)) \ + && !defined(BOOST_MSVC) + #define BOOST_RANDOM_NO_STREAM_OPERATORS +#endif diff --git a/Boost/boost/random/detail/const_mod.hpp b/Boost/boost/random/detail/const_mod.hpp new file mode 100644 index 000000000..07f754506 --- /dev/null +++ b/Boost/boost/random/detail/const_mod.hpp @@ -0,0 +1,363 @@ +/* boost random/detail/const_mod.hpp header file + * + * Copyright Jens Maurer 2000-2001 + * 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: const_mod.hpp 58649 2010-01-02 21:23:17Z steven_watanabe $ + * + * Revision history + * 2001-02-18 moved to individual header files + */ + +#ifndef BOOST_RANDOM_CONST_MOD_HPP +#define BOOST_RANDOM_CONST_MOD_HPP + +#include +#include +#include +#include +#include + +#include + +namespace boost { +namespace random { + +/* + * Some random number generators require modular arithmetic. Put + * everything we need here. + * IntType must be an integral type. + */ + +namespace detail { + + template + struct do_add + { }; + + template<> + struct do_add + { + template + static IntType add(IntType m, IntType x, IntType c) + { + if (x < m - c) + return x + c; + else + return x - (m-c); + } + }; + + template<> + struct do_add + { + template + static IntType add(IntType, IntType, IntType) + { + // difficult + assert(!"const_mod::add with c too large"); + return 0; + } + }; +} // namespace detail + +#if !(defined(__BORLANDC__) && (__BORLANDC__ == 0x560)) + +template +class const_mod +{ +public: + static IntType add(IntType x, IntType c) + { + if(c == 0) + return x; + else if(c <= traits::const_max - m) // i.e. m+c < max + return add_small(x, c); + else + return detail::do_add::add(m, x, c); + } + + static IntType mult(IntType a, IntType x) + { + if(a == 1) + return x; + else if(m <= traits::const_max/a) // i.e. a*m <= max + return mult_small(a, x); + else if(traits::is_signed && (m%a < m/a)) + return mult_schrage(a, x); + else { + // difficult + assert(!"const_mod::mult with a too large"); + return 0; + } + } + + static IntType mult_add(IntType a, IntType x, IntType c) + { + if(m <= (traits::const_max-c)/a) // i.e. a*m+c <= max + return (a*x+c) % m; + else + return add(mult(a, x), c); + } + + static IntType invert(IntType x) + { return x == 0 ? 0 : invert_euclidian(x); } + +private: + typedef integer_traits traits; + + const_mod(); // don't instantiate + + static IntType add_small(IntType x, IntType c) + { + x += c; + if(x >= m) + x -= m; + return x; + } + + static IntType mult_small(IntType a, IntType x) + { + return a*x % m; + } + + static IntType mult_schrage(IntType a, IntType value) + { + const IntType q = m / a; + const IntType r = m % a; + + assert(r < q); // check that overflow cannot happen + + value = a*(value%q) - r*(value/q); + // An optimizer bug in the SGI MIPSpro 7.3.1.x compiler requires this + // convoluted formulation of the loop (Synge Todo) + for(;;) { + if (value > 0) + break; + value += m; + } + return value; + } + + // invert c in the finite field (mod m) (m must be prime) + static IntType invert_euclidian(IntType c) + { + // we are interested in the gcd factor for c, because this is our inverse + BOOST_STATIC_ASSERT(m > 0); +#if BOOST_WORKAROUND(__MWERKS__, BOOST_TESTED_AT(0x3003)) + assert(boost::integer_traits::is_signed); +#elif !defined(BOOST_NO_LIMITS_COMPILE_TIME_CONSTANTS) + BOOST_STATIC_ASSERT(boost::integer_traits::is_signed); +#endif + assert(c > 0); + IntType l1 = 0; + IntType l2 = 1; + IntType n = c; + IntType p = m; + for(;;) { + IntType q = p / n; + l1 -= q * l2; // this requires a signed IntType! + p -= q * n; + if(p == 0) + return (l2 < 1 ? l2 + m : l2); + IntType q2 = n / p; + l2 -= q2 * l1; + n -= q2 * p; + if(n == 0) + return (l1 < 1 ? l1 + m : l1); + } + } +}; + +// The modulus is exactly the word size: rely on machine overflow handling. +// Due to a GCC bug, we cannot partially specialize in the presence of +// template value parameters. +template<> +class const_mod +{ + typedef unsigned int IntType; +public: + static IntType add(IntType x, IntType c) { return x+c; } + static IntType mult(IntType a, IntType x) { return a*x; } + static IntType mult_add(IntType a, IntType x, IntType c) { return a*x+c; } + + // m is not prime, thus invert is not useful +private: // don't instantiate + const_mod(); +}; + +template<> +class const_mod +{ + typedef unsigned long IntType; +public: + static IntType add(IntType x, IntType c) { return x+c; } + static IntType mult(IntType a, IntType x) { return a*x; } + static IntType mult_add(IntType a, IntType x, IntType c) { return a*x+c; } + + // m is not prime, thus invert is not useful +private: // don't instantiate + const_mod(); +}; + +// the modulus is some power of 2: rely partly on machine overflow handling +// we only specialize for rand48 at the moment +#ifndef BOOST_NO_INT64_T +template<> +class const_mod +{ + typedef uint64_t IntType; +public: + static IntType add(IntType x, IntType c) { return c == 0 ? x : mod(x+c); } + static IntType mult(IntType a, IntType x) { return mod(a*x); } + static IntType mult_add(IntType a, IntType x, IntType c) + { return mod(a*x+c); } + static IntType mod(IntType x) { return x &= ((uint64_t(1) << 48)-1); } + + // m is not prime, thus invert is not useful +private: // don't instantiate + const_mod(); +}; +#endif /* !BOOST_NO_INT64_T */ + +#else + +// +// for some reason Borland C++ Builder 6 has problems with +// the full specialisations of const_mod, define a generic version +// instead, the compiler will optimise away the const-if statements: +// + +template +class const_mod +{ +public: + static IntType add(IntType x, IntType c) + { + if(0 == m) + { + return x+c; + } + else + { + if(c == 0) + return x; + else if(c <= traits::const_max - m) // i.e. m+c < max + return add_small(x, c); + else + return detail::do_add::add(m, x, c); + } + } + + static IntType mult(IntType a, IntType x) + { + if(x == 0) + { + return a*x; + } + else + { + if(a == 1) + return x; + else if(m <= traits::const_max/a) // i.e. a*m <= max + return mult_small(a, x); + else if(traits::is_signed && (m%a < m/a)) + return mult_schrage(a, x); + else { + // difficult + assert(!"const_mod::mult with a too large"); + return 0; + } + } + } + + static IntType mult_add(IntType a, IntType x, IntType c) + { + if(m == 0) + { + return a*x+c; + } + else + { + if(m <= (traits::const_max-c)/a) // i.e. a*m+c <= max + return (a*x+c) % m; + else + return add(mult(a, x), c); + } + } + + static IntType invert(IntType x) + { return x == 0 ? 0 : invert_euclidian(x); } + +private: + typedef integer_traits traits; + + const_mod(); // don't instantiate + + static IntType add_small(IntType x, IntType c) + { + x += c; + if(x >= m) + x -= m; + return x; + } + + static IntType mult_small(IntType a, IntType x) + { + return a*x % m; + } + + static IntType mult_schrage(IntType a, IntType value) + { + const IntType q = m / a; + const IntType r = m % a; + + assert(r < q); // check that overflow cannot happen + + value = a*(value%q) - r*(value/q); + while(value <= 0) + value += m; + return value; + } + + // invert c in the finite field (mod m) (m must be prime) + static IntType invert_euclidian(IntType c) + { + // we are interested in the gcd factor for c, because this is our inverse + BOOST_STATIC_ASSERT(m > 0); +#ifndef BOOST_NO_LIMITS_COMPILE_TIME_CONSTANTS + BOOST_STATIC_ASSERT(boost::integer_traits::is_signed); +#endif + assert(c > 0); + IntType l1 = 0; + IntType l2 = 1; + IntType n = c; + IntType p = m; + for(;;) { + IntType q = p / n; + l1 -= q * l2; // this requires a signed IntType! + p -= q * n; + if(p == 0) + return (l2 < 1 ? l2 + m : l2); + IntType q2 = n / p; + l2 -= q2 * l1; + n -= q2 * p; + if(n == 0) + return (l1 < 1 ? l1 + m : l1); + } + } +}; + + +#endif + +} // namespace random +} // namespace boost + +#include + +#endif // BOOST_RANDOM_CONST_MOD_HPP diff --git a/Boost/boost/random/detail/disable_warnings.hpp b/Boost/boost/random/detail/disable_warnings.hpp new file mode 100644 index 000000000..ab129e0b9 --- /dev/null +++ b/Boost/boost/random/detail/disable_warnings.hpp @@ -0,0 +1,23 @@ +/* boost random/detail/disable_warnings.hpp header file + * + * Copyright Steven Watanabe 2009 + * 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: disable_warnings.hpp 60755 2010-03-22 00:45:06Z steven_watanabe $ + * + */ + +// No #include guard. This header is intended to be included multiple times. + +#include + +#ifdef BOOST_MSVC +#pragma warning(push) +#pragma warning(disable:4512) +#pragma warning(disable:4127) +#pragma warning(disable:4724) +#endif diff --git a/Boost/boost/random/detail/enable_warnings.hpp b/Boost/boost/random/detail/enable_warnings.hpp new file mode 100644 index 000000000..12f0840ef --- /dev/null +++ b/Boost/boost/random/detail/enable_warnings.hpp @@ -0,0 +1,18 @@ +/* boost random/detail/enable_warnings.hpp header file + * + * Copyright Steven Watanabe 2009 + * 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: enable_warnings.hpp 58649 2010-01-02 21:23:17Z steven_watanabe $ + * + */ + +// No #include guard. This header is intended to be included multiple times. + +#ifdef BOOST_MSVC +#pragma warning(pop) +#endif diff --git a/Boost/boost/random/detail/iterator_mixin.hpp b/Boost/boost/random/detail/iterator_mixin.hpp new file mode 100644 index 000000000..ca67ebb5a --- /dev/null +++ b/Boost/boost/random/detail/iterator_mixin.hpp @@ -0,0 +1,45 @@ +/* boost random/detail/iterator_mixin.hpp header file + * + * Copyright Jens Maurer 2000-2001 + * 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. + * + * Revision history + */ + +#ifndef BOOST_ITERATOR_MIXIN_HPP +#define BOOST_ITERATOR_MIXIN_HPP + +#include + +namespace boost { + +// must be in boost namespace, otherwise the inline friend trick fails +template +class generator_iterator_mixin_adapter + : incrementable, equality_comparable +{ +public: + typedef std::input_iterator_tag iterator_category; + typedef ResultType value_type; + typedef std::ptrdiff_t difference_type; + typedef const value_type * pointer; + typedef const value_type & reference; + Generator& operator++() { v = cast()(); return cast(); } + const value_type& operator*() const { return v; } + +protected: + // instantiate from derived classes only + generator_iterator_mixin_adapter() { } + void iterator_init() { operator++(); } +private: + Generator & cast() { return static_cast(*this); } + value_type v; +}; + +} // namespace boost + +#endif // BOOST_ITERATOR_MIXIN_HPP diff --git a/Boost/boost/random/detail/pass_through_engine.hpp b/Boost/boost/random/detail/pass_through_engine.hpp new file mode 100644 index 000000000..8627025ef --- /dev/null +++ b/Boost/boost/random/detail/pass_through_engine.hpp @@ -0,0 +1,100 @@ +/* boost random/detail/uniform_int_float.hpp header file + * + * Copyright Jens Maurer 2000-2001 + * 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: pass_through_engine.hpp 58649 2010-01-02 21:23:17Z steven_watanabe $ + * + */ + +#ifndef BOOST_RANDOM_DETAIL_PASS_THROUGH_ENGINE_HPP +#define BOOST_RANDOM_DETAIL_PASS_THROUGH_ENGINE_HPP + +#include +#include +#include + +namespace boost { +namespace random { +namespace detail { + +template +class pass_through_engine +{ +private: + typedef ptr_helper helper_type; + +public: + typedef typename helper_type::value_type base_type; + typedef typename base_type::result_type result_type; + + explicit pass_through_engine(UniformRandomNumberGenerator rng) + // make argument an rvalue to avoid matching Generator& constructor + : _rng(static_cast(rng)) + { } + + result_type min BOOST_PREVENT_MACRO_SUBSTITUTION () const { return (base().min)(); } + result_type max BOOST_PREVENT_MACRO_SUBSTITUTION () const { return (base().max)(); } + base_type& base() { return helper_type::ref(_rng); } + const base_type& base() const { return helper_type::ref(_rng); } + + result_type operator()() { return base()(); } + +private: + UniformRandomNumberGenerator _rng; +}; + +#ifndef BOOST_NO_STD_LOCALE + +template +std::basic_ostream& +operator<<( + std::basic_ostream& os + , const pass_through_engine& ud + ) +{ + return os << ud.base(); +} + +template +std::basic_istream& +operator>>( + std::basic_istream& is + , const pass_through_engine& ud + ) +{ + return is >> ud.base(); +} + +#else // no new streams + +template +inline std::ostream& +operator<<(std::ostream& os, + const pass_through_engine& ud) +{ + return os << ud.base(); +} + +template +inline std::istream& +operator>>(std::istream& is, + const pass_through_engine& ud) +{ + return is >> ud.base(); +} + +#endif + +} // namespace detail +} // namespace random +} // namespace boost + +#include + +#endif // BOOST_RANDOM_DETAIL_PASS_THROUGH_ENGINE_HPP + diff --git a/Boost/boost/random/detail/ptr_helper.hpp b/Boost/boost/random/detail/ptr_helper.hpp new file mode 100644 index 000000000..d583afe8e --- /dev/null +++ b/Boost/boost/random/detail/ptr_helper.hpp @@ -0,0 +1,94 @@ +/* boost random/detail/ptr_helper.hpp header file + * + * Copyright Jens Maurer 2002 + * 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: ptr_helper.hpp 24096 2004-07-27 03:43:34Z dgregor $ + * + */ + +#ifndef BOOST_RANDOM_DETAIL_PTR_HELPER_HPP +#define BOOST_RANDOM_DETAIL_PTR_HELPER_HPP + +#include + + +namespace boost { +namespace random { +namespace detail { + +// type_traits could help here, but I don't want to depend on type_traits. +template +struct ptr_helper +{ + typedef T value_type; + typedef T& reference_type; + typedef const T& rvalue_type; + static reference_type ref(T& r) { return r; } + static const T& ref(const T& r) { return r; } +}; + +#ifndef BOOST_NO_TEMPLATE_PARTIAL_SPECIALIZATION +template +struct ptr_helper +{ + typedef T value_type; + typedef T& reference_type; + typedef T& rvalue_type; + static reference_type ref(T& r) { return r; } + static const T& ref(const T& r) { return r; } +}; + +template +struct ptr_helper +{ + typedef T value_type; + typedef T& reference_type; + typedef T* rvalue_type; + static reference_type ref(T * p) { return *p; } + static const T& ref(const T * p) { return *p; } +}; +#endif + +} // namespace detail +} // namespace random +} // namespace boost + +// +// BOOST_RANDOM_PTR_HELPER_SPEC -- +// +// Helper macro for broken compilers defines specializations of +// ptr_helper. +// +#ifdef BOOST_NO_TEMPLATE_PARTIAL_SPECIALIZATION +# define BOOST_RANDOM_PTR_HELPER_SPEC(T) \ +namespace boost { namespace random { namespace detail { \ +template<> \ +struct ptr_helper \ +{ \ + typedef T value_type; \ + typedef T& reference_type; \ + typedef T& rvalue_type; \ + static reference_type ref(T& r) { return r; } \ + static const T& ref(const T& r) { return r; } \ +}; \ + \ +template<> \ +struct ptr_helper \ +{ \ + typedef T value_type; \ + typedef T& reference_type; \ + typedef T* rvalue_type; \ + static reference_type ref(T * p) { return *p; } \ + static const T& ref(const T * p) { return *p; } \ +}; \ +}}} +#else +# define BOOST_RANDOM_PTR_HELPER_SPEC(T) +#endif + +#endif // BOOST_RANDOM_DETAIL_PTR_HELPER_HPP diff --git a/Boost/boost/random/detail/seed.hpp b/Boost/boost/random/detail/seed.hpp new file mode 100644 index 000000000..9b0f744b8 --- /dev/null +++ b/Boost/boost/random/detail/seed.hpp @@ -0,0 +1,89 @@ +/* boost random/detail/seed.hpp header file + * + * Copyright Steven Watanabe 2009 + * 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: seed.hpp 60755 2010-03-22 00:45:06Z steven_watanabe $ + */ + +#ifndef BOOST_RANDOM_DETAIL_SEED_HPP +#define BOOST_RANDOM_DETAIL_SEED_HPP + +#include + +// Sun seems to have trouble with the use of SFINAE for the +// templated constructor. So does Borland. +#if !defined(BOOST_NO_SFINAE) && !defined(__SUNPRO_CC) && !defined(__BORLANDC__) + +#include +#include + +namespace boost { +namespace random { +namespace detail { + +template +struct disable_seed : boost::disable_if > {}; + +template +struct disable_constructor : disable_seed {}; + +template +struct disable_constructor {}; + +#define BOOST_RANDOM_DETAIL_GENERATOR_CONSTRUCTOR(Self, Generator, gen) \ + template \ + explicit Self(Generator& gen, typename ::boost::random::detail::disable_constructor::type* = 0) + +#define BOOST_RANDOM_DETAIL_GENERATOR_SEED(Self, Generator, gen) \ + template \ + void seed(Generator& gen, typename ::boost::random::detail::disable_seed::type* = 0) + +#define BOOST_RANDOM_DETAIL_ARITHMETIC_CONSTRUCTOR(Self, T, x) \ + explicit Self(const T& x) + +#define BOOST_RANDOM_DETAIL_ARITHMETIC_SEED(Self, T, x) \ + void seed(const T& x) + +} +} +} + +#else + +#include +#include + +#define BOOST_RANDOM_DETAIL_GENERATOR_CONSTRUCTOR(Self, Generator, gen) \ + Self(Self& other) { *this = other; } \ + Self(const Self& other) { *this = other; } \ + template \ + explicit Self(Generator& gen) { \ + boost_random_constructor_impl(gen, ::boost::is_arithmetic());\ + } \ + template \ + void boost_random_constructor_impl(Generator& gen, ::boost::mpl::false_) + +#define BOOST_RANDOM_DETAIL_GENERATOR_SEED(Self, Generator, gen) \ + template \ + void seed(Generator& gen) { \ + boost_random_seed_impl(gen, ::boost::is_arithmetic());\ + }\ + template\ + void boost_random_seed_impl(Generator& gen, ::boost::mpl::false_) + +#define BOOST_RANDOM_DETAIL_ARITHMETIC_CONSTRUCTOR(Self, T, x) \ + explicit Self(const T& x) { boost_random_constructor_impl(x, ::boost::mpl::true_()); }\ + void boost_random_constructor_impl(const T& x, ::boost::mpl::true_) + +#define BOOST_RANDOM_DETAIL_ARITHMETIC_SEED(Self, T, x) \ + void seed(const T& x) { boost_random_seed_impl(x, ::boost::mpl::true_()); }\ + void boost_random_seed_impl(const T& x, ::boost::mpl::true_) + +#endif + +#endif diff --git a/Boost/boost/random/detail/signed_unsigned_tools.hpp b/Boost/boost/random/detail/signed_unsigned_tools.hpp new file mode 100644 index 000000000..c1344a3de --- /dev/null +++ b/Boost/boost/random/detail/signed_unsigned_tools.hpp @@ -0,0 +1,89 @@ +/* boost random/detail/signed_unsigned_tools.hpp header file + * + * Copyright Jens Maurer 2006 + * 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. + */ + +#ifndef BOOST_RANDOM_DETAIL_SIGNED_UNSIGNED_TOOLS +#define BOOST_RANDOM_DETAIL_SIGNED_UNSIGNED_TOOLS + +#include +#include +#include + +namespace boost { +namespace random { +namespace detail { + + +/* + * Compute x - y, we know that x >= y, return an unsigned value. + */ + +template::is_signed> +struct subtract { }; + +template +struct subtract +{ + typedef T result_type; + result_type operator()(T x, T y) { return x - y; } +}; + +template +struct subtract +{ + typedef typename make_unsigned::type result_type; + result_type operator()(T x, T y) + { + if (y >= 0) // because x >= y, it follows that x >= 0, too + return result_type(x) - result_type(y); + if (x >= 0) // y < 0 + // avoid the nasty two's complement case for y == min() + return result_type(x) + result_type(-(y+1)) + 1; + // both x and y are negative: no signed overflow + return result_type(x - y); + } +}; + +/* + * Compute x + y, x is unsigned, result fits in type of "y". + */ + +template::is_signed> +struct add { }; + +template +struct add +{ + typedef T2 result_type; + result_type operator()(T1 x, T2 y) { return T2(x) + y; } +}; + +template +struct add +{ + typedef T2 result_type; + result_type operator()(T1 x, T2 y) + { + if (y >= 0) + return T2(x) + y; + // y < 0 + if (x >= T1(-(y+1))) // result >= 0 after subtraction + // avoid the nasty two's complement edge case for y == min() + return T2(x - T1(-(y+1)) - 1); + // abs(x) < abs(y), thus T2 able to represent x + return T2(x) + y; + } +}; + +} // namespace detail +} // namespace random +} // namespace boost + +#endif // BOOST_RANDOM_DETAIL_SIGNED_UNSIGNED_TOOLS + diff --git a/Boost/boost/random/detail/uniform_int_float.hpp b/Boost/boost/random/detail/uniform_int_float.hpp new file mode 100644 index 000000000..8475f1bf7 --- /dev/null +++ b/Boost/boost/random/detail/uniform_int_float.hpp @@ -0,0 +1,85 @@ +/* boost random/detail/uniform_int_float.hpp header file + * + * Copyright Jens Maurer 2000-2001 + * 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: uniform_int_float.hpp 52492 2009-04-19 14:55:57Z steven_watanabe $ + * + */ + +#ifndef BOOST_RANDOM_DETAIL_UNIFORM_INT_FLOAT_HPP +#define BOOST_RANDOM_DETAIL_UNIFORM_INT_FLOAT_HPP + +#include +#include +#include + + +namespace boost { +namespace random { +namespace detail { + +template +class uniform_int_float +{ +public: + typedef UniformRandomNumberGenerator base_type; + typedef IntType result_type; + + uniform_int_float(base_type rng, IntType min_arg = 0, IntType max_arg = 0xffffffff) + : _rng(rng), _min(min_arg), _max(max_arg) + { + init(); + } + + result_type min BOOST_PREVENT_MACRO_SUBSTITUTION () const { return _min; } + result_type max BOOST_PREVENT_MACRO_SUBSTITUTION () const { return _max; } + base_type& base() { return _rng.base(); } + const base_type& base() const { return _rng.base(); } + + result_type operator()() + { + return static_cast(_rng() * _range) + _min; + } + +#ifndef BOOST_RANDOM_NO_STREAM_OPERATORS + template + friend std::basic_ostream& + operator<<(std::basic_ostream& os, const uniform_int_float& ud) + { + os << ud._min << " " << ud._max; + return os; + } + + template + friend std::basic_istream& + operator>>(std::basic_istream& is, uniform_int_float& ud) + { + is >> std::ws >> ud._min >> std::ws >> ud._max; + ud.init(); + return is; + } +#endif + +private: + void init() + { + _range = static_cast(_max-_min)+1; + } + + typedef typename base_type::result_type base_result; + uniform_01 _rng; + result_type _min, _max; + base_result _range; +}; + + +} // namespace detail +} // namespace random +} // namespace boost + +#endif // BOOST_RANDOM_DETAIL_UNIFORM_INT_FLOAT_HPP diff --git a/Boost/boost/random/discard_block.hpp b/Boost/boost/random/discard_block.hpp new file mode 100644 index 000000000..313363d68 --- /dev/null +++ b/Boost/boost/random/discard_block.hpp @@ -0,0 +1,132 @@ +/* boost random/discard_block.hpp header file + * + * Copyright Jens Maurer 2002 + * 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: discard_block.hpp 60755 2010-03-22 00:45:06Z steven_watanabe $ + * + * Revision history + * 2001-03-02 created + */ + +#ifndef BOOST_RANDOM_DISCARD_BLOCK_HPP +#define BOOST_RANDOM_DISCARD_BLOCK_HPP + +#include +#include +#include +#include +#include + + +namespace boost { +namespace random { + +/** + * The class template \discard_block is a model of + * \pseudo_random_number_generator. It modifies + * another generator by discarding parts of its output. + * Out of every block of @c r results, the first @c p + * will be returned and the rest discarded. + * + * Requires: 0 < p <= r + */ +template +class discard_block +{ +public: + typedef UniformRandomNumberGenerator base_type; + typedef typename base_type::result_type result_type; + + BOOST_STATIC_CONSTANT(bool, has_fixed_range = false); + BOOST_STATIC_CONSTANT(unsigned int, total_block = p); + BOOST_STATIC_CONSTANT(unsigned int, returned_block = r); + +#ifndef BOOST_NO_LIMITS_COMPILE_TIME_CONSTANTS + BOOST_STATIC_ASSERT(total_block >= returned_block); +#endif + + discard_block() : _rng(), _n(0) { } + explicit discard_block(const base_type & rng) : _rng(rng), _n(0) { } + template explicit discard_block(T s) : _rng(s), _n(0) {} + template discard_block(It& first, It last) + : _rng(first, last), _n(0) { } + void seed() { _rng.seed(); _n = 0; } + template void seed(T s) { _rng.seed(s); _n = 0; } + template void seed(It& first, It last) + { _n = 0; _rng.seed(first, last); } + + const base_type& base() const { return _rng; } + + result_type operator()() + { + if(_n >= returned_block) { + // discard values of random number generator + for( ; _n < total_block; ++_n) + _rng(); + _n = 0; + } + ++_n; + return _rng(); + } + + result_type min BOOST_PREVENT_MACRO_SUBSTITUTION () const { return (_rng.min)(); } + result_type max BOOST_PREVENT_MACRO_SUBSTITUTION () const { return (_rng.max)(); } + static bool validation(result_type x) { return true; } // dummy + +#ifndef BOOST_NO_OPERATORS_IN_NAMESPACE + +#ifndef BOOST_RANDOM_NO_STREAM_OPERATORS + template + friend std::basic_ostream& + operator<<(std::basic_ostream& os, const discard_block& s) + { + os << s._rng << " " << s._n << " "; + return os; + } + + template + friend std::basic_istream& + operator>>(std::basic_istream& is, discard_block& s) + { + is >> s._rng >> std::ws >> s._n >> std::ws; + return is; + } +#endif + + friend bool operator==(const discard_block& x, const discard_block& y) + { return x._rng == y._rng && x._n == y._n; } + friend bool operator!=(const discard_block& x, const discard_block& y) + { return !(x == y); } +#else + // Use a member function; Streamable concept not supported. + bool operator==(const discard_block& rhs) const + { return _rng == rhs._rng && _n == rhs._n; } + bool operator!=(const discard_block& rhs) const + { return !(*this == rhs); } +#endif + +private: + base_type _rng; + unsigned int _n; +}; + +#ifndef BOOST_NO_INCLASS_MEMBER_INITIALIZATION +// A definition is required even for integral static constants +template +const bool discard_block::has_fixed_range; +template +const unsigned int discard_block::total_block; +template +const unsigned int discard_block::returned_block; +#endif + +} // namespace random + +} // namespace boost + +#endif // BOOST_RANDOM_DISCARD_BLOCK_HPP diff --git a/Boost/boost/random/exponential_distribution.hpp b/Boost/boost/random/exponential_distribution.hpp new file mode 100644 index 000000000..63abb0402 --- /dev/null +++ b/Boost/boost/random/exponential_distribution.hpp @@ -0,0 +1,86 @@ +/* boost random/exponential_distribution.hpp header file + * + * Copyright Jens Maurer 2000-2001 + * 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: exponential_distribution.hpp 60755 2010-03-22 00:45:06Z steven_watanabe $ + * + * Revision history + * 2001-02-18 moved to individual header files + */ + +#ifndef BOOST_RANDOM_EXPONENTIAL_DISTRIBUTION_HPP +#define BOOST_RANDOM_EXPONENTIAL_DISTRIBUTION_HPP + +#include +#include +#include +#include +#include +#include + +namespace boost { + +/** + * The exponential distribution has a single parameter lambda. + * + * It has \f$p(x) = \lambda e^{-\lambda x}\f$ + */ +template +class exponential_distribution +{ +public: + typedef RealType input_type; + typedef RealType result_type; + +#if !defined(BOOST_NO_LIMITS_COMPILE_TIME_CONSTANTS) && !(defined(BOOST_MSVC) && BOOST_MSVC <= 1300) + BOOST_STATIC_ASSERT(!std::numeric_limits::is_integer); +#endif + + explicit exponential_distribution(result_type lambda_arg = result_type(1)) + : _lambda(lambda_arg) { assert(_lambda > result_type(0)); } + + // compiler-generated copy ctor and assignment operator are fine + + result_type lambda() const { return _lambda; } + + void reset() { } + + template + result_type operator()(Engine& eng) + { +#ifndef BOOST_NO_STDC_NAMESPACE + using std::log; +#endif + return -result_type(1) / _lambda * log(result_type(1)-eng()); + } + +#ifndef BOOST_RANDOM_NO_STREAM_OPERATORS + template + friend std::basic_ostream& + operator<<(std::basic_ostream& os, const exponential_distribution& ed) + { + os << ed._lambda; + return os; + } + + template + friend std::basic_istream& + operator>>(std::basic_istream& is, exponential_distribution& ed) + { + is >> std::ws >> ed._lambda; + return is; + } +#endif + +private: + result_type _lambda; +}; + +} // namespace boost + +#endif // BOOST_RANDOM_EXPONENTIAL_DISTRIBUTION_HPP diff --git a/Boost/boost/random/gamma_distribution.hpp b/Boost/boost/random/gamma_distribution.hpp new file mode 100644 index 000000000..c211286ee --- /dev/null +++ b/Boost/boost/random/gamma_distribution.hpp @@ -0,0 +1,143 @@ +/* boost random/gamma_distribution.hpp header file + * + * Copyright Jens Maurer 2002 + * 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: gamma_distribution.hpp 60755 2010-03-22 00:45:06Z steven_watanabe $ + * + */ + +#ifndef BOOST_RANDOM_GAMMA_DISTRIBUTION_HPP +#define BOOST_RANDOM_GAMMA_DISTRIBUTION_HPP + +#include +#include +#include +#include +#include +#include + +namespace boost { + +// The algorithm is taken from Knuth + +/** + * The gamma distribution is a continuous distribution with a single + * parameter alpha. + * + * It has \f$p(x) = x^{\alpha-1}\frac{e^{-x}}{\Gamma(\alpha)}\f$. + */ +template +class gamma_distribution +{ +public: + typedef RealType input_type; + typedef RealType result_type; + +#ifndef BOOST_NO_LIMITS_COMPILE_TIME_CONSTANTS + BOOST_STATIC_ASSERT(!std::numeric_limits::is_integer); +#endif + + explicit gamma_distribution(const result_type& alpha_arg = result_type(1)) + : _exp(result_type(1)), _alpha(alpha_arg) + { + assert(_alpha > result_type(0)); + init(); + } + + // compiler-generated copy ctor and assignment operator are fine + + RealType alpha() const { return _alpha; } + + void reset() { _exp.reset(); } + + template + result_type operator()(Engine& eng) + { +#ifndef BOOST_NO_STDC_NAMESPACE + // allow for Koenig lookup + using std::tan; using std::sqrt; using std::exp; using std::log; + using std::pow; +#endif + if(_alpha == result_type(1)) { + return _exp(eng); + } else if(_alpha > result_type(1)) { + // Can we have a boost::mathconst please? + const result_type pi = result_type(3.14159265358979323846); + for(;;) { + result_type y = tan(pi * eng()); + result_type x = sqrt(result_type(2)*_alpha-result_type(1))*y + + _alpha-result_type(1); + if(x <= result_type(0)) + continue; + if(eng() > + (result_type(1)+y*y) * exp((_alpha-result_type(1)) + *log(x/(_alpha-result_type(1))) + - sqrt(result_type(2)*_alpha + -result_type(1))*y)) + continue; + return x; + } + } else /* alpha < 1.0 */ { + for(;;) { + result_type u = eng(); + result_type y = _exp(eng); + result_type x, q; + if(u < _p) { + x = exp(-y/_alpha); + q = _p*exp(-x); + } else { + x = result_type(1)+y; + q = _p + (result_type(1)-_p) * pow(x, _alpha-result_type(1)); + } + if(u >= q) + continue; + return x; + } + } + } + +#ifndef BOOST_RANDOM_NO_STREAM_OPERATORS + template + friend std::basic_ostream& + operator<<(std::basic_ostream& os, const gamma_distribution& gd) + { + os << gd._alpha; + return os; + } + + template + friend std::basic_istream& + operator>>(std::basic_istream& is, gamma_distribution& gd) + { + is >> std::ws >> gd._alpha; + gd.init(); + return is; + } +#endif + +private: + /// \cond hide_private_members + void init() + { +#ifndef BOOST_NO_STDC_NAMESPACE + // allow for Koenig lookup + using std::exp; +#endif + _p = exp(result_type(1)) / (_alpha + exp(result_type(1))); + } + /// \endcond + + exponential_distribution _exp; + result_type _alpha; + // some data precomputed from the parameters + result_type _p; +}; + +} // namespace boost + +#endif // BOOST_RANDOM_GAMMA_DISTRIBUTION_HPP diff --git a/Boost/boost/random/geometric_distribution.hpp b/Boost/boost/random/geometric_distribution.hpp new file mode 100644 index 000000000..6f34287fd --- /dev/null +++ b/Boost/boost/random/geometric_distribution.hpp @@ -0,0 +1,118 @@ +/* boost random/geometric_distribution.hpp header file + * + * Copyright Jens Maurer 2000-2001 + * 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: geometric_distribution.hpp 60755 2010-03-22 00:45:06Z steven_watanabe $ + * + * Revision history + * 2001-02-18 moved to individual header files + */ + +#ifndef BOOST_RANDOM_GEOMETRIC_DISTRIBUTION_HPP +#define BOOST_RANDOM_GEOMETRIC_DISTRIBUTION_HPP + +#include // std::log +#include +#include +#include +#include + +namespace boost { + +#if defined(__GNUC__) && (__GNUC__ < 3) +// Special gcc workaround: gcc 2.95.x ignores using-declarations +// in template classes (confirmed by gcc author Martin v. Loewis) + using std::log; +#endif + +/** + * An instantiation of the class template @c geometric_distribution models + * a \random_distribution. The distribution produces positive + * integers which are the number of bernoulli trials + * with probability @c p required to get one that fails. + * + * For the geometric distribution, \f$p(i) = (1-p) p^{i-1}\f$. + */ +template +class geometric_distribution +{ +public: + typedef RealType input_type; + typedef IntType result_type; + + /** + * Contructs a new geometric_distribution with the paramter @c p. + * + * Requires: 0 < p < 1 + */ + explicit geometric_distribution(const RealType& p = RealType(0.5)) + : _p(p) + { + assert(RealType(0) < _p && _p < RealType(1)); + init(); + } + + // compiler-generated copy ctor and assignment operator are fine + + /** + * Returns: the distribution parameter @c p + */ + RealType p() const { return _p; } + void reset() { } + + template + result_type operator()(Engine& eng) + { +#ifndef BOOST_NO_STDC_NAMESPACE + using std::log; + using std::floor; +#endif + return IntType(floor(log(RealType(1)-eng()) / _log_p)) + IntType(1); + } + +#ifndef BOOST_RANDOM_NO_STREAM_OPERATORS + template + friend std::basic_ostream& + operator<<(std::basic_ostream& os, const geometric_distribution& gd) + { + os << gd._p; + return os; + } + + template + friend std::basic_istream& + operator>>(std::basic_istream& is, geometric_distribution& gd) + { + is >> std::ws >> gd._p; + gd.init(); + return is; + } +#endif + +private: + + /// \cond hide_private_functions + + void init() + { +#ifndef BOOST_NO_STDC_NAMESPACE + using std::log; +#endif + _log_p = log(_p); + } + + /// \endcond + + RealType _p; + RealType _log_p; +}; + +} // namespace boost + +#endif // BOOST_RANDOM_GEOMETRIC_DISTRIBUTION_HPP + diff --git a/Boost/boost/random/inversive_congruential.hpp b/Boost/boost/random/inversive_congruential.hpp new file mode 100644 index 000000000..c901ff7e2 --- /dev/null +++ b/Boost/boost/random/inversive_congruential.hpp @@ -0,0 +1,173 @@ +/* boost random/inversive_congruential.hpp header file + * + * Copyright Jens Maurer 2000-2001 + * 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: inversive_congruential.hpp 60755 2010-03-22 00:45:06Z steven_watanabe $ + * + * Revision history + * 2001-02-18 moved to individual header files + */ + +#ifndef BOOST_RANDOM_INVERSIVE_CONGRUENTIAL_HPP +#define BOOST_RANDOM_INVERSIVE_CONGRUENTIAL_HPP + +#include +#include +#include +#include +#include +#include +#include + +namespace boost { +namespace random { + +// Eichenauer and Lehn 1986 +/** + * Instantiations of class template @c inversive_congruential model a + * \pseudo_random_number_generator. It uses the inversive congruential + * algorithm (ICG) described in + * + * @blockquote + * "Inversive pseudorandom number generators: concepts, results and links", + * Peter Hellekalek, In: "Proceedings of the 1995 Winter Simulation + * Conference", C. Alexopoulos, K. Kang, W.R. Lilegdon, and D. Goldsman + * (editors), 1995, pp. 255-262. ftp://random.mat.sbg.ac.at/pub/data/wsc95.ps + * @endblockquote + * + * The output sequence is defined by x(n+1) = (a*inv(x(n)) - b) (mod p), + * where x(0), a, b, and the prime number p are parameters of the generator. + * The expression inv(k) denotes the multiplicative inverse of k in the + * field of integer numbers modulo p, with inv(0) := 0. + * + * The template parameter IntType shall denote a signed integral type large + * enough to hold p; a, b, and p are the parameters of the generators. The + * template parameter val is the validation value checked by validation. + * + * @xmlnote + * The implementation currently uses the Euclidian Algorithm to compute + * the multiplicative inverse. Therefore, the inversive generators are about + * 10-20 times slower than the others (see section"performance"). However, + * the paper talks of only 3x slowdown, so the Euclidian Algorithm is probably + * not optimal for calculating the multiplicative inverse. + * @endxmlnote + */ +template +class inversive_congruential +{ +public: + typedef IntType result_type; +#ifndef BOOST_NO_INCLASS_MEMBER_INITIALIZATION + static const bool has_fixed_range = true; + static const result_type min_value = (b == 0 ? 1 : 0); + static const result_type max_value = p-1; +#else + BOOST_STATIC_CONSTANT(bool, has_fixed_range = false); +#endif + BOOST_STATIC_CONSTANT(result_type, multiplier = a); + BOOST_STATIC_CONSTANT(result_type, increment = b); + BOOST_STATIC_CONSTANT(result_type, modulus = p); + + result_type min BOOST_PREVENT_MACRO_SUBSTITUTION () const { return b == 0 ? 1 : 0; } + result_type max BOOST_PREVENT_MACRO_SUBSTITUTION () const { return p-1; } + + /** + * Constructs an inversive_congruential generator with + * @c y0 as the initial state. + */ + explicit inversive_congruential(IntType y0 = 1) : value(y0) + { + BOOST_STATIC_ASSERT(b >= 0); + BOOST_STATIC_ASSERT(p > 1); + BOOST_STATIC_ASSERT(a >= 1); + if(b == 0) + assert(y0 > 0); + } + template inversive_congruential(It& first, It last) + { seed(first, last); } + + /** Changes the current state to y0. */ + void seed(IntType y0 = 1) { value = y0; if(b == 0) assert(y0 > 0); } + template void seed(It& first, It last) + { + if(first == last) + throw std::invalid_argument("inversive_congruential::seed"); + value = *first++; + } + IntType operator()() + { + typedef const_mod do_mod; + value = do_mod::mult_add(a, do_mod::invert(value), b); + return value; + } + + static bool validation(result_type x) { return val == x; } + +#ifndef BOOST_NO_OPERATORS_IN_NAMESPACE + +#ifndef BOOST_RANDOM_NO_STREAM_OPERATORS + template + friend std::basic_ostream& + operator<<(std::basic_ostream& os, inversive_congruential x) + { os << x.value; return os; } + + template + friend std::basic_istream& + operator>>(std::basic_istream& is, inversive_congruential& x) + { is >> x.value; return is; } +#endif + + friend bool operator==(inversive_congruential x, inversive_congruential y) + { return x.value == y.value; } + friend bool operator!=(inversive_congruential x, inversive_congruential y) + { return !(x == y); } +#else + // Use a member function; Streamable concept not supported. + bool operator==(inversive_congruential rhs) const + { return value == rhs.value; } + bool operator!=(inversive_congruential rhs) const + { return !(*this == rhs); } +#endif +private: + IntType value; +}; + +#ifndef BOOST_NO_INCLASS_MEMBER_INITIALIZATION +// A definition is required even for integral static constants +template +const bool inversive_congruential::has_fixed_range; +template +const typename inversive_congruential::result_type inversive_congruential::min_value; +template +const typename inversive_congruential::result_type inversive_congruential::max_value; +template +const typename inversive_congruential::result_type inversive_congruential::multiplier; +template +const typename inversive_congruential::result_type inversive_congruential::increment; +template +const typename inversive_congruential::result_type inversive_congruential::modulus; +#endif + +} // namespace random + +/** + * The specialization hellekalek1995 was suggested in + * + * @blockquote + * "Inversive pseudorandom number generators: concepts, results and links", + * Peter Hellekalek, In: "Proceedings of the 1995 Winter Simulation + * Conference", C. Alexopoulos, K. Kang, W.R. Lilegdon, and D. Goldsman + * (editors), 1995, pp. 255-262. ftp://random.mat.sbg.ac.at/pub/data/wsc95.ps + * @endblockquote + */ +typedef random::inversive_congruential hellekalek1995; + +} // namespace boost + +#endif // BOOST_RANDOM_INVERSIVE_CONGRUENTIAL_HPP diff --git a/Boost/boost/random/lagged_fibonacci.hpp b/Boost/boost/random/lagged_fibonacci.hpp new file mode 100644 index 000000000..50db47d88 --- /dev/null +++ b/Boost/boost/random/lagged_fibonacci.hpp @@ -0,0 +1,605 @@ +/* boost random/lagged_fibonacci.hpp header file + * + * Copyright Jens Maurer 2000-2001 + * 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: lagged_fibonacci.hpp 60755 2010-03-22 00:45:06Z steven_watanabe $ + * + * Revision history + * 2001-02-18 moved to individual header files + */ + +#ifndef BOOST_RANDOM_LAGGED_FIBONACCI_HPP +#define BOOST_RANDOM_LAGGED_FIBONACCI_HPP + +#include +#include +#include // std::max +#include +#include // std::pow +#include +#include +#include +#include +#include +#include +#include +#include +#include + +namespace boost { +namespace random { + +#if BOOST_WORKAROUND(_MSC_FULL_VER, BOOST_TESTED_AT(13102292)) && BOOST_MSVC > 1300 +# define BOOST_RANDOM_EXTRACT_LF +#endif + +#if defined(__APPLE_CC__) && defined(__GNUC__) && (__GNUC__ == 3) && (__GNUC_MINOR__ <= 3) +# define BOOST_RANDOM_EXTRACT_LF +#endif + +# ifdef BOOST_RANDOM_EXTRACT_LF +namespace detail +{ + template + IStream& + extract_lagged_fibonacci_01( + IStream& is + , F const& f + , unsigned int& i + , RealType* x + , RealType modulus) + { + is >> i >> std::ws; + for(unsigned int i = 0; i < f.long_lag; ++i) + { + RealType value; + is >> value >> std::ws; + x[i] = value / modulus; + } + return is; + } + + template + IStream& + extract_lagged_fibonacci( + IStream& is + , F const& f + , unsigned int& i + , UIntType* x) + { + is >> i >> std::ws; + for(unsigned int i = 0; i < f.long_lag; ++i) + is >> x[i] >> std::ws; + return is; + } +} +# endif + +/** + * Instantiations of class template \lagged_fibonacci model a + * \pseudo_random_number_generator. It uses a lagged Fibonacci + * algorithm with two lags @c p and @c q: + * x(i) = x(i-p) + x(i-q) (mod 2w) with p > q. + */ +template +class lagged_fibonacci +{ +public: + typedef UIntType result_type; + BOOST_STATIC_CONSTANT(bool, has_fixed_range = false); + BOOST_STATIC_CONSTANT(int, word_size = w); + BOOST_STATIC_CONSTANT(unsigned int, long_lag = p); + BOOST_STATIC_CONSTANT(unsigned int, short_lag = q); + + /** + * Returns: the smallest value that the generator can produce + */ + result_type min BOOST_PREVENT_MACRO_SUBSTITUTION () const { return 0; } + /** + * Returns: the largest value that the generator can produce + */ + result_type max BOOST_PREVENT_MACRO_SUBSTITUTION () const { return wordmask; } + + /** + * Creates a new @c lagged_fibonacci generator and calls @c seed() + */ + lagged_fibonacci() { init_wordmask(); seed(); } + /** + * Creates a new @c lagged_fibonacci generator and calls @c seed(value) + */ + explicit lagged_fibonacci(uint32_t value) { init_wordmask(); seed(value); } + /** + * Creates a new @c lagged_fibonacci generator and calls @c seed(first, last) + */ + template lagged_fibonacci(It& first, It last) + { init_wordmask(); seed(first, last); } + // compiler-generated copy ctor and assignment operator are fine + +private: + /// \cond hide_private_members + void init_wordmask() + { + wordmask = 0; + for(int j = 0; j < w; ++j) + wordmask |= (1u << j); + } + /// \endcond + +public: + /** + * Sets the state of the generator to values produced by + * a \minstd_rand generator. + */ + void seed(uint32_t value = 331u) + { + minstd_rand0 gen(value); + for(unsigned int j = 0; j < long_lag; ++j) + x[j] = gen() & wordmask; + i = long_lag; + } + + /** + * Sets the state of the generator to values from the iterator + * range [first, last). If there are not enough elements in the + * range [first, last) throws @c std::invalid_argument. + */ + template + void seed(It& first, It last) + { + // word size could be smaller than the seed values + unsigned int j; + for(j = 0; j < long_lag && first != last; ++j, ++first) + x[j] = *first & wordmask; + i = long_lag; + if(first == last && j < long_lag) + throw std::invalid_argument("lagged_fibonacci::seed"); + } + + /** + * Returns: the next value of the generator + */ + result_type operator()() + { + if(i >= long_lag) + fill(); + return x[i++]; + } + + static bool validation(result_type x) + { + return x == val; + } + +#ifndef BOOST_NO_OPERATORS_IN_NAMESPACE + +#ifndef BOOST_RANDOM_NO_STREAM_OPERATORS + template + friend std::basic_ostream& + operator<<(std::basic_ostream& os, const lagged_fibonacci& f) + { + os << f.i << " "; + for(unsigned int i = 0; i < f.long_lag; ++i) + os << f.x[i] << " "; + return os; + } + + template + friend std::basic_istream& + operator>>(std::basic_istream& is, lagged_fibonacci& f) + { +# ifdef BOOST_RANDOM_EXTRACT_LF + return detail::extract_lagged_fibonacci(is, f, f.i, f.x); +# else + is >> f.i >> std::ws; + for(unsigned int i = 0; i < f.long_lag; ++i) + is >> f.x[i] >> std::ws; + return is; +# endif + } +#endif + + friend bool operator==(const lagged_fibonacci& x, const lagged_fibonacci& y) + { return x.i == y.i && std::equal(x.x, x.x+long_lag, y.x); } + friend bool operator!=(const lagged_fibonacci& x, + const lagged_fibonacci& y) + { return !(x == y); } +#else + // Use a member function; Streamable concept not supported. + bool operator==(const lagged_fibonacci& rhs) const + { return i == rhs.i && std::equal(x, x+long_lag, rhs.x); } + bool operator!=(const lagged_fibonacci& rhs) const + { return !(*this == rhs); } +#endif + +private: + /// \cond hide_private_members + void fill(); + /// \endcond + + UIntType wordmask; + unsigned int i; + UIntType x[long_lag]; +}; + +#ifndef BOOST_NO_INCLASS_MEMBER_INITIALIZATION +// A definition is required even for integral static constants +template +const bool lagged_fibonacci::has_fixed_range; +template +const unsigned int lagged_fibonacci::long_lag; +template +const unsigned int lagged_fibonacci::short_lag; +#endif + +/// \cond hide_private_members + +template +void lagged_fibonacci::fill() +{ + // two loops to avoid costly modulo operations + { // extra scope for MSVC brokenness w.r.t. for scope + for(unsigned int j = 0; j < short_lag; ++j) + x[j] = (x[j] + x[j+(long_lag-short_lag)]) & wordmask; + } + for(unsigned int j = short_lag; j < long_lag; ++j) + x[j] = (x[j] + x[j-short_lag]) & wordmask; + i = 0; +} + + + +// lagged Fibonacci generator for the range [0..1) +// contributed by Matthias Troyer +// for p=55, q=24 originally by G. J. Mitchell and D. P. Moore 1958 + +template +struct fibonacci_validation +{ + BOOST_STATIC_CONSTANT(bool, is_specialized = false); + static T value() { return 0; } + static T tolerance() { return 0; } +}; + +#ifndef BOOST_NO_INCLASS_MEMBER_INITIALIZATION +// A definition is required even for integral static constants +template +const bool fibonacci_validation::is_specialized; +#endif + +#define BOOST_RANDOM_FIBONACCI_VAL(T,P,Q,V,E) \ +template<> \ +struct fibonacci_validation \ +{ \ + BOOST_STATIC_CONSTANT(bool, is_specialized = true); \ + static T value() { return V; } \ + static T tolerance() \ +{ return (std::max)(E, static_cast(5*std::numeric_limits::epsilon())); } \ +}; +// (The extra static_cast in the std::max call above is actually +// unnecessary except for HP aCC 1.30, which claims that +// numeric_limits::epsilon() doesn't actually return a double.) + +BOOST_RANDOM_FIBONACCI_VAL(double, 607, 273, 0.4293817707235914, 1e-14) +BOOST_RANDOM_FIBONACCI_VAL(double, 1279, 418, 0.9421630240437659, 1e-14) +BOOST_RANDOM_FIBONACCI_VAL(double, 2281, 1252, 0.1768114046909004, 1e-14) +BOOST_RANDOM_FIBONACCI_VAL(double, 3217, 576, 0.1956232694868209, 1e-14) +BOOST_RANDOM_FIBONACCI_VAL(double, 4423, 2098, 0.9499762202147172, 1e-14) +BOOST_RANDOM_FIBONACCI_VAL(double, 9689, 5502, 0.05737836943695162, 1e-14) +BOOST_RANDOM_FIBONACCI_VAL(double, 19937, 9842, 0.5076528587449834, 1e-14) +BOOST_RANDOM_FIBONACCI_VAL(double, 23209, 13470, 0.5414473810619185, 1e-14) +BOOST_RANDOM_FIBONACCI_VAL(double, 44497,21034, 0.254135073399297, 1e-14) + +#undef BOOST_RANDOM_FIBONACCI_VAL + +/// \endcond + +/** + * Instantiations of class template @c lagged_fibonacci_01 model a + * \pseudo_random_number_generator. It uses a lagged Fibonacci + * algorithm with two lags @c p and @c q, evaluated in floating-point + * arithmetic: x(i) = x(i-p) + x(i-q) (mod 1) with p > q. See + * + * @blockquote + * "Uniform random number generators for supercomputers", Richard Brent, + * Proc. of Fifth Australian Supercomputer Conference, Melbourne, + * Dec. 1992, pp. 704-706. + * @endblockquote + * + * @xmlnote + * The quality of the generator crucially depends on the choice + * of the parameters. User code should employ one of the sensibly + * parameterized generators such as \lagged_fibonacci607 instead. + * @endxmlnote + * + * The generator requires considerable amounts of memory for the storage + * of its state array. For example, \lagged_fibonacci607 requires about + * 4856 bytes and \lagged_fibonacci44497 requires about 350 KBytes. + */ +template +class lagged_fibonacci_01 +{ +public: + typedef RealType result_type; + BOOST_STATIC_CONSTANT(bool, has_fixed_range = false); + BOOST_STATIC_CONSTANT(int, word_size = w); + BOOST_STATIC_CONSTANT(unsigned int, long_lag = p); + BOOST_STATIC_CONSTANT(unsigned int, short_lag = q); + + /** Constructs a @c lagged_fibonacci_01 generator and calls @c seed(). */ + lagged_fibonacci_01() { init_modulus(); seed(); } + /** Constructs a @c lagged_fibonacci_01 generator and calls @c seed(value). */ + BOOST_RANDOM_DETAIL_ARITHMETIC_CONSTRUCTOR(lagged_fibonacci_01, uint32_t, value) + { init_modulus(); seed(value); } + /** Constructs a @c lagged_fibonacci_01 generator and calls @c seed(gen). */ + BOOST_RANDOM_DETAIL_GENERATOR_CONSTRUCTOR(lagged_fibonacci_01, Generator, gen) + { init_modulus(); seed(gen); } + template lagged_fibonacci_01(It& first, It last) + { init_modulus(); seed(first, last); } + // compiler-generated copy ctor and assignment operator are fine + +private: + /// \cond hide_private_members + void init_modulus() + { +#ifndef BOOST_NO_STDC_NAMESPACE + // allow for Koenig lookup + using std::pow; +#endif + _modulus = pow(RealType(2), word_size); + } + /// \endcond + +public: + /** Calls seed(331u). */ + void seed() { seed(331u); } + /** + * Constructs a \minstd_rand0 generator with the constructor parameter + * value and calls seed with it. Distinct seeds in the range + * [1, 2147483647) will produce generators with different states. Other + * seeds will be equivalent to some seed within this range. See + * \linear_congruential for details. + */ + BOOST_RANDOM_DETAIL_ARITHMETIC_SEED(lagged_fibonacci_01, uint32_t, value) + { + minstd_rand0 intgen(value); + seed(intgen); + } + + /** + * Sets the state of this @c lagged_fibonacci_01 to the values returned + * by p invocations of \uniform_01\()(gen). + * + * Complexity: Exactly p invocations of gen. + */ + BOOST_RANDOM_DETAIL_GENERATOR_SEED(lagged_fibonacci, Generator, gen) + { + // use pass-by-reference, but wrap argument in pass_through_engine + typedef detail::pass_through_engine ref_gen; + uniform_01 gen01 = + uniform_01(ref_gen(gen)); + // I could have used std::generate_n, but it takes "gen" by value + for(unsigned int j = 0; j < long_lag; ++j) + x[j] = gen01(); + i = long_lag; + } + + template + void seed(It& first, It last) + { +#ifndef BOOST_NO_STDC_NAMESPACE + // allow for Koenig lookup + using std::fmod; + using std::pow; +#endif + unsigned long mask = ~((~0u) << (w%32)); // now lowest w bits set + RealType two32 = pow(RealType(2), 32); + unsigned int j; + for(j = 0; j < long_lag && first != last; ++j) { + x[j] = RealType(0); + for(int k = 0; k < w/32 && first != last; ++k, ++first) + x[j] += *first / pow(two32,k+1); + if(first != last && mask != 0) { + x[j] += fmod((*first & mask) / _modulus, RealType(1)); + ++first; + } + } + i = long_lag; + if(first == last && j < long_lag) + throw std::invalid_argument("lagged_fibonacci_01::seed"); + } + + result_type min BOOST_PREVENT_MACRO_SUBSTITUTION () const { return result_type(0); } + result_type max BOOST_PREVENT_MACRO_SUBSTITUTION () const { return result_type(1); } + + result_type operator()() + { + if(i >= long_lag) + fill(); + return x[i++]; + } + + static bool validation(result_type x) + { + result_type v = fibonacci_validation::value(); + result_type epsilon = fibonacci_validation::tolerance(); + // std::abs is a source of trouble: sometimes, it's not overloaded + // for double, plus the usual namespace std noncompliance -> avoid it + // using std::abs; + // return abs(x - v) < 5 * epsilon + return x > v - epsilon && x < v + epsilon; + } + +#ifndef BOOST_NO_OPERATORS_IN_NAMESPACE + +#ifndef BOOST_RANDOM_NO_STREAM_OPERATORS + template + friend std::basic_ostream& + operator<<(std::basic_ostream& os, const lagged_fibonacci_01&f) + { +#ifndef BOOST_NO_STDC_NAMESPACE + // allow for Koenig lookup + using std::pow; +#endif + os << f.i << " "; + std::ios_base::fmtflags oldflags = os.flags(os.dec | os.fixed | os.left); + for(unsigned int i = 0; i < f.long_lag; ++i) + os << f.x[i] * f._modulus << " "; + os.flags(oldflags); + return os; + } + + template + friend std::basic_istream& + operator>>(std::basic_istream& is, lagged_fibonacci_01& f) + { +# ifdef BOOST_RANDOM_EXTRACT_LF + return detail::extract_lagged_fibonacci_01(is, f, f.i, f.x, f._modulus); +# else + is >> f.i >> std::ws; + for(unsigned int i = 0; i < f.long_lag; ++i) { + typename lagged_fibonacci_01::result_type value; + is >> value >> std::ws; + f.x[i] = value / f._modulus; + } + return is; +# endif + } +#endif + + friend bool operator==(const lagged_fibonacci_01& x, + const lagged_fibonacci_01& y) + { return x.i == y.i && std::equal(x.x, x.x+long_lag, y.x); } + friend bool operator!=(const lagged_fibonacci_01& x, + const lagged_fibonacci_01& y) + { return !(x == y); } +#else + // Use a member function; Streamable concept not supported. + bool operator==(const lagged_fibonacci_01& rhs) const + { return i == rhs.i && std::equal(x, x+long_lag, rhs.x); } + bool operator!=(const lagged_fibonacci_01& rhs) const + { return !(*this == rhs); } +#endif + +private: + /// \cond hide_private_members + void fill(); + /// \endcond + unsigned int i; + RealType x[long_lag]; + RealType _modulus; +}; + +#ifndef BOOST_NO_INCLASS_MEMBER_INITIALIZATION +// A definition is required even for integral static constants +template +const bool lagged_fibonacci_01::has_fixed_range; +template +const unsigned int lagged_fibonacci_01::long_lag; +template +const unsigned int lagged_fibonacci_01::short_lag; +template +const int lagged_fibonacci_01::word_size; + +#endif + +/// \cond hide_private_members +template +void lagged_fibonacci_01::fill() +{ + // two loops to avoid costly modulo operations + { // extra scope for MSVC brokenness w.r.t. for scope + for(unsigned int j = 0; j < short_lag; ++j) { + RealType t = x[j] + x[j+(long_lag-short_lag)]; + if(t >= RealType(1)) + t -= RealType(1); + x[j] = t; + } + } + for(unsigned int j = short_lag; j < long_lag; ++j) { + RealType t = x[j] + x[j-short_lag]; + if(t >= RealType(1)) + t -= RealType(1); + x[j] = t; + } + i = 0; +} +/// \endcond + +} // namespace random + +#ifdef BOOST_RANDOM_DOXYGEN +namespace detail { +/** + * The specializations lagged_fibonacci607 ... lagged_fibonacci44497 + * use well tested lags. + * + * See + * + * @blockquote + * "On the Periods of Generalized Fibonacci Recurrences", Richard P. Brent + * Computer Sciences Laboratory Australian National University, December 1992 + * @endblockquote + * + * The lags used here can be found in + * + * @blockquote + * "Uniform random number generators for supercomputers", Richard Brent, + * Proc. of Fifth Australian Supercomputer Conference, Melbourne, + * Dec. 1992, pp. 704-706. + * @endblockquote + */ +struct lagged_fibonacci_doc {}; +} +#endif + +/** + * @copydoc boost::detail::lagged_fibonacci_doc + */ +typedef random::lagged_fibonacci_01 lagged_fibonacci607; +/** + * @copydoc boost::detail::lagged_fibonacci_doc + */ +typedef random::lagged_fibonacci_01 lagged_fibonacci1279; +/** + * @copydoc boost::detail::lagged_fibonacci_doc + */ +typedef random::lagged_fibonacci_01 lagged_fibonacci2281; +/** + * @copydoc boost::detail::lagged_fibonacci_doc + */ +typedef random::lagged_fibonacci_01 lagged_fibonacci3217; +/** + * @copydoc boost::detail::lagged_fibonacci_doc + */ +typedef random::lagged_fibonacci_01 lagged_fibonacci4423; +/** + * @copydoc boost::detail::lagged_fibonacci_doc + */ +typedef random::lagged_fibonacci_01 lagged_fibonacci9689; +/** + * @copydoc boost::detail::lagged_fibonacci_doc + */ +typedef random::lagged_fibonacci_01 lagged_fibonacci19937; +/** + * @copydoc boost::detail::lagged_fibonacci_doc + */ +typedef random::lagged_fibonacci_01 lagged_fibonacci23209; +/** + * @copydoc boost::detail::lagged_fibonacci_doc + */ +typedef random::lagged_fibonacci_01 lagged_fibonacci44497; + + +// It is possible to partially specialize uniform_01<> on lagged_fibonacci_01<> +// to help the compiler generate efficient code. For GCC, this seems useless, +// because GCC optimizes (x-0)/(1-0) to (x-0). This is good enough for now. + +} // namespace boost + +#endif // BOOST_RANDOM_LAGGED_FIBONACCI_HPP diff --git a/Boost/boost/random/linear_congruential.hpp b/Boost/boost/random/linear_congruential.hpp new file mode 100644 index 000000000..d7277f156 --- /dev/null +++ b/Boost/boost/random/linear_congruential.hpp @@ -0,0 +1,403 @@ +/* boost random/linear_congruential.hpp header file + * + * Copyright Jens Maurer 2000-2001 + * 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: linear_congruential.hpp 60755 2010-03-22 00:45:06Z steven_watanabe $ + * + * Revision history + * 2001-02-18 moved to individual header files + */ + +#ifndef BOOST_RANDOM_LINEAR_CONGRUENTIAL_HPP +#define BOOST_RANDOM_LINEAR_CONGRUENTIAL_HPP + +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include + +namespace boost { +namespace random { + +/** + * Instantiations of class template linear_congruential model a + * \pseudo_random_number_generator. Linear congruential pseudo-random + * number generators are described in: + * + * "Mathematical methods in large-scale computing units", D. H. Lehmer, + * Proc. 2nd Symposium on Large-Scale Digital Calculating Machines, + * Harvard University Press, 1951, pp. 141-146 + * + * Let x(n) denote the sequence of numbers returned by some pseudo-random + * number generator. Then for the linear congruential generator, + * x(n+1) := (a * x(n) + c) mod m. Parameters for the generator are + * x(0), a, c, m. The template parameter IntType shall denote an integral + * type. It must be large enough to hold values a, c, and m. The template + * parameters a and c must be smaller than m. + * + * Note: The quality of the generator crucially depends on the choice of + * the parameters. User code should use one of the sensibly parameterized + * generators such as minstd_rand instead. + */ +template +class linear_congruential +{ +public: + typedef IntType result_type; +#ifndef BOOST_NO_INCLASS_MEMBER_INITIALIZATION + static const bool has_fixed_range = true; + static const result_type min_value = ( c == 0 ? 1 : 0 ); + static const result_type max_value = m-1; +#else + BOOST_STATIC_CONSTANT(bool, has_fixed_range = false); +#endif + BOOST_STATIC_CONSTANT(IntType, multiplier = a); + BOOST_STATIC_CONSTANT(IntType, increment = c); + BOOST_STATIC_CONSTANT(IntType, modulus = m); + + // MSVC 6 and possibly others crash when encountering complicated integral + // constant expressions. Avoid the check for now. + // BOOST_STATIC_ASSERT(m == 0 || a < m); + // BOOST_STATIC_ASSERT(m == 0 || c < m); + + /** + * Constructs a linear_congruential generator, seeding it with @c x0. + */ + explicit linear_congruential(IntType x0 = 1) + { + seed(x0); + + // MSVC fails BOOST_STATIC_ASSERT with std::numeric_limits at class scope +#ifndef BOOST_NO_LIMITS_COMPILE_TIME_CONSTANTS + BOOST_STATIC_ASSERT(std::numeric_limits::is_integer); +#endif + } + + /** + * Constructs a @c linear_congruential generator and seeds it + * with values taken from the itrator range [first, last) + * and adjusts first to point to the element after the last one + * used. If there are not enough elements, throws @c std::invalid_argument. + * + * first and last must be input iterators. + */ + template + linear_congruential(It& first, It last) + { + seed(first, last); + } + + // compiler-generated copy constructor and assignment operator are fine + + /** + * If c mod m is zero and x0 mod m is zero, changes the current value of + * the generator to 1. Otherwise, changes it to x0 mod m. If c is zero, + * distinct seeds in the range [1,m) will leave the generator in distinct + * states. If c is not zero, the range is [0,m). + */ + void seed(IntType x0 = 1) + { + // wrap _x if it doesn't fit in the destination + if(modulus == 0) { + _x = x0; + } else { + _x = x0 % modulus; + } + // handle negative seeds + if(_x <= 0 && _x != 0) { + _x += modulus; + } + // adjust to the correct range + if(increment == 0 && _x == 0) { + _x = 1; + } + assert(_x >= (min)()); + assert(_x <= (max)()); + } + + /** + * seeds a @c linear_congruential generator with values taken + * from the itrator range [first, last) and adjusts @c first to + * point to the element after the last one used. If there are + * not enough elements, throws @c std::invalid_argument. + * + * @c first and @c last must be input iterators. + */ + template + void seed(It& first, It last) + { + if(first == last) + throw std::invalid_argument("linear_congruential::seed"); + seed(*first++); + } + + /** + * Returns the smallest value that the @c linear_congruential generator + * can produce. + */ + result_type min BOOST_PREVENT_MACRO_SUBSTITUTION () const { return c == 0 ? 1 : 0; } + /** + * Returns the largest value that the @c linear_congruential generator + * can produce. + */ + result_type max BOOST_PREVENT_MACRO_SUBSTITUTION () const { return modulus-1; } + + /** Returns the next value of the @c linear_congruential generator. */ + IntType operator()() + { + _x = const_mod::mult_add(a, _x, c); + return _x; + } + + static bool validation(IntType x) { return val == x; } + +#ifdef BOOST_NO_OPERATORS_IN_NAMESPACE + + // Use a member function; Streamable concept not supported. + bool operator==(const linear_congruential& rhs) const + { return _x == rhs._x; } + bool operator!=(const linear_congruential& rhs) const + { return !(*this == rhs); } + +#else + friend bool operator==(const linear_congruential& x, + const linear_congruential& y) + { return x._x == y._x; } + friend bool operator!=(const linear_congruential& x, + const linear_congruential& y) + { return !(x == y); } + +#if !defined(BOOST_RANDOM_NO_STREAM_OPERATORS) && !BOOST_WORKAROUND(__BORLANDC__, BOOST_TESTED_AT(0x551)) + template + friend std::basic_ostream& + operator<<(std::basic_ostream& os, + const linear_congruential& lcg) + { + return os << lcg._x; + } + + template + friend std::basic_istream& + operator>>(std::basic_istream& is, + linear_congruential& lcg) + { + return is >> lcg._x; + } + +private: +#endif +#endif + + IntType _x; +}; + +// probably needs the "no native streams" caveat for STLPort +#if !defined(__SGI_STL_PORT) && BOOST_WORKAROUND(__GNUC__, == 2) +template +std::ostream& +operator<<(std::ostream& os, + const linear_congruential& lcg) +{ + return os << lcg._x; +} + +template +std::istream& +operator>>(std::istream& is, + linear_congruential& lcg) +{ + return is >> lcg._x; +} +#elif defined(BOOST_RANDOM_NO_STREAM_OPERATORS) || BOOST_WORKAROUND(__BORLANDC__, BOOST_TESTED_AT(0x551)) +template +std::basic_ostream& +operator<<(std::basic_ostream& os, + const linear_congruential& lcg) +{ + return os << lcg._x; +} + +template +std::basic_istream& +operator>>(std::basic_istream& is, + linear_congruential& lcg) +{ + return is >> lcg._x; +} +#endif + +#ifndef BOOST_NO_INCLASS_MEMBER_INITIALIZATION +// A definition is required even for integral static constants +template +const bool linear_congruential::has_fixed_range; +template +const typename linear_congruential::result_type linear_congruential::min_value; +template +const typename linear_congruential::result_type linear_congruential::max_value; +template +const IntType linear_congruential::modulus; +#endif + +} // namespace random + +// validation values from the publications +/** + * The specialization \minstd_rand0 was originally suggested in + * + * @blockquote + * A pseudo-random number generator for the System/360, P.A. Lewis, + * A.S. Goodman, J.M. Miller, IBM Systems Journal, Vol. 8, No. 2, + * 1969, pp. 136-146 + * @endblockquote + * + * It is examined more closely together with \minstd_rand in + * + * @blockquote + * "Random Number Generators: Good ones are hard to find", + * Stephen K. Park and Keith W. Miller, Communications of + * the ACM, Vol. 31, No. 10, October 1988, pp. 1192-1201 + * @endblockquote + */ +typedef random::linear_congruential minstd_rand0; + +/** The specialization \minstd_rand was suggested in + * + * @blockquote + * "Random Number Generators: Good ones are hard to find", + * Stephen K. Park and Keith W. Miller, Communications of + * the ACM, Vol. 31, No. 10, October 1988, pp. 1192-1201 + * @endblockquote + */ +typedef random::linear_congruential minstd_rand; + + +#if !defined(BOOST_NO_INT64_T) && !defined(BOOST_NO_INTEGRAL_INT64_T) +/** Class @c rand48 models a \pseudo_random_number_generator. It uses + * the linear congruential algorithm with the parameters a = 0x5DEECE66D, + * c = 0xB, m = 2**48. It delivers identical results to the @c lrand48() + * function available on some systems (assuming lcong48 has not been called). + * + * It is only available on systems where @c uint64_t is provided as an + * integral type, so that for example static in-class constants and/or + * enum definitions with large @c uint64_t numbers work. + */ +class rand48 +{ +public: + typedef int32_t result_type; +#ifndef BOOST_NO_INCLASS_MEMBER_INITIALIZATION + static const bool has_fixed_range = true; + static const int32_t min_value = 0; + static const int32_t max_value = integer_traits::const_max; +#else + enum { has_fixed_range = false }; +#endif + /** + * Returns the smallest value that the generator can produce + */ + int32_t min BOOST_PREVENT_MACRO_SUBSTITUTION () const { return 0; } + /** + * Returns the largest value that the generator can produce + */ + int32_t max BOOST_PREVENT_MACRO_SUBSTITUTION () const { return std::numeric_limits::max BOOST_PREVENT_MACRO_SUBSTITUTION (); } + +#ifdef BOOST_RANDOM_DOXYGEN + /** + * If T is an integral type smaller than int46_t, constructs + * a \rand48 generator with x(0) := (x0 << 16) | 0x330e. Otherwise + * constructs a \rand48 generator with x(0) = x0. + */ + template explicit rand48(T x0 = 1); +#else + rand48() : lcf(cnv(static_cast(1))) {} + template explicit rand48(T x0) : lcf(cnv(x0)) { } +#endif + template rand48(It& first, It last) : lcf(first, last) { } + + // compiler-generated copy ctor and assignment operator are fine + +#ifdef BOOST_RANDOM_DOXYGEN + /** + * If T is an integral type smaller than int46_t, changes + * the current value x(n) of the generator to (x0 << 16) | 0x330e. + * Otherwise changes the current value x(n) to x0. + */ + template void seed(T x0 = 1); +#else + void seed() { seed(static_cast(1)); } + template void seed(T x0) { lcf.seed(cnv(x0)); } +#endif + template void seed(It& first, It last) { lcf.seed(first,last); } + + /** + * Returns the next value of the generator. + */ + int32_t operator()() { return static_cast(lcf() >> 17); } + // by experiment from lrand48() + static bool validation(int32_t x) { return x == 1993516219; } + +#ifndef BOOST_NO_OPERATORS_IN_NAMESPACE + +#ifndef BOOST_RANDOM_NO_STREAM_OPERATORS + template + friend std::basic_ostream& + operator<<(std::basic_ostream& os, const rand48& r) + { os << r.lcf; return os; } + + template + friend std::basic_istream& + operator>>(std::basic_istream& is, rand48& r) + { is >> r.lcf; return is; } +#endif + + friend bool operator==(const rand48& x, const rand48& y) + { return x.lcf == y.lcf; } + friend bool operator!=(const rand48& x, const rand48& y) + { return !(x == y); } +#else + // Use a member function; Streamable concept not supported. + bool operator==(const rand48& rhs) const + { return lcf == rhs.lcf; } + bool operator!=(const rand48& rhs) const + { return !(*this == rhs); } +#endif +private: + /// \cond hide_private_members + random::linear_congruential lcf; + template + static uint64_t cnv(T x) + { + if(sizeof(T) < sizeof(uint64_t)) { + return (static_cast(x) << 16) | 0x330e; + } else { + return(static_cast(x)); + } + } + static uint64_t cnv(float x) { return(static_cast(x)); } + static uint64_t cnv(double x) { return(static_cast(x)); } + static uint64_t cnv(long double x) { return(static_cast(x)); } + /// \endcond +}; +#endif /* !BOOST_NO_INT64_T && !BOOST_NO_INTEGRAL_INT64_T */ + +} // namespace boost + +#include + +#endif // BOOST_RANDOM_LINEAR_CONGRUENTIAL_HPP diff --git a/Boost/boost/random/linear_feedback_shift.hpp b/Boost/boost/random/linear_feedback_shift.hpp new file mode 100644 index 000000000..6b3a87826 --- /dev/null +++ b/Boost/boost/random/linear_feedback_shift.hpp @@ -0,0 +1,160 @@ +/* boost random/tausworthe.hpp header file + * + * Copyright Jens Maurer 2002 + * 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: linear_feedback_shift.hpp 60755 2010-03-22 00:45:06Z steven_watanabe $ + * + */ + +#ifndef BOOST_RANDOM_LINEAR_FEEDBACK_SHIFT_HPP +#define BOOST_RANDOM_LINEAR_FEEDBACK_SHIFT_HPP + +#include +#include +#include +#include +#include +#include +#include + +namespace boost { +namespace random { + +/** + * Instatiation of @c linear_feedback_shift model a + * \pseudo_random_number_generator. It was originally + * proposed in + * + * @blockquote + * "Random numbers generated by linear recurrence modulo two.", + * Tausworthe, R. C.(1965), Mathematics of Computation 19, 201-209. + * @endblockquote + */ +template +class linear_feedback_shift +{ +public: + typedef UIntType result_type; + // avoid the warning trouble when using (1< 0); + // BOOST_STATIC_ASSERT(q > 0); + // BOOST_STATIC_ASSERT(k < w); + // BOOST_STATIC_ASSERT(0 < 2*q && 2*q < k); + // BOOST_STATIC_ASSERT(0 < s && s <= k-q); + + explicit linear_feedback_shift(UIntType s0 = 341) : wordmask(0) + { + // MSVC fails BOOST_STATIC_ASSERT with std::numeric_limits at class scope +#ifndef BOOST_NO_LIMITS_COMPILE_TIME_CONSTANTS + BOOST_STATIC_ASSERT(std::numeric_limits::is_integer); + BOOST_STATIC_ASSERT(!std::numeric_limits::is_signed); +#endif + + // avoid "left shift count >= with of type" warning + for(int i = 0; i < w; ++i) + wordmask |= (1u << i); + seed(s0); + } + + template linear_feedback_shift(It& first, It last) : wordmask(0) + { + // MSVC fails BOOST_STATIC_ASSERT with std::numeric_limits at class scope +#ifndef BOOST_NO_LIMITS_COMPILE_TIME_CONSTANTS + BOOST_STATIC_ASSERT(std::numeric_limits::is_integer); + BOOST_STATIC_ASSERT(!std::numeric_limits::is_signed); +#endif + + // avoid "left shift count >= with of type" warning + for(int i = 0; i < w; ++i) + wordmask |= (1u << i); + seed(first, last); + } + + void seed(UIntType s0 = 341) { + if(s0 < (1 << (w-k))) { + s0 += 1 << (w-k); + } + value = s0; + } + template void seed(It& first, It last) + { + if(first == last) + throw std::invalid_argument("linear_feedback_shift::seed"); + value = *first++; + assert(value >= (1 << (w-k))); + } + + result_type operator()() + { + const UIntType b = (((value << q) ^ value) & wordmask) >> (k-s); + const UIntType mask = ( (~static_cast(0)) << (w-k) ) & wordmask; + value = ((value & mask) << s) ^ b; + return value; + } + static bool validation(result_type x) { return val == x; } + +#ifndef BOOST_NO_OPERATORS_IN_NAMESPACE + +#ifndef BOOST_RANDOM_NO_STREAM_OPERATORS + template + friend std::basic_ostream& + operator<<(std::basic_ostream& os, linear_feedback_shift x) + { os << x.value; return os; } + + template + friend std::basic_istream& + operator>>(std::basic_istream& is, linear_feedback_shift& x) + { is >> x.value; return is; } +#endif + + friend bool operator==(linear_feedback_shift x, linear_feedback_shift y) + { return x.value == y.value; } + friend bool operator!=(linear_feedback_shift x, linear_feedback_shift y) + { return !(x == y); } +#else + // Use a member function; Streamable concept not supported. + bool operator==(linear_feedback_shift rhs) const + { return value == rhs.value; } + bool operator!=(linear_feedback_shift rhs) const + { return !(*this == rhs); } +#endif + +private: + UIntType wordmask; // avoid "left shift count >= width of type" warnings + UIntType value; +}; + +#ifndef BOOST_NO_INCLASS_MEMBER_INITIALIZATION +// A definition is required even for integral static constants +template +const bool linear_feedback_shift::has_fixed_range; +template +const int linear_feedback_shift::word_size; +template +const int linear_feedback_shift::exponent1; +template +const int linear_feedback_shift::exponent2; +template +const int linear_feedback_shift::step_size; +#endif + +} // namespace random +} // namespace boost + +#endif // BOOST_RANDOM_LINEAR_FEEDBACK_SHIFT_HPP diff --git a/Boost/boost/random/lognormal_distribution.hpp b/Boost/boost/random/lognormal_distribution.hpp new file mode 100644 index 000000000..bb54829e2 --- /dev/null +++ b/Boost/boost/random/lognormal_distribution.hpp @@ -0,0 +1,129 @@ +/* boost random/lognormal_distribution.hpp header file + * + * Copyright Jens Maurer 2000-2001 + * 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: lognormal_distribution.hpp 60755 2010-03-22 00:45:06Z steven_watanabe $ + * + * Revision history + * 2001-02-18 moved to individual header files + */ + +#ifndef BOOST_RANDOM_LOGNORMAL_DISTRIBUTION_HPP +#define BOOST_RANDOM_LOGNORMAL_DISTRIBUTION_HPP + +#include // std::exp, std::sqrt +#include +#include +#include +#include +#include +#include + +#ifdef BOOST_NO_STDC_NAMESPACE +namespace std { + using ::log; + using ::sqrt; +} +#endif + +namespace boost { + +#if defined(__GNUC__) && (__GNUC__ < 3) +// Special gcc workaround: gcc 2.95.x ignores using-declarations +// in template classes (confirmed by gcc author Martin v. Loewis) + using std::sqrt; + using std::exp; +#endif + +/** + * Instantiations of class template lognormal_distribution model a + * \random_distribution. Such a distribution produces random numbers + * with \f$p(x) = \frac{1}{x \sigma_N \sqrt{2\pi}} e^{\frac{-\left(\log(x)-\mu_N\right)^2}{2\sigma_N^2}}\f$ + * for x > 0, where \f$\mu_N = \log\left(\frac{\mu^2}{\sqrt{\sigma^2 + \mu^2}}\right)\f$ and + * \f$\sigma_N = \sqrt{\log\left(1 + \frac{\sigma^2}{\mu^2}\right)}\f$. + */ +template +class lognormal_distribution +{ +public: + typedef typename normal_distribution::input_type input_type; + typedef RealType result_type; + +#ifndef BOOST_NO_LIMITS_COMPILE_TIME_CONSTANTS + BOOST_STATIC_ASSERT(!std::numeric_limits::is_integer); +#endif + + /** + * Constructs a lognormal_distribution. @c mean and @c sigma are the + * mean and standard deviation of the lognormal distribution. + */ + explicit lognormal_distribution(result_type mean_arg = result_type(1), + result_type sigma_arg = result_type(1)) + : _mean(mean_arg), _sigma(sigma_arg) + { + assert(_mean > result_type(0)); + init(); + } + + // compiler-generated copy ctor and assignment operator are fine + + RealType mean() const { return _mean; } + RealType sigma() const { return _sigma; } + void reset() { _normal.reset(); } + + template + result_type operator()(Engine& eng) + { +#ifndef BOOST_NO_STDC_NAMESPACE + // allow for Koenig lookup + using std::exp; +#endif + return exp(_normal(eng) * _nsigma + _nmean); + } + +#ifndef BOOST_RANDOM_NO_STREAM_OPERATORS + template + friend std::basic_ostream& + operator<<(std::basic_ostream& os, const lognormal_distribution& ld) + { + os << ld._normal << " " << ld._mean << " " << ld._sigma; + return os; + } + + template + friend std::basic_istream& + operator>>(std::basic_istream& is, lognormal_distribution& ld) + { + is >> std::ws >> ld._normal >> std::ws >> ld._mean >> std::ws >> ld._sigma; + ld.init(); + return is; + } +#endif + +private: + + /// \cond hide_private_members + void init() + { +#ifndef BOOST_NO_STDC_NAMESPACE + // allow for Koenig lookup + using std::exp; using std::log; using std::sqrt; +#endif + _nmean = log(_mean*_mean/sqrt(_sigma*_sigma + _mean*_mean)); + _nsigma = sqrt(log(_sigma*_sigma/_mean/_mean+result_type(1))); + } + /// \endcond + + RealType _mean, _sigma; + RealType _nmean, _nsigma; + normal_distribution _normal; +}; + +} // namespace boost + +#endif // BOOST_RANDOM_LOGNORMAL_DISTRIBUTION_HPP diff --git a/Boost/boost/random/mersenne_twister.hpp b/Boost/boost/random/mersenne_twister.hpp new file mode 100644 index 000000000..fae5cd514 --- /dev/null +++ b/Boost/boost/random/mersenne_twister.hpp @@ -0,0 +1,367 @@ +/* boost random/mersenne_twister.hpp header file + * + * Copyright Jens Maurer 2000-2001 + * 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: mersenne_twister.hpp 60755 2010-03-22 00:45:06Z steven_watanabe $ + * + * Revision history + * 2001-02-18 moved to individual header files + */ + +#ifndef BOOST_RANDOM_MERSENNE_TWISTER_HPP +#define BOOST_RANDOM_MERSENNE_TWISTER_HPP + +#include +#include // std::copy +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +namespace boost { +namespace random { + +/** + * Instantiations of class template mersenne_twister model a + * \pseudo_random_number_generator. It uses the algorithm described in + * + * @blockquote + * "Mersenne Twister: A 623-dimensionally equidistributed uniform + * pseudo-random number generator", Makoto Matsumoto and Takuji Nishimura, + * ACM Transactions on Modeling and Computer Simulation: Special Issue on + * Uniform Random Number Generation, Vol. 8, No. 1, January 1998, pp. 3-30. + * @endblockquote + * + * @xmlnote + * The boost variant has been implemented from scratch and does not + * derive from or use mt19937.c provided on the above WWW site. However, it + * was verified that both produce identical output. + * @endxmlnote + * + * The seeding from an integer was changed in April 2005 to address a + * weakness. + * + * The quality of the generator crucially depends on the choice of the + * parameters. User code should employ one of the sensibly parameterized + * generators such as \mt19937 instead. + * + * The generator requires considerable amounts of memory for the storage of + * its state array. For example, \mt11213b requires about 1408 bytes and + * \mt19937 requires about 2496 bytes. + */ +template +class mersenne_twister +{ +public: + typedef UIntType result_type; + BOOST_STATIC_CONSTANT(int, word_size = w); + BOOST_STATIC_CONSTANT(int, state_size = n); + BOOST_STATIC_CONSTANT(int, shift_size = m); + BOOST_STATIC_CONSTANT(int, mask_bits = r); + BOOST_STATIC_CONSTANT(UIntType, parameter_a = a); + BOOST_STATIC_CONSTANT(int, output_u = u); + BOOST_STATIC_CONSTANT(int, output_s = s); + BOOST_STATIC_CONSTANT(UIntType, output_b = b); + BOOST_STATIC_CONSTANT(int, output_t = t); + BOOST_STATIC_CONSTANT(UIntType, output_c = c); + BOOST_STATIC_CONSTANT(int, output_l = l); + + BOOST_STATIC_CONSTANT(bool, has_fixed_range = false); + + /** + * Constructs a @c mersenne_twister and calls @c seed(). + */ + mersenne_twister() { seed(); } + + /** + * Constructs a @c mersenne_twister and calls @c seed(value). + */ + BOOST_RANDOM_DETAIL_ARITHMETIC_CONSTRUCTOR(mersenne_twister, UIntType, value) + { seed(value); } + template mersenne_twister(It& first, It last) { seed(first,last); } + + /** + * Constructs a mersenne_twister and calls @c seed(gen). + * + * @xmlnote + * The copy constructor will always be preferred over + * the templated constructor. + * @endxmlnote + */ + BOOST_RANDOM_DETAIL_GENERATOR_CONSTRUCTOR(mersenne_twister, Generator, gen) + { seed(gen); } + + // compiler-generated copy ctor and assignment operator are fine + + /** Calls @c seed(result_type(5489)). */ + void seed() { seed(UIntType(5489)); } + + /** + * Sets the state x(0) to v mod 2w. Then, iteratively, + * sets x(i) to (i + 1812433253 * (x(i-1) xor (x(i-1) rshift w-2))) mod 2w + * for i = 1 .. n-1. x(n) is the first value to be returned by operator(). + */ + BOOST_RANDOM_DETAIL_ARITHMETIC_SEED(mersenne_twister, UIntType, value) + { + // New seeding algorithm from + // http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/MT2002/emt19937ar.html + // In the previous versions, MSBs of the seed affected only MSBs of the + // state x[]. + const UIntType mask = ~0u; + x[0] = value & mask; + for (i = 1; i < n; i++) { + // See Knuth "The Art of Computer Programming" Vol. 2, 3rd ed., page 106 + x[i] = (1812433253UL * (x[i-1] ^ (x[i-1] >> (w-2))) + i) & mask; + } + } + + /** + * Sets the state of this mersenne_twister to the values + * returned by n invocations of gen. + * + * Complexity: Exactly n invocations of gen. + */ + BOOST_RANDOM_DETAIL_GENERATOR_SEED(mersenne_twister, Generator, gen) + { +#ifndef BOOST_NO_LIMITS_COMPILE_TIME_CONSTANTS + BOOST_STATIC_ASSERT(!std::numeric_limits::is_signed); +#endif + // I could have used std::generate_n, but it takes "gen" by value + for(int j = 0; j < n; j++) + x[j] = gen(); + i = n; + } + + template + void seed(It& first, It last) + { + int j; + for(j = 0; j < n && first != last; ++j, ++first) + x[j] = *first; + i = n; + if(first == last && j < n) + throw std::invalid_argument("mersenne_twister::seed"); + } + + result_type min BOOST_PREVENT_MACRO_SUBSTITUTION () const { return 0; } + result_type max BOOST_PREVENT_MACRO_SUBSTITUTION () const + { + // avoid "left shift count >= with of type" warning + result_type res = 0; + for(int j = 0; j < w; ++j) + res |= (1u << j); + return res; + } + + result_type operator()(); + static bool validation(result_type v) { return val == v; } + +#ifndef BOOST_NO_OPERATORS_IN_NAMESPACE + +#ifndef BOOST_RANDOM_NO_STREAM_OPERATORS + template + friend std::basic_ostream& + operator<<(std::basic_ostream& os, const mersenne_twister& mt) + { + for(int j = 0; j < mt.state_size; ++j) + os << mt.compute(j) << " "; + return os; + } + + template + friend std::basic_istream& + operator>>(std::basic_istream& is, mersenne_twister& mt) + { + for(int j = 0; j < mt.state_size; ++j) + is >> mt.x[j] >> std::ws; + // MSVC (up to 7.1) and Borland (up to 5.64) don't handle the template + // value parameter "n" available from the class template scope, so use + // the static constant with the same value + mt.i = mt.state_size; + return is; + } +#endif + + friend bool operator==(const mersenne_twister& x, const mersenne_twister& y) + { + for(int j = 0; j < state_size; ++j) + if(x.compute(j) != y.compute(j)) + return false; + return true; + } + + friend bool operator!=(const mersenne_twister& x, const mersenne_twister& y) + { return !(x == y); } +#else + // Use a member function; Streamable concept not supported. + bool operator==(const mersenne_twister& rhs) const + { + for(int j = 0; j < state_size; ++j) + if(compute(j) != rhs.compute(j)) + return false; + return true; + } + + bool operator!=(const mersenne_twister& rhs) const + { return !(*this == rhs); } +#endif + +private: + /// \cond hide_private_members + // returns x(i-n+index), where index is in 0..n-1 + UIntType compute(unsigned int index) const + { + // equivalent to (i-n+index) % 2n, but doesn't produce negative numbers + return x[ (i + n + index) % (2*n) ]; + } + void twist(int block); + /// \endcond + + // state representation: next output is o(x(i)) + // x[0] ... x[k] x[k+1] ... x[n-1] x[n] ... x[2*n-1] represents + // x(i-k) ... x(i) x(i+1) ... x(i-k+n-1) x(i-k-n) ... x[i(i-k-1)] + // The goal is to always have x(i-n) ... x(i-1) available for + // operator== and save/restore. + + UIntType x[2*n]; + int i; +}; + +#ifndef BOOST_NO_INCLASS_MEMBER_INITIALIZATION +// A definition is required even for integral static constants +template +const bool mersenne_twister::has_fixed_range; +template +const int mersenne_twister::state_size; +template +const int mersenne_twister::shift_size; +template +const int mersenne_twister::mask_bits; +template +const UIntType mersenne_twister::parameter_a; +template +const int mersenne_twister::output_u; +template +const int mersenne_twister::output_s; +template +const UIntType mersenne_twister::output_b; +template +const int mersenne_twister::output_t; +template +const UIntType mersenne_twister::output_c; +template +const int mersenne_twister::output_l; +#endif + +/// \cond hide_private_members +template +void mersenne_twister::twist(int block) +{ + const UIntType upper_mask = (~0u) << r; + const UIntType lower_mask = ~upper_mask; + + if(block == 0) { + for(int j = n; j < 2*n; j++) { + UIntType y = (x[j-n] & upper_mask) | (x[j-(n-1)] & lower_mask); + x[j] = x[j-(n-m)] ^ (y >> 1) ^ (y&1 ? a : 0); + } + } else if (block == 1) { + // split loop to avoid costly modulo operations + { // extra scope for MSVC brokenness w.r.t. for scope + for(int j = 0; j < n-m; j++) { + UIntType y = (x[j+n] & upper_mask) | (x[j+n+1] & lower_mask); + x[j] = x[j+n+m] ^ (y >> 1) ^ (y&1 ? a : 0); + } + } + + for(int j = n-m; j < n-1; j++) { + UIntType y = (x[j+n] & upper_mask) | (x[j+n+1] & lower_mask); + x[j] = x[j-(n-m)] ^ (y >> 1) ^ (y&1 ? a : 0); + } + // last iteration + UIntType y = (x[2*n-1] & upper_mask) | (x[0] & lower_mask); + x[n-1] = x[m-1] ^ (y >> 1) ^ (y&1 ? a : 0); + i = 0; + } +} +/// \endcond + +template +inline typename mersenne_twister::result_type +mersenne_twister::operator()() +{ + if(i == n) + twist(0); + else if(i >= 2*n) + twist(1); + // Step 4 + UIntType z = x[i]; + ++i; + z ^= (z >> u); + z ^= ((z << s) & b); + z ^= ((z << t) & c); + z ^= (z >> l); + return z; +} + +} // namespace random + +/** + * The specializations \mt11213b and \mt19937 are from + * + * @blockquote + * "Mersenne Twister: A 623-dimensionally equidistributed + * uniform pseudo-random number generator", Makoto Matsumoto + * and Takuji Nishimura, ACM Transactions on Modeling and + * Computer Simulation: Special Issue on Uniform Random Number + * Generation, Vol. 8, No. 1, January 1998, pp. 3-30. + * @endblockquote + */ +typedef random::mersenne_twister mt11213b; + +/** + * The specializations \mt11213b and \mt19937 are from + * + * @blockquote + * "Mersenne Twister: A 623-dimensionally equidistributed + * uniform pseudo-random number generator", Makoto Matsumoto + * and Takuji Nishimura, ACM Transactions on Modeling and + * Computer Simulation: Special Issue on Uniform Random Number + * Generation, Vol. 8, No. 1, January 1998, pp. 3-30. + * @endblockquote + */ +typedef random::mersenne_twister mt19937; + +} // namespace boost + +BOOST_RANDOM_PTR_HELPER_SPEC(boost::mt19937) + +#endif // BOOST_RANDOM_MERSENNE_TWISTER_HPP diff --git a/Boost/boost/random/normal_distribution.hpp b/Boost/boost/random/normal_distribution.hpp new file mode 100644 index 000000000..f6317460c --- /dev/null +++ b/Boost/boost/random/normal_distribution.hpp @@ -0,0 +1,131 @@ +/* boost random/normal_distribution.hpp header file + * + * Copyright Jens Maurer 2000-2001 + * 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: normal_distribution.hpp 60755 2010-03-22 00:45:06Z steven_watanabe $ + * + * Revision history + * 2001-02-18 moved to individual header files + */ + +#ifndef BOOST_RANDOM_NORMAL_DISTRIBUTION_HPP +#define BOOST_RANDOM_NORMAL_DISTRIBUTION_HPP + +#include +#include +#include +#include +#include +#include + +namespace boost { + +/** + * Instantiations of class template normal_distribution model a + * \random_distribution. Such a distribution produces random numbers + * @c x distributed with probability density function + * \f$p(x) = \frac{1}{\sqrt{2\pi\sigma}} e^{-\frac{(x-\mu)^2}{2\sigma^2}}\f$, + * where mean and sigma are the parameters of the distribution. + */ +// deterministic Box-Muller method, uses trigonometric functions +template +class normal_distribution +{ +public: + typedef RealType input_type; + typedef RealType result_type; + +#if !defined(BOOST_NO_LIMITS_COMPILE_TIME_CONSTANTS) && !(defined(BOOST_MSVC) && BOOST_MSVC <= 1300) + BOOST_STATIC_ASSERT(!std::numeric_limits::is_integer); +#endif + + /** + * Constructs a normal_distribution object. @c mean and @c sigma are + * the parameters for the distribution. + * + * Requires: sigma > 0 + */ + explicit normal_distribution(const result_type& mean_arg = result_type(0), + const result_type& sigma_arg = result_type(1)) + : _mean(mean_arg), _sigma(sigma_arg), _valid(false) + { + assert(_sigma >= result_type(0)); + } + + // compiler-generated copy constructor is NOT fine, need to purge cache + normal_distribution(const normal_distribution& other) + : _mean(other._mean), _sigma(other._sigma), _valid(false) + { + } + + // compiler-generated copy ctor and assignment operator are fine + + /** + * Returns: The "mean" parameter of the distribution. + */ + RealType mean() const { return _mean; } + /** + * Returns: The "sigma" parameter of the distribution. + */ + RealType sigma() const { return _sigma; } + + void reset() { _valid = false; } + + template + result_type operator()(Engine& eng) + { +#ifndef BOOST_NO_STDC_NAMESPACE + // allow for Koenig lookup + using std::sqrt; using std::log; using std::sin; using std::cos; +#endif + if(!_valid) { + _r1 = eng(); + _r2 = eng(); + _cached_rho = sqrt(-result_type(2) * log(result_type(1)-_r2)); + _valid = true; + } else { + _valid = false; + } + // Can we have a boost::mathconst please? + const result_type pi = result_type(3.14159265358979323846); + + return _cached_rho * (_valid ? + cos(result_type(2)*pi*_r1) : + sin(result_type(2)*pi*_r1)) + * _sigma + _mean; + } + +#ifndef BOOST_RANDOM_NO_STREAM_OPERATORS + template + friend std::basic_ostream& + operator<<(std::basic_ostream& os, const normal_distribution& nd) + { + os << nd._mean << " " << nd._sigma << " " + << nd._valid << " " << nd._cached_rho << " " << nd._r1; + return os; + } + + template + friend std::basic_istream& + operator>>(std::basic_istream& is, normal_distribution& nd) + { + is >> std::ws >> nd._mean >> std::ws >> nd._sigma + >> std::ws >> nd._valid >> std::ws >> nd._cached_rho + >> std::ws >> nd._r1; + return is; + } +#endif +private: + result_type _mean, _sigma; + result_type _r1, _r2, _cached_rho; + bool _valid; +}; + +} // namespace boost + +#endif // BOOST_RANDOM_NORMAL_DISTRIBUTION_HPP diff --git a/Boost/boost/random/poisson_distribution.hpp b/Boost/boost/random/poisson_distribution.hpp new file mode 100644 index 000000000..ae72ee8a1 --- /dev/null +++ b/Boost/boost/random/poisson_distribution.hpp @@ -0,0 +1,116 @@ +/* boost random/poisson_distribution.hpp header file + * + * Copyright Jens Maurer 2002 + * 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: poisson_distribution.hpp 60755 2010-03-22 00:45:06Z steven_watanabe $ + * + */ + +#ifndef BOOST_RANDOM_POISSON_DISTRIBUTION_HPP +#define BOOST_RANDOM_POISSON_DISTRIBUTION_HPP + +#include +#include +#include +#include +#include +#include + +namespace boost { + +// Knuth + +/** + * An instantiation of the class template @c poisson_distribution is a + * model of \random_distribution. The poisson distribution has + * \f$p(i) = \frac{e^{-\lambda}\lambda^i}{i!}\f$ + */ +template +class poisson_distribution +{ +public: + typedef RealType input_type; + typedef IntType result_type; + + /** + * Constructs a @c poisson_distribution with the parameter @c mean. + * + * Requires: mean > 0 + */ + explicit poisson_distribution(const RealType& mean_arg = RealType(1)) + : _mean(mean_arg) + { +#ifndef BOOST_NO_LIMITS_COMPILE_TIME_CONSTANTS + // MSVC fails BOOST_STATIC_ASSERT with std::numeric_limits at class scope + BOOST_STATIC_ASSERT(std::numeric_limits::is_integer); + BOOST_STATIC_ASSERT(!std::numeric_limits::is_integer); +#endif + + assert(_mean > RealType(0)); + init(); + } + + // compiler-generated copy ctor and assignment operator are fine + + /** + * Returns: the "mean" parameter of the distribution. + */ + RealType mean() const { return _mean; } + void reset() { } + + template + result_type operator()(Engine& eng) + { + // TODO: This is O(_mean), but it should be O(log(_mean)) for large _mean + RealType product = RealType(1); + for(result_type m = 0; ; ++m) { + product *= eng(); + if(product <= _exp_mean) + return m; + } + } + +#ifndef BOOST_RANDOM_NO_STREAM_OPERATORS + template + friend std::basic_ostream& + operator<<(std::basic_ostream& os, const poisson_distribution& pd) + { + os << pd._mean; + return os; + } + + template + friend std::basic_istream& + operator>>(std::basic_istream& is, poisson_distribution& pd) + { + is >> std::ws >> pd._mean; + pd.init(); + return is; + } +#endif + +private: + /// \cond hide_private_members + void init() + { +#ifndef BOOST_NO_STDC_NAMESPACE + // allow for Koenig lookup + using std::exp; +#endif + _exp_mean = exp(-_mean); + } + /// \endcond + + RealType _mean; + // some precomputed data from the parameters + RealType _exp_mean; +}; + +} // namespace boost + +#endif // BOOST_RANDOM_POISSON_DISTRIBUTION_HPP diff --git a/Boost/boost/random/random_number_generator.hpp b/Boost/boost/random/random_number_generator.hpp new file mode 100644 index 000000000..7fb9733d2 --- /dev/null +++ b/Boost/boost/random/random_number_generator.hpp @@ -0,0 +1,71 @@ +/* boost random/random_number_generator.hpp header file + * + * Copyright Jens Maurer 2000-2001 + * 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: random_number_generator.hpp 60755 2010-03-22 00:45:06Z steven_watanabe $ + * + * Revision history + * 2001-02-18 moved to individual header files + */ + +#ifndef BOOST_RANDOM_RANDOM_NUMBER_GENERATOR_HPP +#define BOOST_RANDOM_RANDOM_NUMBER_GENERATOR_HPP + +#include +#include +#include +#include +#include + +namespace boost { + +/** + * Instantiations of class template random_number_generator model a + * RandomNumberGenerator (std:25.2.11 [lib.alg.random.shuffle]). On + * each invocation, it returns a uniformly distributed integer in + * the range [0..n). + * + * The template parameter IntType shall denote some integer-like value type. + */ +template +class random_number_generator +{ +public: + typedef UniformRandomNumberGenerator base_type; + typedef IntType argument_type; + typedef IntType result_type; + /** + * Constructs a random_number_generator functor with the given + * \uniform_random_number_generator as the underlying source of + * random numbers. + */ + random_number_generator(base_type& rng) : _rng(rng) + { +#ifndef BOOST_NO_LIMITS_COMPILE_TIME_CONSTANTS + BOOST_STATIC_ASSERT(std::numeric_limits::is_integer); +#endif + } + // compiler-generated copy ctor is fine + // assignment is disallowed because there is a reference member + + /** + * Returns a value in the range [0, n) + */ + result_type operator()(argument_type n) + { + typedef uniform_int dist_type; + return variate_generator(_rng, dist_type(0, n-1))(); + } + +private: + base_type& _rng; +}; + +} // namespace boost + +#endif // BOOST_RANDOM_RANDOM_NUMBER_GENERATOR_HPP diff --git a/Boost/boost/random/ranlux.hpp b/Boost/boost/random/ranlux.hpp new file mode 100644 index 000000000..39eb32b3d --- /dev/null +++ b/Boost/boost/random/ranlux.hpp @@ -0,0 +1,81 @@ +/* boost random/ranlux.hpp header file + * + * Copyright Jens Maurer 2002 + * 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: ranlux.hpp 60755 2010-03-22 00:45:06Z steven_watanabe $ + * + * Revision history + * 2001-02-18 created + */ + +#ifndef BOOST_RANDOM_RANLUX_HPP +#define BOOST_RANDOM_RANLUX_HPP + +#include +#include +#include + +namespace boost { + +namespace random { + typedef subtract_with_carry ranlux_base; + typedef subtract_with_carry_01 ranlux_base_01; + typedef subtract_with_carry_01 ranlux64_base_01; +} + +namespace random { +namespace detail { +/** + * The ranlux family of generators are described in + * + * @blockquote + * "A portable high-quality random number generator for lattice field theory + * calculations", M. Luescher, Computer Physics Communications, 79 (1994) + * pp 100-110. + * @endblockquote + * + * The levels are given in + * + * @blockquote + * "RANLUX: A Fortran implementation ofthe high-quality + * pseudorandom number generator of Luescher", F. James, + * Computer Physics Communications 79 (1994) 111-114 + * @endblockquote + */ +class ranlux_documentation {}; +} +} + +/** @copydoc boost::random::detail::ranlux_documentation */ +typedef random::discard_block ranlux3; +/** @copydoc boost::random::detail::ranlux_documentation */ +typedef random::discard_block ranlux4; + +/** @copydoc boost::random::detail::ranlux_documentation */ +typedef random::discard_block ranlux3_01; +/** @copydoc boost::random::detail::ranlux_documentation */ +typedef random::discard_block ranlux4_01; + +/** @copydoc boost::random::detail::ranlux_documentation */ +typedef random::discard_block ranlux64_3_01; +/** @copydoc boost::random::detail::ranlux_documentation */ +typedef random::discard_block ranlux64_4_01; + +#if !defined(BOOST_NO_INT64_T) && !defined(BOOST_NO_INTEGRAL_INT64_T) +namespace random { + typedef random::subtract_with_carry ranlux64_base; +} +/** @copydoc boost::random::detail::ranlux_documentation */ +typedef random::discard_block ranlux64_3; +/** @copydoc boost::random::detail::ranlux_documentation */ +typedef random::discard_block ranlux64_4; +#endif /* !BOOST_NO_INT64_T && !BOOST_NO_INTEGRAL_INT64_T */ + +} // namespace boost + +#endif // BOOST_RANDOM_LINEAR_CONGRUENTIAL_HPP diff --git a/Boost/boost/random/shuffle_output.hpp b/Boost/boost/random/shuffle_output.hpp new file mode 100644 index 000000000..9346a191d --- /dev/null +++ b/Boost/boost/random/shuffle_output.hpp @@ -0,0 +1,234 @@ +/* boost random/shuffle_output.hpp header file + * + * Copyright Jens Maurer 2000-2001 + * 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: shuffle_output.hpp 60755 2010-03-22 00:45:06Z steven_watanabe $ + * + * Revision history + * 2001-02-18 moved to individual header files + */ + +#ifndef BOOST_RANDOM_SHUFFLE_OUTPUT_HPP +#define BOOST_RANDOM_SHUFFLE_OUTPUT_HPP + +#include +#include // std::copy +#include +#include +#include +#include +#include +#include + +namespace boost { +namespace random { + +/** + * Instatiations of class template shuffle_output model a + * \pseudo_random_number_generator. It mixes the output + * of some (usually \linear_congruential) \uniform_random_number_generator + * to get better statistical properties. + * The algorithm is described in + * + * @blockquote + * "Improving a poor random number generator", Carter Bays + * and S.D. Durham, ACM Transactions on Mathematical Software, + * Vol 2, No. 1, March 1976, pp. 59-64. + * http://doi.acm.org/10.1145/355666.355670 + * @endblockquote + * + * The output of the base generator is buffered in an array of + * length k. Every output X(n) has a second role: It gives an + * index into the array where X(n+1) will be retrieved. Used + * array elements are replaced with fresh output from the base + * generator. + * + * Template parameters are the base generator and the array + * length k, which should be around 100. The template parameter + * val is the validation value checked by validation. + */ +template +class shuffle_output +{ +public: + typedef UniformRandomNumberGenerator base_type; + typedef typename base_type::result_type result_type; + + BOOST_STATIC_CONSTANT(bool, has_fixed_range = false); + BOOST_STATIC_CONSTANT(int, buffer_size = k); + + /** + * Constructs a @c shuffle_output generator by invoking the + * default constructor of the base generator. + * + * Complexity: Exactly k+1 invocations of the base generator. + */ + shuffle_output() : _rng() { init(); } +#if defined(BOOST_MSVC) && _MSC_VER < 1300 + // MSVC does not implicitly generate the copy constructor here + shuffle_output(const shuffle_output & x) + : _rng(x._rng), y(x.y) { std::copy(x.v, x.v+k, v); } +#endif + /** + * Constructs a shuffle_output generator by invoking the one-argument + * constructor of the base generator with the parameter seed. + * + * Complexity: Exactly k+1 invocations of the base generator. + */ + template + explicit shuffle_output(T s) : _rng(s) { init(); } + /** + * Constructs a shuffle_output generator by using a copy + * of the provided generator. + * + * Precondition: The template argument UniformRandomNumberGenerator + * shall denote a CopyConstructible type. + * + * Complexity: Exactly k+1 invocations of the base generator. + */ + explicit shuffle_output(const base_type & rng) : _rng(rng) { init(); } + template shuffle_output(It& first, It last) + : _rng(first, last) { init(); } + void seed() { _rng.seed(); init(); } + /** + * Invokes the one-argument seed method of the base generator + * with the parameter seed and re-initializes the internal buffer array. + * + * Complexity: Exactly k+1 invocations of the base generator. + */ + template + void seed(T s) { _rng.seed(s); init(); } + template void seed(It& first, It last) + { + _rng.seed(first, last); + init(); + } + + const base_type& base() const { return _rng; } + + result_type operator()() { + // calculating the range every time may seem wasteful. However, this + // makes the information locally available for the optimizer. + result_type range = (max)()-(min)()+1; + int j = k*(y-(min)())/range; + // assert(0 <= j && j < k); + y = v[j]; + v[j] = _rng(); + return y; + } + + result_type min BOOST_PREVENT_MACRO_SUBSTITUTION () const { return (_rng.min)(); } + result_type max BOOST_PREVENT_MACRO_SUBSTITUTION () const { return (_rng.max)(); } + static bool validation(result_type x) { return val == x; } + +#ifndef BOOST_NO_OPERATORS_IN_NAMESPACE + +#ifndef BOOST_RANDOM_NO_STREAM_OPERATORS + template + friend std::basic_ostream& + operator<<(std::basic_ostream& os, const shuffle_output& s) + { + os << s._rng << " " << s.y << " "; + for(int i = 0; i < s.buffer_size; ++i) + os << s.v[i] << " "; + return os; + } + + template + friend std::basic_istream& + operator>>(std::basic_istream& is, shuffle_output& s) + { + is >> s._rng >> std::ws >> s.y >> std::ws; + for(int i = 0; i < s.buffer_size; ++i) + is >> s.v[i] >> std::ws; + return is; + } +#endif + + friend bool operator==(const shuffle_output& x, const shuffle_output& y) + { return x._rng == y._rng && x.y == y.y && std::equal(x.v, x.v+k, y.v); } + friend bool operator!=(const shuffle_output& x, const shuffle_output& y) + { return !(x == y); } +#else + // Use a member function; Streamable concept not supported. + bool operator==(const shuffle_output& rhs) const + { return _rng == rhs._rng && y == rhs.y && std::equal(v, v+k, rhs.v); } + bool operator!=(const shuffle_output& rhs) const + { return !(*this == rhs); } +#endif +private: + void init() + { +#ifndef BOOST_NO_LIMITS_COMPILE_TIME_CONSTANTS + BOOST_STATIC_ASSERT(std::numeric_limits::is_integer); +#endif + result_type range = (max)()-(min)(); + assert(range > 0); // otherwise there would be little choice + if(static_cast(k * range) < + static_cast(range)) // not a sufficient condition + // likely overflow with bucket number computation + assert(!"overflow will occur"); + + // we cannot use std::generate, because it uses pass-by-value for _rng + for(result_type * p = v; p != v+k; ++p) + *p = _rng(); + y = _rng(); + } + + base_type _rng; + result_type v[k]; + result_type y; +}; + +#ifndef BOOST_NO_INCLASS_MEMBER_INITIALIZATION +// A definition is required even for integral static constants +template +const bool shuffle_output::has_fixed_range; + +template +const int shuffle_output::buffer_size; +#endif + +} // namespace random + +// validation by experiment from Harry Erwin's generator.h (private e-mail) +/** + * According to Harry Erwin (private e-mail), the specialization + * @c kreutzer1986 was suggested in: + * + * @blockquote + * "System Simulation: Programming Styles and Languages (International + * Computer Science Series)", Wolfgang Kreutzer, Addison-Wesley, December 1986. + * @endblockquote + */ +typedef random::shuffle_output< + random::linear_congruential, + 97, 139726> kreutzer1986; + + +} // namespace boost + +#endif // BOOST_RANDOM_SHUFFLE_OUTPUT_HPP diff --git a/Boost/boost/random/subtract_with_carry.hpp b/Boost/boost/random/subtract_with_carry.hpp new file mode 100644 index 000000000..bdb1401da --- /dev/null +++ b/Boost/boost/random/subtract_with_carry.hpp @@ -0,0 +1,465 @@ +/* boost random/subtract_with_carry.hpp header file + * + * Copyright Jens Maurer 2002 + * 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: subtract_with_carry.hpp 60755 2010-03-22 00:45:06Z steven_watanabe $ + * + * Revision history + * 2002-03-02 created + */ + +#ifndef BOOST_RANDOM_SUBTRACT_WITH_CARRY_HPP +#define BOOST_RANDOM_SUBTRACT_WITH_CARRY_HPP + +#include +#include +#include // std::equal +#include +#include // std::pow +#include +#include +#include +#include +#include +#include +#include +#include + + +namespace boost { +namespace random { + +#if BOOST_WORKAROUND(_MSC_FULL_VER, BOOST_TESTED_AT(13102292)) && BOOST_MSVC > 1300 +# define BOOST_RANDOM_EXTRACT_SWC_01 +#endif + +#if defined(__APPLE_CC__) && defined(__GNUC__) && (__GNUC__ == 3) && (__GNUC_MINOR__ <= 3) +# define BOOST_RANDOM_EXTRACT_SWC_01 +#endif + +# ifdef BOOST_RANDOM_EXTRACT_SWC_01 +namespace detail +{ + template + void extract_subtract_with_carry_01( + IStream& is + , SubtractWithCarry& f + , RealType& carry + , RealType* x + , RealType modulus) + { + RealType value; + for(unsigned int j = 0; j < f.long_lag; ++j) { + is >> value >> std::ws; + x[j] = value / modulus; + } + is >> value >> std::ws; + carry = value / modulus; + } +} +# endif + +/** + * Instantiations of @c subtract_with_carry model a + * \pseudo_random_number_generator. The algorithm is + * described in + * + * @blockquote + * "A New Class of Random Number Generators", George + * Marsaglia and Arif Zaman, Annals of Applied Probability, + * Volume 1, Number 3 (1991), 462-480. + * @endblockquote + */ +template +class subtract_with_carry +{ +public: + typedef IntType result_type; + BOOST_STATIC_CONSTANT(bool, has_fixed_range = true); + BOOST_STATIC_CONSTANT(result_type, min_value = 0); + BOOST_STATIC_CONSTANT(result_type, max_value = m-1); + BOOST_STATIC_CONSTANT(result_type, modulus = m); + BOOST_STATIC_CONSTANT(unsigned int, long_lag = r); + BOOST_STATIC_CONSTANT(unsigned int, short_lag = s); + + subtract_with_carry() { + // MSVC fails BOOST_STATIC_ASSERT with std::numeric_limits at class scope +#ifndef BOOST_NO_LIMITS_COMPILE_TIME_CONSTANTS + BOOST_STATIC_ASSERT(std::numeric_limits::is_signed); + BOOST_STATIC_ASSERT(std::numeric_limits::is_integer); +#endif + seed(); + } + BOOST_RANDOM_DETAIL_ARITHMETIC_CONSTRUCTOR(subtract_with_carry, uint32_t, value) + { seed(value); } + BOOST_RANDOM_DETAIL_GENERATOR_CONSTRUCTOR(subtract_with_carry, Generator, gen) + { seed(gen); } + template subtract_with_carry(It& first, It last) { seed(first,last); } + + // compiler-generated copy ctor and assignment operator are fine + + void seed() { seed(19780503u); } + BOOST_RANDOM_DETAIL_ARITHMETIC_SEED(subtract_with_carry, uint32_t, value) + { + random::linear_congruential intgen(value); + seed(intgen); + } + + // For GCC, moving this function out-of-line prevents inlining, which may + // reduce overall object code size. However, MSVC does not grok + // out-of-line template member functions. + BOOST_RANDOM_DETAIL_GENERATOR_SEED(subtract_with_carry, Generator, gen) + { + // I could have used std::generate_n, but it takes "gen" by value + for(unsigned int j = 0; j < long_lag; ++j) + x[j] = gen() % modulus; + carry = (x[long_lag-1] == 0); + k = 0; + } + + template + void seed(It& first, It last) + { + unsigned int j; + for(j = 0; j < long_lag && first != last; ++j, ++first) + x[j] = *first % modulus; + if(first == last && j < long_lag) + throw std::invalid_argument("subtract_with_carry::seed"); + carry = (x[long_lag-1] == 0); + k = 0; + } + + result_type min BOOST_PREVENT_MACRO_SUBSTITUTION () const { return min_value; } + result_type max BOOST_PREVENT_MACRO_SUBSTITUTION () const { return max_value; } + + result_type operator()() + { + int short_index = k - short_lag; + if(short_index < 0) + short_index += long_lag; + IntType delta; + if (x[short_index] >= x[k] + carry) { + // x(n) >= 0 + delta = x[short_index] - (x[k] + carry); + carry = 0; + } else { + // x(n) < 0 + delta = modulus - x[k] - carry + x[short_index]; + carry = 1; + } + x[k] = delta; + ++k; + if(k >= long_lag) + k = 0; + return delta; + } + +public: + static bool validation(result_type x) { return x == val; } + +#ifndef BOOST_NO_OPERATORS_IN_NAMESPACE + +#ifndef BOOST_RANDOM_NO_STREAM_OPERATORS + template + friend std::basic_ostream& + operator<<(std::basic_ostream& os, + const subtract_with_carry& f) + { + for(unsigned int j = 0; j < f.long_lag; ++j) + os << f.compute(j) << " "; + os << f.carry << " "; + return os; + } + + template + friend std::basic_istream& + operator>>(std::basic_istream& is, subtract_with_carry& f) + { + for(unsigned int j = 0; j < f.long_lag; ++j) + is >> f.x[j] >> std::ws; + is >> f.carry >> std::ws; + f.k = 0; + return is; + } +#endif + + friend bool operator==(const subtract_with_carry& x, const subtract_with_carry& y) + { + for(unsigned int j = 0; j < r; ++j) + if(x.compute(j) != y.compute(j)) + return false; + return true; + } + + friend bool operator!=(const subtract_with_carry& x, const subtract_with_carry& y) + { return !(x == y); } +#else + // Use a member function; Streamable concept not supported. + bool operator==(const subtract_with_carry& rhs) const + { + for(unsigned int j = 0; j < r; ++j) + if(compute(j) != rhs.compute(j)) + return false; + return true; + } + + bool operator!=(const subtract_with_carry& rhs) const + { return !(*this == rhs); } +#endif + +private: + /// \cond hide_private_members + // returns x(i-r+index), where index is in 0..r-1 + IntType compute(unsigned int index) const + { + return x[(k+index) % long_lag]; + } + /// \endcond + + // state representation; next output (state) is x(i) + // x[0] ... x[k] x[k+1] ... x[long_lag-1] represents + // x(i-k) ... x(i) x(i+1) ... x(i-k+long_lag-1) + // speed: base: 20-25 nsec + // ranlux_4: 230 nsec, ranlux_7: 430 nsec, ranlux_14: 810 nsec + // This state representation makes operator== and save/restore more + // difficult, because we've already computed "too much" and thus + // have to undo some steps to get at x(i-r) etc. + + // state representation: next output (state) is x(i) + // x[0] ... x[k] x[k+1] ... x[long_lag-1] represents + // x(i-k) ... x(i) x(i-long_lag+1) ... x(i-k-1) + // speed: base 28 nsec + // ranlux_4: 370 nsec, ranlux_7: 688 nsec, ranlux_14: 1343 nsec + IntType x[long_lag]; + unsigned int k; + int carry; +}; + +#ifndef BOOST_NO_INCLASS_MEMBER_INITIALIZATION +// A definition is required even for integral static constants +template +const bool subtract_with_carry::has_fixed_range; +template +const IntType subtract_with_carry::min_value; +template +const IntType subtract_with_carry::max_value; +template +const IntType subtract_with_carry::modulus; +template +const unsigned int subtract_with_carry::long_lag; +template +const unsigned int subtract_with_carry::short_lag; +#endif + + +// use a floating-point representation to produce values in [0..1) +/** @copydoc boost::random::subtract_with_carry */ +template +class subtract_with_carry_01 +{ +public: + typedef RealType result_type; + BOOST_STATIC_CONSTANT(bool, has_fixed_range = false); + BOOST_STATIC_CONSTANT(int, word_size = w); + BOOST_STATIC_CONSTANT(unsigned int, long_lag = r); + BOOST_STATIC_CONSTANT(unsigned int, short_lag = s); + +#ifndef BOOST_NO_LIMITS_COMPILE_TIME_CONSTANTS + BOOST_STATIC_ASSERT(!std::numeric_limits::is_integer); +#endif + + subtract_with_carry_01() { init_modulus(); seed(); } + explicit subtract_with_carry_01(uint32_t value) + { init_modulus(); seed(value); } + template subtract_with_carry_01(It& first, It last) + { init_modulus(); seed(first,last); } + +private: + /// \cond hide_private_members + void init_modulus() + { +#ifndef BOOST_NO_STDC_NAMESPACE + // allow for Koenig lookup + using std::pow; +#endif + _modulus = pow(RealType(2), word_size); + } + /// \endcond hide_private_members + +public: + // compiler-generated copy ctor and assignment operator are fine + + void seed(uint32_t value = 19780503u) + { +#ifndef BOOST_NO_STDC_NAMESPACE + // allow for Koenig lookup + using std::fmod; +#endif + random::linear_congruential gen(value); + unsigned long array[(w+31)/32 * long_lag]; + for(unsigned int j = 0; j < sizeof(array)/sizeof(unsigned long); ++j) + array[j] = gen(); + unsigned long * start = array; + seed(start, array + sizeof(array)/sizeof(unsigned long)); + } + + template + void seed(It& first, It last) + { +#ifndef BOOST_NO_STDC_NAMESPACE + // allow for Koenig lookup + using std::fmod; + using std::pow; +#endif + unsigned long mask = ~((~0u) << (w%32)); // now lowest (w%32) bits set + RealType two32 = pow(RealType(2), 32); + unsigned int j; + for(j = 0; j < long_lag && first != last; ++j) { + x[j] = RealType(0); + for(int i = 0; i < w/32 && first != last; ++i, ++first) + x[j] += *first / pow(two32,i+1); + if(first != last && mask != 0) { + x[j] += fmod((*first & mask) / _modulus, RealType(1)); + ++first; + } + } + if(first == last && j < long_lag) + throw std::invalid_argument("subtract_with_carry_01::seed"); + carry = (x[long_lag-1] ? 0 : 1 / _modulus); + k = 0; + } + + result_type min BOOST_PREVENT_MACRO_SUBSTITUTION () const { return result_type(0); } + result_type max BOOST_PREVENT_MACRO_SUBSTITUTION () const { return result_type(1); } + + result_type operator()() + { + int short_index = k - short_lag; + if(short_index < 0) + short_index += long_lag; + RealType delta = x[short_index] - x[k] - carry; + if(delta < 0) { + delta += RealType(1); + carry = RealType(1)/_modulus; + } else { + carry = 0; + } + x[k] = delta; + ++k; + if(k >= long_lag) + k = 0; + return delta; + } + + static bool validation(result_type x) + { return x == val/pow(RealType(2), word_size); } + +#ifndef BOOST_NO_OPERATORS_IN_NAMESPACE + +#ifndef BOOST_RANDOM_NO_STREAM_OPERATORS + template + friend std::basic_ostream& + operator<<(std::basic_ostream& os, + const subtract_with_carry_01& f) + { +#ifndef BOOST_NO_STDC_NAMESPACE + // allow for Koenig lookup + using std::pow; +#endif + std::ios_base::fmtflags oldflags = os.flags(os.dec | os.fixed | os.left); + for(unsigned int j = 0; j < f.long_lag; ++j) + os << (f.compute(j) * f._modulus) << " "; + os << (f.carry * f._modulus); + os.flags(oldflags); + return os; + } + + template + friend std::basic_istream& + operator>>(std::basic_istream& is, subtract_with_carry_01& f) + { +# ifdef BOOST_RANDOM_EXTRACT_SWC_01 + detail::extract_subtract_with_carry_01(is, f, f.carry, f.x, f._modulus); +# else + // MSVC (up to 7.1) and Borland (up to 5.64) don't handle the template type + // parameter "RealType" available from the class template scope, so use + // the member typedef + typename subtract_with_carry_01::result_type value; + for(unsigned int j = 0; j < long_lag; ++j) { + is >> value >> std::ws; + f.x[j] = value / f._modulus; + } + is >> value >> std::ws; + f.carry = value / f._modulus; +# endif + f.k = 0; + return is; + } +#endif + + friend bool operator==(const subtract_with_carry_01& x, + const subtract_with_carry_01& y) + { + for(unsigned int j = 0; j < r; ++j) + if(x.compute(j) != y.compute(j)) + return false; + return true; + } + + friend bool operator!=(const subtract_with_carry_01& x, + const subtract_with_carry_01& y) + { return !(x == y); } +#else + // Use a member function; Streamable concept not supported. + bool operator==(const subtract_with_carry_01& rhs) const + { + for(unsigned int j = 0; j < r; ++j) + if(compute(j) != rhs.compute(j)) + return false; + return true; + } + + bool operator!=(const subtract_with_carry_01& rhs) const + { return !(*this == rhs); } +#endif + +private: + /// \cond hide_private_members + RealType compute(unsigned int index) const; + /// \endcond + unsigned int k; + RealType carry; + RealType x[long_lag]; + RealType _modulus; +}; + +#ifndef BOOST_NO_INCLASS_MEMBER_INITIALIZATION +// A definition is required even for integral static constants +template +const bool subtract_with_carry_01::has_fixed_range; +template +const int subtract_with_carry_01::word_size; +template +const unsigned int subtract_with_carry_01::long_lag; +template +const unsigned int subtract_with_carry_01::short_lag; +#endif + +/// \cond hide_private_members +template +RealType subtract_with_carry_01::compute(unsigned int index) const +{ + return x[(k+index) % long_lag]; +} +/// \endcond + +} // namespace random +} // namespace boost + +#endif // BOOST_RANDOM_SUBTRACT_WITH_CARRY_HPP diff --git a/Boost/boost/random/triangle_distribution.hpp b/Boost/boost/random/triangle_distribution.hpp new file mode 100644 index 000000000..f4e139e00 --- /dev/null +++ b/Boost/boost/random/triangle_distribution.hpp @@ -0,0 +1,118 @@ +/* boost random/triangle_distribution.hpp header file + * + * Copyright Jens Maurer 2000-2001 + * 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: triangle_distribution.hpp 60755 2010-03-22 00:45:06Z steven_watanabe $ + * + * Revision history + * 2001-02-18 moved to individual header files + */ + +#ifndef BOOST_RANDOM_TRIANGLE_DISTRIBUTION_HPP +#define BOOST_RANDOM_TRIANGLE_DISTRIBUTION_HPP + +#include +#include +#include +#include + +namespace boost { + +/** + * Instantiations of @c triangle_distribution model a \random_distribution. + * A @c triangle_distribution has three parameters, @c a, @c b, and @c c, + * which are the smallest, the most probable and the largest values of + * the distribution respectively. + */ +template +class triangle_distribution +{ +public: + typedef RealType input_type; + typedef RealType result_type; + + /** + * Constructs a @c triangle_distribution with the parameters + * @c a, @c b, and @c c. + * + * Preconditions: a <= b <= c. + */ + explicit triangle_distribution(result_type a_arg = result_type(0), + result_type b_arg = result_type(0.5), + result_type c_arg = result_type(1)) + : _a(a_arg), _b(b_arg), _c(c_arg) + { + assert(_a <= _b && _b <= _c); + init(); + } + + // compiler-generated copy ctor and assignment operator are fine + + /** Returns the @c a parameter of the distribution */ + result_type a() const { return _a; } + /** Returns the @c b parameter of the distribution */ + result_type b() const { return _b; } + /** Returns the @c c parameter of the distribution */ + result_type c() const { return _c; } + + void reset() { } + + template + result_type operator()(Engine& eng) + { +#ifndef BOOST_NO_STDC_NAMESPACE + using std::sqrt; +#endif + result_type u = eng(); + if( u <= q1 ) + return _a + p1*sqrt(u); + else + return _c - d3*sqrt(d2*u-d1); + } + +#ifndef BOOST_RANDOM_NO_STREAM_OPERATORS + template + friend std::basic_ostream& + operator<<(std::basic_ostream& os, const triangle_distribution& td) + { + os << td._a << " " << td._b << " " << td._c; + return os; + } + + template + friend std::basic_istream& + operator>>(std::basic_istream& is, triangle_distribution& td) + { + is >> std::ws >> td._a >> std::ws >> td._b >> std::ws >> td._c; + td.init(); + return is; + } +#endif + +private: + /// \cond hide_private_members + void init() + { +#ifndef BOOST_NO_STDC_NAMESPACE + using std::sqrt; +#endif + d1 = _b - _a; + d2 = _c - _a; + d3 = sqrt(_c - _b); + q1 = d1 / d2; + p1 = sqrt(d1 * d2); + } + /// \endcond + + result_type _a, _b, _c; + result_type d1, d2, d3, q1, p1; +}; + +} // namespace boost + +#endif // BOOST_RANDOM_TRIANGLE_DISTRIBUTION_HPP diff --git a/Boost/boost/random/uniform_01.hpp b/Boost/boost/random/uniform_01.hpp new file mode 100644 index 000000000..87103cb12 --- /dev/null +++ b/Boost/boost/random/uniform_01.hpp @@ -0,0 +1,273 @@ +/* boost random/uniform_01.hpp header file + * + * Copyright Jens Maurer 2000-2001 + * 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: uniform_01.hpp 60755 2010-03-22 00:45:06Z steven_watanabe $ + * + * Revision history + * 2001-02-18 moved to individual header files + */ + +#ifndef BOOST_RANDOM_UNIFORM_01_HPP +#define BOOST_RANDOM_UNIFORM_01_HPP + +#include +#include +#include +#include +#include +#include + +#include + +namespace boost { + +#ifdef BOOST_RANDOM_DOXYGEN + +/** + * The distribution function uniform_01 models a \random_distribution. + * On each invocation, it returns a random floating-point value + * uniformly distributed in the range [0..1). + * + * The template parameter RealType shall denote a float-like value type + * with support for binary operators +, -, and /. + * + * Note: The current implementation is buggy, because it may not fill + * all of the mantissa with random bits. I'm unsure how to fill a + * (to-be-invented) @c boost::bigfloat class with random bits efficiently. + * It's probably time for a traits class. + */ +template +class uniform_01 +{ +public: + typedef RealType input_type; + typedef RealType result_type; + result_type min BOOST_PREVENT_MACRO_SUBSTITUTION () const; + result_type max BOOST_PREVENT_MACRO_SUBSTITUTION () const; + void reset(); + + template + result_type operator()(Engine& eng); + +#ifndef BOOST_RANDOM_NO_STREAM_OPERATORS + template + friend std::basic_ostream& + operator<<(std::basic_ostream& os, const new_uniform_01&) + { + return os; + } + + template + friend std::basic_istream& + operator>>(std::basic_istream& is, new_uniform_01&) + { + return is; + } +#endif +}; + +#else + +namespace detail { + +template +class new_uniform_01 +{ +public: + typedef RealType input_type; + typedef RealType result_type; + // compiler-generated copy ctor and copy assignment are fine + result_type min BOOST_PREVENT_MACRO_SUBSTITUTION () const { return result_type(0); } + result_type max BOOST_PREVENT_MACRO_SUBSTITUTION () const { return result_type(1); } + void reset() { } + + template + result_type operator()(Engine& eng) { + for (;;) { + typedef typename Engine::result_type base_result; + result_type factor = result_type(1) / + (result_type((eng.max)()-(eng.min)()) + + result_type(std::numeric_limits::is_integer ? 1 : 0)); + result_type result = result_type(eng() - (eng.min)()) * factor; + if (result < result_type(1)) + return result; + } + } + +#ifndef BOOST_RANDOM_NO_STREAM_OPERATORS + template + friend std::basic_ostream& + operator<<(std::basic_ostream& os, const new_uniform_01&) + { + return os; + } + + template + friend std::basic_istream& + operator>>(std::basic_istream& is, new_uniform_01&) + { + return is; + } +#endif +}; + +template +class backward_compatible_uniform_01 +{ + typedef boost::random::detail::ptr_helper traits; + typedef boost::random::detail::pass_through_engine internal_engine_type; +public: + typedef UniformRandomNumberGenerator base_type; + typedef RealType result_type; + + BOOST_STATIC_CONSTANT(bool, has_fixed_range = false); + +#if !defined(BOOST_NO_LIMITS_COMPILE_TIME_CONSTANTS) && !(defined(BOOST_MSVC) && BOOST_MSVC <= 1300) + BOOST_STATIC_ASSERT(!std::numeric_limits::is_integer); +#endif + + explicit backward_compatible_uniform_01(typename traits::rvalue_type rng) + : _rng(rng), + _factor(result_type(1) / + (result_type((_rng.max)()-(_rng.min)()) + + result_type(std::numeric_limits::is_integer ? 1 : 0))) + { + } + // compiler-generated copy ctor and copy assignment are fine + + result_type min BOOST_PREVENT_MACRO_SUBSTITUTION () const { return result_type(0); } + result_type max BOOST_PREVENT_MACRO_SUBSTITUTION () const { return result_type(1); } + typename traits::value_type& base() { return _rng.base(); } + const typename traits::value_type& base() const { return _rng.base(); } + void reset() { } + + result_type operator()() { + for (;;) { + result_type result = result_type(_rng() - (_rng.min)()) * _factor; + if (result < result_type(1)) + return result; + } + } + +#if !defined(BOOST_NO_OPERATORS_IN_NAMESPACE) && !defined(BOOST_NO_MEMBER_TEMPLATE_FRIENDS) + template + friend std::basic_ostream& + operator<<(std::basic_ostream& os, const backward_compatible_uniform_01& u) + { + os << u._rng; + return os; + } + + template + friend std::basic_istream& + operator>>(std::basic_istream& is, backward_compatible_uniform_01& u) + { + is >> u._rng; + return is; + } +#endif + +private: + typedef typename internal_engine_type::result_type base_result; + internal_engine_type _rng; + result_type _factor; +}; + +#ifndef BOOST_NO_INCLASS_MEMBER_INITIALIZATION +// A definition is required even for integral static constants +template +const bool backward_compatible_uniform_01::has_fixed_range; +#endif + +template +struct select_uniform_01 +{ + template + struct apply + { + typedef backward_compatible_uniform_01 type; + }; +}; + +template<> +struct select_uniform_01 +{ + template + struct apply + { + typedef new_uniform_01 type; + }; +}; + +template<> +struct select_uniform_01 +{ + template + struct apply + { + typedef new_uniform_01 type; + }; +}; + +template<> +struct select_uniform_01 +{ + template + struct apply + { + typedef new_uniform_01 type; + }; +}; + +} + +// Because it is so commonly used: uniform distribution on the real [0..1) +// range. This allows for specializations to avoid a costly int -> float +// conversion plus float multiplication +template +class uniform_01 + : public detail::select_uniform_01::BOOST_NESTED_TEMPLATE apply::type +{ + typedef typename detail::select_uniform_01::BOOST_NESTED_TEMPLATE apply::type impl_type; + typedef boost::random::detail::ptr_helper traits; +public: + + uniform_01() {} + + explicit uniform_01(typename traits::rvalue_type rng) + : impl_type(rng) + { + } + +#if !defined(BOOST_NO_OPERATORS_IN_NAMESPACE) && !defined(BOOST_NO_MEMBER_TEMPLATE_FRIENDS) + template + friend std::basic_ostream& + operator<<(std::basic_ostream& os, const uniform_01& u) + { + os << static_cast(u); + return os; + } + + template + friend std::basic_istream& + operator>>(std::basic_istream& is, uniform_01& u) + { + is >> static_cast(u); + return is; + } +#endif +}; + +#endif + +} // namespace boost + +#include + +#endif // BOOST_RANDOM_UNIFORM_01_HPP diff --git a/Boost/boost/random/uniform_int.hpp b/Boost/boost/random/uniform_int.hpp new file mode 100644 index 000000000..ef2f4f96a --- /dev/null +++ b/Boost/boost/random/uniform_int.hpp @@ -0,0 +1,300 @@ +/* boost random/uniform_int.hpp header file + * + * Copyright Jens Maurer 2000-2001 + * 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: uniform_int.hpp 60755 2010-03-22 00:45:06Z steven_watanabe $ + * + * Revision history + * 2001-04-08 added min +#include +#include +#include +#include +#include +#include +#include +#include + +namespace boost { + +/** + * The distribution function uniform_int models a \random_distribution. + * On each invocation, it returns a random integer value uniformly + * distributed in the set of integer numbers {min, min+1, min+2, ..., max}. + * + * The template parameter IntType shall denote an integer-like value type. + */ +template +class uniform_int +{ +public: + typedef IntType input_type; + typedef IntType result_type; + + /// \cond hide_private_members + typedef typename make_unsigned::type range_type; + /// \endcond + + /** + * Constructs a uniform_int object. @c min and @c max are + * the parameters of the distribution. + * + * Requires: min <= max + */ + explicit uniform_int(IntType min_arg = 0, IntType max_arg = 9) + : _min(min_arg), _max(max_arg) + { +#ifndef BOOST_NO_LIMITS_COMPILE_TIME_CONSTANTS + // MSVC fails BOOST_STATIC_ASSERT with std::numeric_limits at class scope + BOOST_STATIC_ASSERT(std::numeric_limits::is_integer); +#endif + assert(min_arg <= max_arg); + init(); + } + + /** + * Returns: The "min" parameter of the distribution + */ + result_type min BOOST_PREVENT_MACRO_SUBSTITUTION () const { return _min; } + /** + * Returns: The "max" parameter of the distribution + */ + result_type max BOOST_PREVENT_MACRO_SUBSTITUTION () const { return _max; } + void reset() { } + + // can't have member function templates out-of-line due to MSVC bugs + template + result_type operator()(Engine& eng) + { + return generate(eng, _min, _max, _range); + } + + template + result_type operator()(Engine& eng, result_type n) + { + assert(n > 0); + + if (n == 1) + { + return 0; + } + + return generate(eng, 0, n - 1, n - 1); + } + +#ifndef BOOST_RANDOM_NO_STREAM_OPERATORS + template + friend std::basic_ostream& + operator<<(std::basic_ostream& os, const uniform_int& ud) + { + os << ud._min << " " << ud._max; + return os; + } + + template + friend std::basic_istream& + operator>>(std::basic_istream& is, uniform_int& ud) + { + is >> std::ws >> ud._min >> std::ws >> ud._max; + ud.init(); + return is; + } +#endif + +private: + +#ifdef BOOST_MSVC +#pragma warning(push) +// disable division by zero warning, since we can't +// actually divide by zero. +#pragma warning(disable:4723) +#endif + + /// \cond hide_private_members + template + static result_type generate(Engine& eng, result_type min_value, result_type /*max_value*/, range_type range) + { + typedef typename Engine::result_type base_result; + // ranges are always unsigned + typedef typename make_unsigned::type base_unsigned; + const base_result bmin = (eng.min)(); + const base_unsigned brange = + random::detail::subtract()((eng.max)(), (eng.min)()); + + if(range == 0) { + return min_value; + } else if(brange == range) { + // this will probably never happen in real life + // basically nothing to do; just take care we don't overflow / underflow + base_unsigned v = random::detail::subtract()(eng(), bmin); + return random::detail::add()(v, min_value); + } else if(brange < range) { + // use rejection method to handle things like 0..3 --> 0..4 + for(;;) { + // concatenate several invocations of the base RNG + // take extra care to avoid overflows + + // limit == floor((range+1)/(brange+1)) + // Therefore limit*(brange+1) <= range+1 + range_type limit; + if(range == (std::numeric_limits::max)()) { + limit = range/(range_type(brange)+1); + if(range % (range_type(brange)+1) == range_type(brange)) + ++limit; + } else { + limit = (range+1)/(range_type(brange)+1); + } + + // We consider "result" as expressed to base (brange+1): + // For every power of (brange+1), we determine a random factor + range_type result = range_type(0); + range_type mult = range_type(1); + + // loop invariants: + // result < mult + // mult <= range + while(mult <= limit) { + // Postcondition: result <= range, thus no overflow + // + // limit*(brange+1)<=range+1 def. of limit (1) + // eng()-bmin<=brange eng() post. (2) + // and mult<=limit. loop condition (3) + // Therefore mult*(eng()-bmin+1)<=range+1 by (1),(2),(3) (4) + // Therefore mult*(eng()-bmin)+mult<=range+1 rearranging (4) (5) + // result(random::detail::subtract()(eng(), bmin) * mult); + + // equivalent to (mult * (brange+1)) == range+1, but avoids overflow. + if(mult * range_type(brange) == range - mult + 1) { + // The destination range is an integer power of + // the generator's range. + return(result); + } + + // Postcondition: mult <= range + // + // limit*(brange+1)<=range+1 def. of limit (1) + // mult<=limit loop condition (2) + // Therefore mult*(brange+1)<=range+1 by (1), (2) (3) + // mult*(brange+1)!=range+1 preceding if (4) + // Therefore mult*(brange+1) limit loop condition (1) + // Suppose range/mult >= brange+1 Assumption (2) + // range >= mult*(brange+1) by (2) (3) + // range+1 > mult*(brange+1) by (3) (4) + // range+1 > (limit+1)*(brange+1) by (1), (4) (5) + // (range+1)/(brange+1) > limit+1 by (5) (6) + // limit < floor((range+1)/(brange+1)) by (6) (7) + // limit==floor((range+1)/(brange+1)) def. of limit (8) + // not (2) reductio (9) + // + // loop postcondition: (range/mult)*mult+(mult-1) >= range + // + // (range/mult)*mult + range%mult == range identity (1) + // range%mult < mult def. of % (2) + // (range/mult)*mult+mult > range by (1), (2) (3) + // (range/mult)*mult+(mult-1) >= range by (3) (4) + // + // Note that the maximum value of result at this point is (mult-1), + // so after this final step, we generate numbers that can be + // at least as large as range. We have to really careful to avoid + // overflow in this final addition and in the rejection. Anything + // that overflows is larger than range and can thus be rejected. + + // range/mult < brange+1 -> no endless loop + range_type result_increment = uniform_int(0, range/mult)(eng); + if((std::numeric_limits::max)() / mult < result_increment) { + // The multiplcation would overflow. Reject immediately. + continue; + } + result_increment *= mult; + // unsigned integers are guaranteed to wrap on overflow. + result += result_increment; + if(result < result_increment) { + // The addition overflowed. Reject. + continue; + } + if(result > range) { + // Too big. Reject. + continue; + } + return random::detail::add()(result, min_value); + } + } else { // brange > range + base_unsigned bucket_size; + // it's safe to add 1 to range, as long as we cast it first, + // because we know that it is less than brange. However, + // we do need to be careful not to cause overflow by adding 1 + // to brange. + if(brange == (std::numeric_limits::max)()) { + bucket_size = brange / (static_cast(range)+1); + if(brange % (static_cast(range)+1) == static_cast(range)) { + ++bucket_size; + } + } else { + bucket_size = (brange+1) / (static_cast(range)+1); + } + for(;;) { + base_unsigned result = + random::detail::subtract()(eng(), bmin); + result /= bucket_size; + // result and range are non-negative, and result is possibly larger + // than range, so the cast is safe + if(result <= static_cast(range)) + return random::detail::add()(result, min_value); + } + } + } + +#ifdef BOOST_MSVC +#pragma warning(pop) +#endif + + void init() + { + _range = random::detail::subtract()(_max, _min); + } + + /// \endcond + + // The result_type may be signed or unsigned, but the _range is always + // unsigned. + result_type _min, _max; + range_type _range; +}; + +} // namespace boost + +#endif // BOOST_RANDOM_UNIFORM_INT_HPP diff --git a/Boost/boost/random/uniform_on_sphere.hpp b/Boost/boost/random/uniform_on_sphere.hpp new file mode 100644 index 000000000..5dae3439a --- /dev/null +++ b/Boost/boost/random/uniform_on_sphere.hpp @@ -0,0 +1,101 @@ +/* boost random/uniform_on_sphere.hpp header file + * + * Copyright Jens Maurer 2000-2001 + * 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: uniform_on_sphere.hpp 60755 2010-03-22 00:45:06Z steven_watanabe $ + * + * Revision history + * 2001-02-18 moved to individual header files + */ + +#ifndef BOOST_RANDOM_UNIFORM_ON_SPHERE_HPP +#define BOOST_RANDOM_UNIFORM_ON_SPHERE_HPP + +#include +#include // std::transform +#include // std::bind2nd, std::divides +#include +#include + +namespace boost { + +/** + * Instantiations of class template uniform_on_sphere model a + * \random_distribution. Such a distribution produces random + * numbers uniformly distributed on the unit sphere of arbitrary + * dimension @c dim. The @c Cont template parameter must be a STL-like + * container type with begin and end operations returning non-const + * ForwardIterators of type @c Cont::iterator. Each invocation of the + * @c UniformRandomNumberGenerator shall result in a floating-point + * value in the range [0,1). + */ +template > +class uniform_on_sphere +{ +public: + typedef RealType input_type; + typedef Cont result_type; + + /** + * Constructs a @c uniform_on_sphere distribution. + * @c dim is the dimension of the sphere. + */ + explicit uniform_on_sphere(int dim = 2) : _container(dim), _dim(dim) { } + + // compiler-generated copy ctor and assignment operator are fine + + void reset() { _normal.reset(); } + + template + const result_type & operator()(Engine& eng) + { + RealType sqsum = 0; + for(typename Cont::iterator it = _container.begin(); + it != _container.end(); + ++it) { + RealType val = _normal(eng); + *it = val; + sqsum += val * val; + } +#ifndef BOOST_NO_STDC_NAMESPACE + using std::sqrt; +#endif + // for all i: result[i] /= sqrt(sqsum) + std::transform(_container.begin(), _container.end(), _container.begin(), + std::bind2nd(std::divides(), sqrt(sqsum))); + return _container; + } + +#ifndef BOOST_RANDOM_NO_STREAM_OPERATORS + template + friend std::basic_ostream& + operator<<(std::basic_ostream& os, const uniform_on_sphere& sd) + { + os << sd._dim; + return os; + } + + template + friend std::basic_istream& + operator>>(std::basic_istream& is, uniform_on_sphere& sd) + { + is >> std::ws >> sd._dim; + sd._container.resize(sd._dim); + return is; + } +#endif + +private: + normal_distribution _normal; + result_type _container; + int _dim; +}; + +} // namespace boost + +#endif // BOOST_RANDOM_UNIFORM_ON_SPHERE_HPP diff --git a/Boost/boost/random/uniform_real.hpp b/Boost/boost/random/uniform_real.hpp new file mode 100644 index 000000000..0a04540a6 --- /dev/null +++ b/Boost/boost/random/uniform_real.hpp @@ -0,0 +1,108 @@ +/* boost random/uniform_real.hpp header file + * + * Copyright Jens Maurer 2000-2001 + * 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: uniform_real.hpp 60755 2010-03-22 00:45:06Z steven_watanabe $ + * + * Revision history + * 2001-04-08 added min +#include +#include +#include +#include +#include + +namespace boost { + +/** + * The distribution function uniform_real models a random distribution. + * On each invocation, it returns a random floating-point value uniformly + * distributed in the range [min..max). The value is computed using + * std::numeric_limits::digits random binary digits, i.e. + * the mantissa of the floating-point value is completely filled with + * random bits. + * + * Note: The current implementation is buggy, because it may not fill + * all of the mantissa with random bits. + */ +template +class uniform_real +{ +public: + typedef RealType input_type; + typedef RealType result_type; + + /** + * Constructs a uniform_real object. @c min and @c max are the + * parameters of the distribution. + * + * Requires: min <= max + */ + explicit uniform_real(RealType min_arg = RealType(0), + RealType max_arg = RealType(1)) + : _min(min_arg), _max(max_arg) + { +#ifndef BOOST_NO_LIMITS_COMPILE_TIME_CONSTANTS + BOOST_STATIC_ASSERT(!std::numeric_limits::is_integer); +#endif + assert(min_arg <= max_arg); + } + + // compiler-generated copy ctor and assignment operator are fine + + /** + * Returns: The "min" parameter of the distribution + */ + result_type min BOOST_PREVENT_MACRO_SUBSTITUTION () const { return _min; } + /** + * Returns: The "max" parameter of the distribution + */ + result_type max BOOST_PREVENT_MACRO_SUBSTITUTION () const { return _max; } + void reset() { } + + template + result_type operator()(Engine& eng) { + result_type numerator = static_cast(eng() - eng.min BOOST_PREVENT_MACRO_SUBSTITUTION()); + result_type divisor = static_cast(eng.max BOOST_PREVENT_MACRO_SUBSTITUTION() - eng.min BOOST_PREVENT_MACRO_SUBSTITUTION()); + assert(divisor > 0); + assert(numerator >= 0 && numerator <= divisor); + return numerator / divisor * (_max - _min) + _min; + } + +#ifndef BOOST_RANDOM_NO_STREAM_OPERATORS + template + friend std::basic_ostream& + operator<<(std::basic_ostream& os, const uniform_real& ud) + { + os << ud._min << " " << ud._max; + return os; + } + + template + friend std::basic_istream& + operator>>(std::basic_istream& is, uniform_real& ud) + { + is >> std::ws >> ud._min >> std::ws >> ud._max; + return is; + } +#endif + +private: + RealType _min, _max; +}; + +} // namespace boost + +#endif // BOOST_RANDOM_UNIFORM_REAL_HPP diff --git a/Boost/boost/random/uniform_smallint.hpp b/Boost/boost/random/uniform_smallint.hpp new file mode 100644 index 000000000..93c7b63dc --- /dev/null +++ b/Boost/boost/random/uniform_smallint.hpp @@ -0,0 +1,158 @@ +/* boost random/uniform_smallint.hpp header file + * + * Copyright Jens Maurer 2000-2001 + * 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: uniform_smallint.hpp 60755 2010-03-22 00:45:06Z steven_watanabe $ + * + * Revision history + * 2001-04-08 added min +#include +#include +#include +#include +#include +#include +#include + +namespace boost { + +// uniform integer distribution on a small range [min, max] + +/** + * The distribution function uniform_smallint models a \random_distribution. + * On each invocation, it returns a random integer value uniformly distributed + * in the set of integer numbers {min, min+1, min+2, ..., max}. It assumes + * that the desired range (max-min+1) is small compared to the range of the + * underlying source of random numbers and thus makes no attempt to limit + * quantization errors. + * + * Let rout=(max-min+1) the desired range of integer numbers, and + * let rbase be the range of the underlying source of random + * numbers. Then, for the uniform distribution, the theoretical probability + * for any number i in the range rout will be pout(i) = + * 1/rout. Likewise, assume a uniform distribution on rbase for + * the underlying source of random numbers, i.e. pbase(i) = + * 1/rbase. Let pout_s(i) denote the random + * distribution generated by @c uniform_smallint. Then the sum over all + * i in rout of (pout_s(i)/pout(i) - 1)2 + * shall not exceed rout/rbase2 + * (rbase mod rout)(rout - + * rbase mod rout). + * + * The template parameter IntType shall denote an integer-like value type. + * + * Note: The property above is the square sum of the relative differences + * in probabilities between the desired uniform distribution + * pout(i) and the generated distribution pout_s(i). + * The property can be fulfilled with the calculation + * (base_rng mod rout), as follows: Let r = rbase mod + * rout. The base distribution on rbase is folded onto the + * range rout. The numbers i < r have assigned (rbase + * div rout)+1 numbers of the base distribution, the rest has + * only (rbase div rout). Therefore, + * pout_s(i) = ((rbase div rout)+1) / + * rbase for i < r and pout_s(i) = (rbase + * div rout)/rbase otherwise. Substituting this in the + * above sum formula leads to the desired result. + * + * Note: The upper bound for (rbase mod rout) + * (rout - rbase mod rout) is + * rout2/4. Regarding the upper bound for the + * square sum of the relative quantization error of + * rout3/(4*rbase2), it + * seems wise to either choose rbase so that rbase > + * 10*rout2 or ensure that rbase is + * divisible by rout. + */ +template +class uniform_smallint +{ +public: + typedef IntType input_type; + typedef IntType result_type; + + /** + * Constructs a @c uniform_smallint. @c min and @c max are the + * lower and upper bounds of the output range, respectively. + */ + explicit uniform_smallint(IntType min_arg = 0, IntType max_arg = 9) + : _min(min_arg), _max(max_arg) + { +#ifndef BOOST_NO_LIMITS_COMPILE_TIME_CONSTANTS + // MSVC fails BOOST_STATIC_ASSERT with std::numeric_limits at class scope + BOOST_STATIC_ASSERT(std::numeric_limits::is_integer); +#endif + } + + result_type min BOOST_PREVENT_MACRO_SUBSTITUTION () const { return _min; } + result_type max BOOST_PREVENT_MACRO_SUBSTITUTION () const { return _max; } + void reset() { } + + template + result_type operator()(Engine& eng) + { + typedef typename Engine::result_type base_result; + base_result _range = static_cast(_max-_min)+1; + base_result _factor = 1; + + // LCGs get bad when only taking the low bits. + // (probably put this logic into a partial template specialization) + // Check how many low bits we can ignore before we get too much + // quantization error. + base_result r_base = (eng.max)() - (eng.min)(); + if(r_base == (std::numeric_limits::max)()) { + _factor = 2; + r_base /= 2; + } + r_base += 1; + if(r_base % _range == 0) { + // No quantization effects, good + _factor = r_base / _range; + } else { + // carefully avoid overflow; pessimizing here + for( ; r_base/_range/32 >= _range; _factor *= 2) + r_base /= 2; + } + + return static_cast(((eng() - (eng.min)()) / _factor) % _range + _min); + } + +#ifndef BOOST_RANDOM_NO_STREAM_OPERATORS + template + friend std::basic_ostream& + operator<<(std::basic_ostream& os, const uniform_smallint& ud) + { + os << ud._min << " " << ud._max; + return os; + } + + template + friend std::basic_istream& + operator>>(std::basic_istream& is, uniform_smallint& ud) + { + is >> std::ws >> ud._min >> std::ws >> ud._max; + return is; + } +#endif + +private: + + result_type _min; + result_type _max; +}; + +} // namespace boost + +#endif // BOOST_RANDOM_UNIFORM_SMALLINT_HPP diff --git a/Boost/boost/random/variate_generator.hpp b/Boost/boost/random/variate_generator.hpp new file mode 100644 index 000000000..92c44e6ef --- /dev/null +++ b/Boost/boost/random/variate_generator.hpp @@ -0,0 +1,220 @@ +/* boost random/variate_generator.hpp header file + * + * Copyright Jens Maurer 2002 + * 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: variate_generator.hpp 60755 2010-03-22 00:45:06Z steven_watanabe $ + * + */ + +#ifndef BOOST_RANDOM_RANDOM_GENERATOR_HPP +#define BOOST_RANDOM_RANDOM_GENERATOR_HPP + +#include + +// implementation details +#include +#include +#include +#include +#include + +// Borland C++ 5.6.0 has problems using its numeric_limits traits as +// template parameters +#if BOOST_WORKAROUND(__BORLANDC__, <= 0x564) +#include +#endif + +#include + +namespace boost { + +/// \cond hide_private_members + +namespace random { +namespace detail { + +template +struct engine_helper; + +// for consistency, always have two levels of decorations +template<> +struct engine_helper +{ + template + struct impl + { + typedef pass_through_engine type; + }; +}; + +template<> +struct engine_helper +{ + template + struct impl + { + typedef uniform_01 type; + }; +}; + +template<> +struct engine_helper +{ + template + struct impl + { + typedef uniform_01 type; + }; +}; + +template<> +struct engine_helper +{ + template + struct impl + { + typedef uniform_int_float type; + }; +}; + +} // namespace detail +} // namespace random + +///\endcond + +/** + * A random variate generator is used to join a random number + * generator together with a random number distribution. + * Boost.Random provides a vast choice of \generators as well + * as \distributions. + * + * Instantations of class template @c variate_generator model + * a \number_generator. + * + * The argument for the template parameter Engine shall be of + * the form U, U&, or U*, where U models a + * \uniform_random_number_generator. Then, the member + * engine_value_type names U (not the pointer or reference to U). + * + * Specializations of @c variate_generator satisfy the + * requirements of CopyConstructible. They also satisfy the + * requirements of Assignable unless the template parameter + * Engine is of the form U&. + * + * The complexity of all functions specified in this section + * is constant. No function described in this section except + * the constructor throws an exception. + */ +template +class variate_generator +{ +private: + typedef random::detail::pass_through_engine decorated_engine; + +public: + typedef typename decorated_engine::base_type engine_value_type; + typedef Engine engine_type; + typedef Distribution distribution_type; + typedef typename Distribution::result_type result_type; + + /** + * Constructs a @c variate_generator object with the associated + * \uniform_random_number_generator eng and the associated + * \random_distribution d. + * + * Throws: If and what the copy constructor of Engine or + * Distribution throws. + */ + variate_generator(Engine e, Distribution d) + : _eng(decorated_engine(e)), _dist(d) { } + + /** + * Returns: distribution()(e) + * + * Notes: The sequence of numbers produced by the + * \uniform_random_number_generator e, se, is + * obtained from the sequence of numbers produced by the + * associated \uniform_random_number_generator eng, seng, + * as follows: Consider the values of @c numeric_limits::is_integer + * for @c T both @c Distribution::input_type and + * @c engine_value_type::result_type. If the values for both types are + * true, then se is identical to seng. Otherwise, if the + * values for both types are false, then the numbers in seng + * are divided by engine().max()-engine().min() to obtain the numbers + * in se. Otherwise, if the value for + * @c engine_value_type::result_type is true and the value for + * @c Distribution::input_type is false, then the numbers in seng + * are divided by engine().max()-engine().min()+1 to obtain the numbers in + * se. Otherwise, the mapping from seng to + * se is implementation-defined. In all cases, an + * implicit conversion from @c engine_value_type::result_type to + * @c Distribution::input_type is performed. If such a conversion does + * not exist, the program is ill-formed. + */ + result_type operator()() { return _dist(_eng); } + /** + * Returns: distribution()(e, value). + * For the semantics of e, see the description of operator()(). + */ + template + result_type operator()(T value) { return _dist(_eng, value); } + + /** + * Returns: A reference to the associated uniform random number generator. + */ + engine_value_type& engine() { return _eng.base().base(); } + /** + * Returns: A reference to the associated uniform random number generator. + */ + const engine_value_type& engine() const { return _eng.base().base(); } + + /** + * Returns: A reference to the associated random distribution. + */ + distribution_type& distribution() { return _dist; } + /** + * Returns: A reference to the associated random distribution. + */ + const distribution_type& distribution() const { return _dist; } + + /** + * Precondition: distribution().min() is well-formed + * + * Returns: distribution().min() + */ + result_type min BOOST_PREVENT_MACRO_SUBSTITUTION () const { return (distribution().min)(); } + /** + * Precondition: distribution().max() is well-formed + * + * Returns: distribution().max() + */ + result_type max BOOST_PREVENT_MACRO_SUBSTITUTION () const { return (distribution().max)(); } + +private: +#if BOOST_WORKAROUND(__BORLANDC__, <= 0x564) + typedef typename random::detail::engine_helper< + ::boost::is_integral::value, + ::boost::is_integral::value + >::BOOST_NESTED_TEMPLATE impl::type internal_engine_type; +#else + enum { + have_int = std::numeric_limits::is_integer, + want_int = std::numeric_limits::is_integer + }; + typedef typename random::detail::engine_helper::BOOST_NESTED_TEMPLATE impl::type internal_engine_type; +#endif + + internal_engine_type _eng; + distribution_type _dist; +}; + +} // namespace boost + +#include + +#endif // BOOST_RANDOM_RANDOM_GENERATOR_HPP diff --git a/Boost/boost/random/xor_combine.hpp b/Boost/boost/random/xor_combine.hpp new file mode 100644 index 000000000..a8e991723 --- /dev/null +++ b/Boost/boost/random/xor_combine.hpp @@ -0,0 +1,184 @@ +/* boost random/xor_combine.hpp header file + * + * Copyright Jens Maurer 2002 + * 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: xor_combine.hpp 60755 2010-03-22 00:45:06Z steven_watanabe $ + * + */ + +#ifndef BOOST_RANDOM_XOR_COMBINE_HPP +#define BOOST_RANDOM_XOR_COMBINE_HPP + +#include +#include +#include // for std::min and std::max +#include +#include +#include +#include // uint32_t +#include + + +namespace boost { +namespace random { + +/// \cond hide_private_members +#ifndef BOOST_NO_DEPENDENT_TYPES_IN_TEMPLATE_VALUE_PARAMETERS + #define BOOST_RANDOM_VAL_TYPE typename URNG1::result_type +#else + #define BOOST_RANDOM_VAL_TYPE uint32_t +#endif +/// \endcond + +/** + * Instantiations of @c xor_combine model a \pseudo_random_number_generator. + * To produce its output it invokes each of the base generators, shifts + * their results and xors them together. + */ +template +class xor_combine +{ +public: + typedef URNG1 base1_type; + typedef URNG2 base2_type; + typedef typename base1_type::result_type result_type; + + BOOST_STATIC_CONSTANT(bool, has_fixed_range = false); + BOOST_STATIC_CONSTANT(int, shift1 = s1); + BOOST_STATIC_CONSTANT(int, shift2 = s2); + + /** + * Constructors a @c xor_combine by default constructing + * both base generators. + */ + xor_combine() : _rng1(), _rng2() + { } + /** + * Constructs a @c xor_combine by copying two base generators. + */ + xor_combine(const base1_type & rng1, const base2_type & rng2) + : _rng1(rng1), _rng2(rng2) { } + /** + * Constructs a @c xor_combine, seeding both base generators + * with @c v. + */ + xor_combine(const result_type & v) + : _rng1(v), _rng2(v) { } + /** + * Constructs a @c xor_combine, seeding both base generators + * with values from the iterator range [first, last) and changes + * first to point to the element after the last one used. If there + * are not enough elements in the range to seed both generators, + * throws @c std::invalid_argument. + */ + template xor_combine(It& first, It last) + : _rng1(first, last), _rng2( /* advanced by other call */ first, last) { } + /** + * Calls @c seed() for both base generators. + */ + void seed() { _rng1.seed(); _rng2.seed(); } + /** + * @c seeds both base generators with @c v. + */ + void seed(const result_type & v) { _rng1.seed(v); _rng2.seed(v); } + /** + * seeds both base generators with values from the iterator + * range [first, last) and changes first to point to the element + * after the last one used. If there are not enough elements in + * the range to seed both generators, throws @c std::invalid_argument. + */ + template void seed(It& first, It last) + { + _rng1.seed(first, last); + _rng2.seed(first, last); + } + + /** Returns the first base generator. */ + const base1_type& base1() { return _rng1; } + /** Returns the second base generator. */ + const base2_type& base2() { return _rng2; } + + /** Returns the next value of the generator. */ + result_type operator()() + { + // MSVC fails BOOST_STATIC_ASSERT with std::numeric_limits at class scope +#if !defined(BOOST_NO_LIMITS_COMPILE_TIME_CONSTANTS) && !(defined(BOOST_MSVC) && BOOST_MSVC <= 1300) + BOOST_STATIC_ASSERT(std::numeric_limits::is_integer); + BOOST_STATIC_ASSERT(std::numeric_limits::is_integer); + BOOST_STATIC_ASSERT(std::numeric_limits::digits >= std::numeric_limits::digits); +#endif + return (_rng1() << s1) ^ (_rng2() << s2); + } + + /** + * Returns the smallest value that the generator can produce. + */ + result_type min BOOST_PREVENT_MACRO_SUBSTITUTION () const { return std::min BOOST_PREVENT_MACRO_SUBSTITUTION((_rng1.min)(), (_rng2.min)()); } + /** + * Returns the largest value that the generator can produce. + */ + result_type max BOOST_PREVENT_MACRO_SUBSTITUTION () const { return std::max BOOST_PREVENT_MACRO_SUBSTITUTION((_rng1.min)(), (_rng2.max)()); } + static bool validation(result_type x) { return val == x; } + +#ifndef BOOST_NO_OPERATORS_IN_NAMESPACE + +#ifndef BOOST_RANDOM_NO_STREAM_OPERATORS + template + friend std::basic_ostream& + operator<<(std::basic_ostream& os, const xor_combine& s) + { + os << s._rng1 << " " << s._rng2 << " "; + return os; + } + + template + friend std::basic_istream& + operator>>(std::basic_istream& is, xor_combine& s) + { + is >> s._rng1 >> std::ws >> s._rng2 >> std::ws; + return is; + } +#endif + + friend bool operator==(const xor_combine& x, const xor_combine& y) + { return x._rng1 == y._rng1 && x._rng2 == y._rng2; } + friend bool operator!=(const xor_combine& x, const xor_combine& y) + { return !(x == y); } +#else + // Use a member function; Streamable concept not supported. + bool operator==(const xor_combine& rhs) const + { return _rng1 == rhs._rng1 && _rng2 == rhs._rng2; } + bool operator!=(const xor_combine& rhs) const + { return !(*this == rhs); } +#endif + +private: + base1_type _rng1; + base2_type _rng2; +}; + +#ifndef BOOST_NO_INCLASS_MEMBER_INITIALIZATION +// A definition is required even for integral static constants +template +const bool xor_combine::has_fixed_range; +template +const int xor_combine::shift1; +template +const int xor_combine::shift2; +#endif + +#undef BOOST_RANDOM_VAL_TYPE + +} // namespace random +} // namespace boost + +#endif // BOOST_RANDOM_XOR_COMBINE_HPP diff --git a/Boost/lib/libboost_system-vc100-mt-1_44.lib b/Boost/lib/libboost_system-vc100-mt-1_44.lib new file mode 100644 index 000000000..47510e6a5 Binary files /dev/null and b/Boost/lib/libboost_system-vc100-mt-1_44.lib differ diff --git a/Boost/lib/libboost_system-vc100-mt-gd-1_44.lib b/Boost/lib/libboost_system-vc100-mt-gd-1_44.lib new file mode 100644 index 000000000..49609c46b Binary files /dev/null and b/Boost/lib/libboost_system-vc100-mt-gd-1_44.lib differ