[Boost-commit] svn:boost r64580 - in sandbox/SOC/2010/quasi_random/boost/random: . detail

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;