diff options
Diffstat (limited to 'src/third_party/boost-1.69.0/boost/random/faure.hpp')
| -rw-r--r-- | src/third_party/boost-1.69.0/boost/random/faure.hpp | 367 |
1 files changed, 367 insertions, 0 deletions
diff --git a/src/third_party/boost-1.69.0/boost/random/faure.hpp b/src/third_party/boost-1.69.0/boost/random/faure.hpp new file mode 100644 index 00000000000..4301c301b5f --- /dev/null +++ b/src/third_party/boost-1.69.0/boost/random/faure.hpp @@ -0,0 +1,367 @@ +/* boost random/faure.hpp header file + * + * Copyright Justinas Vygintas Daugmaudis 2010-2018 + * 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) + */ + +#ifndef BOOST_RANDOM_FAURE_HPP +#define BOOST_RANDOM_FAURE_HPP + +#include <boost/random/detail/qrng_base.hpp> + +#include <cmath> +#include <vector> +#include <algorithm> + +#include <boost/assert.hpp> + +namespace boost { +namespace random { + +/** @cond */ +namespace detail { + +namespace qrng_tables { + +// There is no particular reason why 187 first primes were chosen +// to be put into this table. The only reason was, perhaps, that +// the number of dimensions for Faure generator would be around +// the same order of magnitude as the number of dimensions supported +// by the Sobol qrng. +struct primes +{ + typedef unsigned short value_type; + + BOOST_STATIC_CONSTANT(int, number_of_primes = 187); + + // A function that returns lower bound prime for a given n + static value_type lower_bound(std::size_t n) + { + static const value_type prim_a[number_of_primes] = { + 2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, + 59, 61, 67, 71, 73, 79, 83, 89, 97, 101, 103, 107, 109, 113, + 127, 131, 137, 139, 149, 151, 157, 163, 167, 173, 179, 181, + 191, 193, 197, 199, 211, 223, 227, 229, 233, 239, 241, 251, + 257, 263, 269, 271, 277, 281, 283, 293, 307, 311, 313, 317, + 331, 337, 347, 349, 353, 359, 367, 373, 379, 383, 389, 397, + 401, 409, 419, 421, 431, 433, 439, 443, 449, 457, 461, 463, + 467, 479, 487, 491, 499, 503, 509, 521, 523, 541, 547, 557, + 563, 569, 571, 577, 587, 593, 599, 601, 607, 613, 617, 619, + 631, 641, 643, 647, 653, 659, 661, 673, 677, 683, 691, 701, + 709, 719, 727, 733, 739, 743, 751, 757, 761, 769, 773, 787, + 797, 809, 811, 821, 823, 827, 829, 839, 853, 857, 859, 863, + 877, 881, 883, 887, 907, 911, 919, 929, 937, 941, 947, 953, + 967, 971, 977, 983, 991, 997, 1009, 1013, 1019, 1021, 1031, + 1033, 1039, 1049, 1051, 1061, 1063, 1069, 1087, 1091, 1093, + 1097, 1103, 1109, 1117 }; + + qrng_detail::dimension_assert("Faure", n, prim_a[number_of_primes - 1]); + + return *std::lower_bound(prim_a, prim_a + number_of_primes, n); + } +}; + +} // namespace qrng_tables +} // namespace detail + +namespace qrng_detail { +namespace fr { + +// Returns the integer part of the logarithm base Base of arg. +// In erroneous situations, e.g., integer_log(base, 0) the function +// returns 0 and does not report the error. This is the intended +// behavior. +template <typename T> +inline T integer_log(T base, T arg) +{ + T ilog = T(); + while (base <= arg) + { + arg /= base; ++ilog; + } + return ilog; +} + +// Perform exponentiation by squaring (potential for code reuse in multiprecision::powm) +template <typename T> +inline T integer_pow(T base, T e) +{ + T result = static_cast<T>(1); + while (e) + { + if (e & static_cast<T>(1)) + result *= base; + e >>= 1; + base *= base; + } + return result; +} + +} // namespace fr + +// Computes a table of binomial coefficients modulo qs. +template<typename RealType, typename SeqSizeT, typename PrimeTable> +struct binomial_coefficients +{ + typedef RealType value_type; + typedef SeqSizeT size_type; + + // Binomial values modulo qs_base will never be bigger than qs_base. + // We can choose an appropriate integer type to hold modulo values and + // shave off memory footprint. + typedef typename PrimeTable::value_type packed_uint_t; + + // default copy c-tor is fine + + explicit binomial_coefficients(std::size_t dimension) + { + resize(dimension); + } + + void resize(std::size_t dimension) + { + qs_base = PrimeTable::lower_bound(dimension); + + // Throw away previously computed coefficients. + // This will trigger recomputation on next update + coeff.clear(); + } + + template <typename Iterator> + void update(size_type seq, Iterator first, Iterator last) + { + if (first != last) + { + const size_type ilog = fr::integer_log(static_cast<size_type>(qs_base), seq); + const size_type hisum = ilog + 1; + if (coeff.size() != size_hint(hisum)) { + ytemp.resize(static_cast<std::size_t>(hisum)); // cast safe because log is small + compute_coefficients(hisum); + qs_pow = fr::integer_pow(static_cast<size_type>(qs_base), ilog); + } + + *first = compute_recip(seq, ytemp.rbegin()); + + // Find other components using the Faure method. + ++first; + for ( ; first != last; ++first) + { + *first = RealType(); + RealType r = static_cast<RealType>(1); + + for (size_type i = 0; i != hisum; ++i) + { + RealType ztemp = ytemp[static_cast<std::size_t>(i)] * upper_element(i, i, hisum); + for (size_type j = i + 1; j != hisum; ++j) + ztemp += ytemp[static_cast<std::size_t>(j)] * upper_element(i, j, hisum); + + // Sum ( J <= I <= HISUM ) ( old ytemp(i) * binom(i,j) ) mod QS. + ytemp[static_cast<std::size_t>(i)] = std::fmod(ztemp, static_cast<RealType>(qs_base)); + r *= static_cast<RealType>(qs_base); + *first += ytemp[static_cast<std::size_t>(i)] / r; + } + } + } + } + +private: + inline static size_type size_hint(size_type n) + { + return n * (n + 1) / 2; + } + + packed_uint_t& upper_element(size_type i, size_type j, size_type dim) + { + BOOST_ASSERT( i < dim ); + BOOST_ASSERT( j < dim ); + BOOST_ASSERT( i <= j ); + return coeff[static_cast<std::size_t>((i * (2 * dim - i + 1)) / 2 + j - i)]; + } + + template<typename Iterator> + RealType compute_recip(size_type seq, Iterator out) const + { + // Here we do + // Sum ( 0 <= J <= HISUM ) YTEMP(J) * QS**J + // Sum ( 0 <= J <= HISUM ) YTEMP(J) / QS**(J+1) + // in one go + RealType r = RealType(); + size_type m, k = qs_pow; + for( ; k != 0; ++out, seq = m, k /= qs_base ) + { + m = seq % k; + RealType v = static_cast<RealType>((seq - m) / k); // RealType <- size type + r += v; + r /= static_cast<RealType>(qs_base); + *out = v; // saves double dereference + } + return r; + } + + void compute_coefficients(const size_type n) + { + // Resize and initialize to zero + coeff.resize(static_cast<std::size_t>(size_hint(n))); + std::fill(coeff.begin(), coeff.end(), packed_uint_t()); + + // The first row and the diagonal is assigned to 1 + upper_element(0, 0, n) = 1; + for (size_type i = 1; i < n; ++i) + { + upper_element(0, i, n) = 1; + upper_element(i, i, n) = 1; + } + + // Computes binomial coefficients MOD qs_base + for (size_type i = 1; i < n; ++i) + { + for (size_type j = i + 1; j < n; ++j) + { + upper_element(i, j, n) = ( upper_element(i, j-1, n) + + upper_element(i-1, j-1, n) ) % qs_base; + } + } + } + +private: + packed_uint_t qs_base; + + // here we cache precomputed data; note that binomial coefficients have + // to be recomputed iff the integer part of the logarithm of seq changes, + // which happens relatively rarely. + std::vector<packed_uint_t> coeff; // packed upper (!) triangular matrix + std::vector<RealType> ytemp; + size_type qs_pow; +}; + +} // namespace qrng_detail + +typedef detail::qrng_tables::primes default_faure_prime_table; + +/** @endcond */ + +//!Instantiations of class template faure_engine model a \quasi_random_number_generator. +//!The faure_engine uses the algorithm described in +//! \blockquote +//!Henri Faure, +//!Discrepance de suites associees a un systeme de numeration (en dimension s), +//!Acta Arithmetica, +//!Volume 41, 1982, pages 337-351. +//! \endblockquote +// +//! \blockquote +//!Bennett Fox, +//!Algorithm 647: +//!Implementation and Relative Efficiency of Quasirandom +//!Sequence Generators, +//!ACM Transactions on Mathematical Software, +//!Volume 12, Number 4, December 1986, pages 362-376. +//! \endblockquote +//! +//!In the following documentation @c X denotes the concrete class of the template +//!faure_engine returning objects of type @c RealType, u and v are the values of @c X. +//! +//!Some member functions may throw exceptions of type @c std::bad_alloc. +template<typename RealType, typename SeqSizeT, typename PrimeTable = default_faure_prime_table> +class faure_engine + : public qrng_detail::qrng_base< + faure_engine<RealType, SeqSizeT, PrimeTable> + , qrng_detail::binomial_coefficients<RealType, SeqSizeT, PrimeTable> + , SeqSizeT + > +{ + typedef faure_engine<RealType, SeqSizeT, PrimeTable> self_t; + + typedef qrng_detail::binomial_coefficients<RealType, SeqSizeT, PrimeTable> lattice_t; + typedef qrng_detail::qrng_base<self_t, lattice_t, SeqSizeT> base_t; + + friend class qrng_detail::qrng_base<self_t, lattice_t, SeqSizeT>; + +public: + typedef RealType result_type; + + /** @copydoc boost::random::niederreiter_base2_engine::min() */ + static BOOST_CONSTEXPR result_type min BOOST_PREVENT_MACRO_SUBSTITUTION () + { return static_cast<result_type>(0); } + + /** @copydoc boost::random::niederreiter_base2_engine::max() */ + static BOOST_CONSTEXPR result_type max BOOST_PREVENT_MACRO_SUBSTITUTION () + { return static_cast<result_type>(1); } + + //!Effects: Constructs the `s`-dimensional default Faure quasi-random number generator. + //! + //!Throws: bad_alloc, invalid_argument. + explicit faure_engine(std::size_t s) + : base_t(s) // initialize the binomial table here + {} + + /** @copydetails boost::random::niederreiter_base2_engine::seed(UIntType) + * Throws: bad_alloc. + */ + void seed(SeqSizeT init = 0) + { + compute_seq(init); + base_t::reset_seq(init); + } + +#ifdef BOOST_RANDOM_DOXYGEN + //=========================Doxygen needs this!============================== + + /** @copydoc boost::random::niederreiter_base2_engine::dimension() */ + std::size_t dimension() const { return base_t::dimension(); } + + /** @copydoc boost::random::niederreiter_base2_engine::operator()() */ + result_type operator()() + { + return base_t::operator()(); + } + + /** @copydoc boost::random::niederreiter_base2_engine::discard(boost::uintmax_t) + * Throws: bad_alloc. + */ + void discard(boost::uintmax_t z) + { + base_t::discard(z); + } + + /** Returns true if the two generators will produce identical sequences of outputs. */ + BOOST_RANDOM_DETAIL_EQUALITY_OPERATOR(faure_engine, x, y) + { return static_cast<const base_t&>(x) == y; } + + /** Returns true if the two generators will produce different sequences of outputs. */ + BOOST_RANDOM_DETAIL_INEQUALITY_OPERATOR(faure_engine) + + /** Writes the textual representation of the generator to a @c std::ostream. */ + BOOST_RANDOM_DETAIL_OSTREAM_OPERATOR(os, faure_engine, s) + { return os << static_cast<const base_t&>(s); } + + /** Reads the textual representation of the generator from a @c std::istream. */ + BOOST_RANDOM_DETAIL_ISTREAM_OPERATOR(is, faure_engine, s) + { return is >> static_cast<base_t&>(s); } + +#endif // BOOST_RANDOM_DOXYGEN + +private: +/** @cond hide_private_members */ + void compute_seq(SeqSizeT seq) + { + qrng_detail::check_seed_sign(seq); + this->lattice.update(seq, this->state_begin(), this->state_end()); + } +/** @endcond */ +}; + +/** + * @attention This specialization of \faure_engine supports up to 1117 dimensions. + * + * However, it is possible to provide your own prime table to \faure_engine should the default one be insufficient. + */ +typedef faure_engine<double, boost::uint_least64_t, default_faure_prime_table> faure; + +} // namespace random + +} // namespace boost + +#endif // BOOST_RANDOM_FAURE_HPP |