Subject: [Boost-commit] svn:boost r64580 - in sandbox/SOC/2010/quasi_random/boost/random: . detail
From: jvd_at_[hidden]
Date: 2010-08-03 15:04:32
Author: qrng
Date: 2010-08-03 15:04:31 EDT (Tue, 03 Aug 2010)
New Revision: 64580
URL: http://svn.boost.org/trac/boost/changeset/64580
Log:
Added heavily optimized integer pow and log functions. Synthetic test pass a little bit faster.
Added:
sandbox/SOC/2010/quasi_random/boost/random/detail/integer_powers.hpp (contents, props changed)
Text files modified:
sandbox/SOC/2010/quasi_random/boost/random/detail/gray_coded_qrng_base.hpp | 18 +++++++++-------
sandbox/SOC/2010/quasi_random/boost/random/faure.hpp | 43 ++-------------------------------------
2 files changed, 13 insertions(+), 48 deletions(-)
Modified: sandbox/SOC/2010/quasi_random/boost/random/detail/gray_coded_qrng_base.hpp
==============================================================================
--- sandbox/SOC/2010/quasi_random/boost/random/detail/gray_coded_qrng_base.hpp (original)
+++ sandbox/SOC/2010/quasi_random/boost/random/detail/gray_coded_qrng_base.hpp 2010-08-03 15:04:31 EDT (Tue, 03 Aug 2010)
@@ -10,6 +10,7 @@
#define BOOST_RANDOM_DETAIL_GRAY_CODED_QRNG_BASE_HPP
#include <boost/random/detail/qrng_base.hpp>
+#include <boost/random/detail/integer_powers.hpp>
//!\file
//!Describes the gray-coded quasi-random number generator base class template.
@@ -88,7 +89,7 @@
this->seq_count = init;
init ^= (init / 2);
- for( int r = 0; init != 0; ++r, init >>= 1 )
+ for( std::size_t r = 0; init != 0; ++r, init >>= 1 )
{
if( init & 1 )
update_quasi(r, msg);
@@ -97,17 +98,18 @@
}
private:
+ // Returns the position of the least significant 1 bit
+ inline static std::size_t lsb(std::size_t v)
+ {
+ return detail::integer_log<2>(v & (std::size_t(0)-v));
+ }
void compute_next(std::size_t seq)
{
// Find the position of the least-significant zero in sequence count.
// This is the bit that changes in the Gray-code representation as
// the count is advanced.
- int r = 0;
- for( ; seq & 1; seq >>= 1 ) {
- ++r;
- }
- update_quasi(r, "compute_next");
+ update_quasi(lsb(~seq), "compute_next");
}
// Initializes currently stored values to zero
@@ -116,9 +118,9 @@
std::fill(this->quasi_state, this->quasi_state + base_t::dimension, result_type());
}
- void update_quasi(int r, const char* msg)
+ void update_quasi(std::size_t r, const char* msg)
{
- if(r >= LatticeT::bit_count)
+ if(r >= (std::size_t)LatticeT::bit_count)
boost::throw_exception( std::overflow_error(msg) );
// Calculate the next state.
Added: sandbox/SOC/2010/quasi_random/boost/random/detail/integer_powers.hpp
==============================================================================
--- (empty file)
+++ sandbox/SOC/2010/quasi_random/boost/random/detail/integer_powers.hpp 2010-08-03 15:04:31 EDT (Tue, 03 Aug 2010)
@@ -0,0 +1,218 @@
+// Copyright Justinas Vygintas Daugmaudis, 2010.
+// Use, modification and distribution is subject to the
+// Boost Software License, Version 1.0. (See accompanying
+// file LICENSE-1.0 or http://www.boost.org/LICENSE-1.0)
+
+#ifndef BOOST_RANDOM_DETAIL_INTEGER_POWERS_HPP
+#define BOOST_RANDOM_DETAIL_INTEGER_POWERS_HPP
+
+#include <boost/static_assert.hpp>
+#include <boost/mpl/identity.hpp>
+#include <boost/mpl/integral_c.hpp>
+#include <boost/mpl/if.hpp>
+#include <boost/mpl/vector_c.hpp>
+#include <boost/mpl/back.hpp>
+#include <boost/mpl/push_back.hpp>
+#include <boost/mpl/at.hpp>
+#include <boost/mpl/size.hpp>
+#include <boost/preprocessor/repetition/repeat.hpp>
+#include <boost/preprocessor/iteration/local.hpp>
+
+namespace boost { namespace random {
+
+namespace detail {
+
+template<typename Powers, std::size_t Base, std::size_t MaxValue>
+struct generate_integer_powers_impl
+{
+ BOOST_STATIC_CONSTANT(std::size_t, last_power =
+ mpl::back<Powers>::type::value );
+
+ typedef typename mpl::push_back<
+ Powers
+ , mpl::integral_c<std::size_t, last_power * Base>
+ >::type powers_t;
+
+ // Computation must be lazy!
+ typedef typename mpl::if_c<
+ (last_power <= MaxValue)
+ , generate_integer_powers_impl<powers_t, Base, MaxValue>
+ , mpl::identity<Powers>
+ >::type selector_t;
+ typedef typename selector_t::type type;
+};
+
+template<std::size_t Base>
+struct generate_integer_powers
+{
+ BOOST_STATIC_ASSERT( Base >= 2 );
+
+ BOOST_STATIC_CONSTANT(std::size_t, max_bin_value = ~std::size_t(0) );
+ BOOST_STATIC_CONSTANT(std::size_t, max_value = max_bin_value / Base );
+
+ typedef boost::mpl::vector_c<std::size_t, 1> base0_power_t;
+ typedef generate_integer_powers_impl<base0_power_t, Base, max_value> gen_powers_t;
+
+ // Here we generate an MPL vector of values
+ // Base^0, ..., Base^m where Base^m <= max_value
+ typedef typename gen_powers_t::type type;
+};
+
+
+// Performs binary search and returns the lower bound of the range
+template<std::size_t Pos>
+struct range_not_found
+{
+ inline static std::size_t apply(std::size_t)
+ { return Pos; }
+};
+
+template<std::size_t Pos>
+struct range_prev
+{
+ BOOST_STATIC_CONSTANT(std::size_t, value = Pos - 1);
+};
+
+template<>
+struct range_prev<0>
+{
+ BOOST_STATIC_CONSTANT(std::size_t, value = 0); // integer_log<Base>(0) == 0!
+};
+
+template<typename Powers, std::size_t Low, std::size_t High>
+struct log_lower_bound
+{
+ BOOST_STATIC_CONSTANT(std::size_t, mid = (Low + High) / 2);
+
+ // Termination conditions are described here
+ typedef log_lower_bound<Powers, Low, mid - 1> lower_search_t;
+ typedef typename mpl::if_c< (Low >= mid)
+ , range_not_found<range_prev<mid>::value>
+ , lower_search_t
+ >::type lower_t;
+ typedef log_lower_bound<Powers, mid + 1, High> upper_search_t;
+ typedef typename mpl::if_c< (mid >= High)
+ , range_not_found<High>
+ , upper_search_t
+ >::type upper_t;
+
+ typedef typename mpl::at_c<Powers, mid>::type value_t;
+
+ inline static std::size_t apply(std::size_t arg)
+ {
+ if( arg > value_t::value )
+ return upper_t::apply(arg);
+ else if( arg < value_t::value )
+ return lower_t::apply(arg);
+ else
+ return mid;
+ }
+};
+
+// 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.
+// The function performs binary range search, so for big args
+// it works substantially faster (measured ~5x) than the trivial
+// std::size_t ilog = 0;
+// while( Base <= arg )
+// {
+// arg /= Base; ++ilog;
+// }
+template<std::size_t Base>
+inline std::size_t integer_log(std::size_t arg)
+{
+ typedef typename generate_integer_powers<Base>::type powers_t;
+ typedef typename mpl::size<powers_t>::type size_t; // always >= 1
+ typedef log_lower_bound<powers_t, 0, size_t::value - 1> binlookup_t;
+ return binlookup_t::apply(arg);
+}
+
+
+// This is more efficient than naively multiplying the base with itself repeatedly.
+// In erroneous situations, e.g., integer_pow<0>(0) the function returns 0
+// and does not report the error. This is the intended behavior.
+
+namespace pow_switch_detail {
+
+// D is the number of cases without the default
+template<std::size_t D>
+struct pow_switch_impl;
+
+// specialize for 0 to appease the compiler,
+// but we know there cannot be size 0 vectors
+template<>
+struct pow_switch_impl<0> {
+ template<class V>
+ inline static std::size_t apply(std::size_t) {
+ return 0;
+ }
+};
+
+#define BOOST_RANDOM_DETAIL_INTEGER_POW_CASE_IMPL(z, N, data) \
+ case N: { \
+ return boost::mpl::at_c<data, N>::type::value; \
+ } \
+/**/
+
+// Here we assume that compilers nowadays are clever enough to optimize this monstrosity
+#define BOOST_RANDOM_DETAIL_INTEGER_POW_SWITCH_IMPL(z, N, data) \
+ template<> \
+ struct pow_switch_impl<N> { \
+ template<typename Seq> \
+ inline static std::size_t apply(std::size_t arg) { \
+ switch( arg ) { \
+ BOOST_PP_REPEAT_##z(N, BOOST_RANDOM_DETAIL_INTEGER_POW_CASE_IMPL, Seq) \
+ default: return 0; \
+ } \
+ } \
+ }; \
+/**/
+
+#define BOOST_PP_LOCAL_LIMITS (1, 64) // up to 2^63
+#define BOOST_PP_LOCAL_MACRO(n) BOOST_RANDOM_DETAIL_INTEGER_POW_SWITCH_IMPL(1, n, ~)
+#include BOOST_PP_LOCAL_ITERATE()
+
+#undef BOOST_RANDOM_DETAIL_INTEGER_POW_SWITCH_IMPL
+#undef BOOST_RANDOM_DETAIL_INTEGER_POW_CASE_IMPL
+
+} // namespace pow_switch_detail
+
+template<std::size_t Base>
+struct integer_pow_func
+{
+ typedef typename generate_integer_powers<Base>::type powers_t;
+ typedef pow_switch_detail::pow_switch_impl<boost::mpl::size<powers_t>::value> switch_t;
+
+ inline static std::size_t apply(std::size_t arg)
+ {
+ return switch_t::template apply<powers_t>(arg);
+ }
+};
+
+template<>
+struct integer_pow_func<1> // 1^p == 1
+{
+ inline static std::size_t apply(std::size_t)
+ { return 1; }
+};
+
+template<>
+struct integer_pow_func<0> // 0^p == 0, also 0^0 == 0!
+{
+ inline static std::size_t apply(std::size_t)
+ { return 0; }
+};
+
+template<std::size_t Base>
+inline std::size_t integer_pow(std::size_t arg)
+{
+ return integer_pow_func<Base>::apply(arg);
+}
+
+} // namespace detail
+
+}} // namespace boost::random
+
+#endif // BOOST_RANDOM_DETAIL_INTEGER_POWERS_HPP
Modified: sandbox/SOC/2010/quasi_random/boost/random/faure.hpp
==============================================================================
--- sandbox/SOC/2010/quasi_random/boost/random/faure.hpp (original)
+++ sandbox/SOC/2010/quasi_random/boost/random/faure.hpp 2010-08-03 15:04:31 EDT (Tue, 03 Aug 2010)
@@ -12,6 +12,7 @@
#include <boost/random/detail/qrng_base.hpp>
#include <boost/random/detail/config.hpp>
#include <boost/random/detail/operators.hpp>
+#include <boost/random/detail/integer_powers.hpp>
#include <vector>
#include <cmath>
@@ -55,44 +56,6 @@
BOOST_STATIC_CONSTANT(int, value = mpl::deref<iter_t>::type::value);
};
-// 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<std::size_t Base>
-inline std::size_t integer_log(std::size_t arg)
-{
- BOOST_STATIC_ASSERT( 2 <= Base );
-
- std::size_t ilog = 0;
- while( Base <= arg )
- {
- arg /= Base; ++ilog;
- }
- return ilog;
-}
-
-// Implements unrolled exponentiation by squaring. For p > ~4 this is computationally
-// more efficient than naively multiplying the base with itself repeatedly.
-// In erroneous situations, e.g., integer_pow(0, 0) the function returns 1
-// and does not report the error. This is the intended behavior.
-inline std::size_t mdelta(std::size_t base, std::size_t p)
-{
- return (p & 1) * base + !(p & 1); // (p & 1) ? base : 1
-}
-
-inline std::size_t integer_pow(std::size_t base, std::size_t p)
-{
- std::size_t result = 1;
- for( ; p != 0; p >>= 1, base *= base ) {
- result *= mdelta(base, p);
- result *= mdelta(base *= base, p >>= 1);
- result *= mdelta(base *= base, p >>= 1);
- result *= mdelta(base *= base, p >>= 1);
- }
- return result;
-}
-
// Computes a table of binomial coefficients modulo qs.
template<typename RealType, std::size_t Dimension>
struct binomial_coefficients
@@ -160,7 +123,7 @@
private:
inline static std::size_t n_elements(std::size_t seq)
{
- return integer_log<qs_base>(seq) + 1;
+ return detail::integer_log<qs_base>(seq) + 1;
}
inline static std::size_t size_hint(std::size_t n)
@@ -184,7 +147,7 @@
// Sum ( 0 <= J <= HISUM ) YTEMP(J) / QS**(J+1)
// in one go
RealType r = RealType();
- std::size_t m, k = integer_pow(qs_base, n - 1);
+ std::size_t m, k = detail::integer_pow<qs_base>(n - 1);
for( ; n != 0; --n, ++out, seq = m, k /= qs_base )
{
m = seq % k;