Boost logo

Boost-Commit :

Subject: [Boost-commit] svn:boost r70047 - in trunk/boost/random: . detail
From: steven_at_[hidden]
Date: 2011-03-16 23:45:40


Author: steven_watanabe
Date: 2011-03-16 23:45:39 EDT (Wed, 16 Mar 2011)
New Revision: 70047
URL: http://svn.boost.org/trac/boost/changeset/70047

Log:
Replace a couple slow iterative algorithms with multi-precision arithmetic.
Added:
   trunk/boost/random/detail/large_arithmetic.hpp (contents, props changed)
Text files modified:
   trunk/boost/random/detail/const_mod.hpp | 54 +++++++++------------------------------
   trunk/boost/random/shuffle_order.hpp | 35 +++++--------------------
   2 files changed, 20 insertions(+), 69 deletions(-)

Modified: trunk/boost/random/detail/const_mod.hpp
==============================================================================
--- trunk/boost/random/detail/const_mod.hpp (original)
+++ trunk/boost/random/detail/const_mod.hpp 2011-03-16 23:45:39 EDT (Wed, 16 Mar 2011)
@@ -16,25 +16,17 @@
 #ifndef BOOST_RANDOM_CONST_MOD_HPP
 #define BOOST_RANDOM_CONST_MOD_HPP
 
-#include <algorithm>
 #include <boost/assert.hpp>
 #include <boost/static_assert.hpp>
-#include <boost/cstdint.hpp>
 #include <boost/integer_traits.hpp>
 #include <boost/type_traits/make_unsigned.hpp>
-#include <boost/detail/workaround.hpp>
+#include <boost/random/detail/large_arithmetic.hpp>
 
 #include <boost/random/detail/disable_warnings.hpp>
 
 namespace boost {
 namespace random {
 
-/*
- * Some random number generators require modular arithmetic. Put
- * everything we need here.
- * IntType must be an integral type.
- */
-
 template<class IntType, IntType m>
 class const_mod
 {
@@ -74,8 +66,6 @@
       return mult_small(a, x);
     else if(traits::is_signed && (m%a < m/a))
       return mult_schrage(a, x);
- else if(x == 1)
- return a;
     else
       return mult_general(a, x);
   }
@@ -133,36 +123,18 @@
     return sub(a*(value%q), r*(value/q));
   }
 
- static IntType mult_general(IntType a, IntType b) {
- IntType q, r;
- IntType value = 0;
-
- IntType supress_warnings = (m == 0);
- BOOST_ASSERT(supress_warnings == 0);
-
- while(true) {
- if(a == 0 || b <= traits::const_max/a)
- return add(value, IntType(a * b % (m + supress_warnings)));
-
- if(b < a) std::swap(a, b);
- q = m / a;
- r = m % a;
-
- value = add(value, IntType(a*(b%q)));
-
- a = r;
- b = IntType(b/q);
- if(a == 0 || b <= traits::const_max/a)
- return sub(value, IntType(a * b % (m + supress_warnings)));
-
- if(b < a) std::swap(a, b);
- q = m / a;
- r = m % a;
-
- value = sub(value, IntType(a*(b%q)));
-
- a = r;
- b = IntType(b/q);
+ static IntType mult_general(IntType a, IntType b)
+ {
+ IntType suppress_warnings = (m == 0);
+ BOOST_ASSERT(suppress_warnings == 0);
+ IntType modulus = m + suppress_warnings;
+ BOOST_ASSERT(modulus == m);
+ if(::boost::uintmax_t(modulus) <=
+ (::std::numeric_limits< ::boost::uintmax_t>::max)() / modulus)
+ {
+ return static_cast<IntType>(boost::uintmax_t(a) * b % modulus);
+ } else {
+ return static_cast<IntType>(detail::mulmod(a, b, modulus));
     }
   }
 

Added: trunk/boost/random/detail/large_arithmetic.hpp
==============================================================================
--- (empty file)
+++ trunk/boost/random/detail/large_arithmetic.hpp 2011-03-16 23:45:39 EDT (Wed, 16 Mar 2011)
@@ -0,0 +1,122 @@
+/* boost random/detail/large_arithmetic.hpp header file
+ *
+ * Copyright Steven Watanabe 2011
+ * Distributed under the Boost Software License, Version 1.0. (See
+ * accompanying file LICENSE_1_0.txt or copy at
+ * http://www.boost.org/LICENSE_1_0.txt)
+ *
+ * See http://www.boost.org for most recent version including documentation.
+ *
+ * $Id$
+ */
+
+#ifndef BOOST_RANDOM_DETAIL_LARGE_ARITHMETIC_HPP
+#define BOOST_RANDOM_DETAIL_LARGE_ARITHMETIC_HPP
+
+#include <boost/cstdint.hpp>
+#include <boost/integer.hpp>
+#include <boost/limits.hpp>
+#include <boost/random/detail/integer_log2.hpp>
+
+#include <boost/random/detail/disable_warnings.hpp>
+
+namespace boost {
+namespace random {
+namespace detail {
+
+struct div_t {
+ boost::uintmax_t quotient;
+ boost::uintmax_t remainder;
+};
+
+div_t muldivmod(boost::uintmax_t a, boost::uintmax_t b, boost::uintmax_t m)
+{
+ static const int bits =
+ ::std::numeric_limits< ::boost::uintmax_t>::digits / 2;
+ static const ::boost::uintmax_t mask = (::boost::uintmax_t(1) << bits) - 1;
+ typedef ::boost::uint_t<bits>::fast digit_t;
+
+ int shift = std::numeric_limits< ::boost::uintmax_t>::digits - 1
+ - detail::integer_log2(m);
+
+ a <<= shift;
+ m <<= shift;
+
+ digit_t product[4] = { 0, 0, 0, 0 };
+ digit_t a_[2] = { digit_t(a & mask), digit_t((a >> bits) & mask) };
+ digit_t b_[2] = { digit_t(b & mask), digit_t((b >> bits) & mask) };
+ digit_t m_[2] = { digit_t(m & mask), digit_t((m >> bits) & mask) };
+
+ // multiply a * b
+ for(int i = 0; i < 2; ++i) {
+ digit_t carry = 0;
+ for(int j = 0; j < 2; ++j) {
+ ::boost::uint64_t temp = ::boost::uintmax_t(a_[i]) * b_[j] +
+ carry + product[i + j];
+ product[i + j] = digit_t(temp & mask);
+ carry = digit_t(temp >> bits);
+ }
+ if(carry != 0) {
+ product[i + 2] += carry;
+ }
+ }
+
+ digit_t quotient[2];
+
+ if(m == 0) {
+ div_t result = {
+ ((::boost::uintmax_t(product[3]) << bits) | product[2]),
+ ((::boost::uintmax_t(product[1]) << bits) | product[0]) >> shift,
+ };
+ return result;
+ }
+
+ // divide product / m
+ for(int i = 3; i >= 2; --i) {
+ ::boost::uintmax_t temp =
+ ::boost::uintmax_t(product[i]) << bits | product[i - 1];
+
+ digit_t q = digit_t((product[i] == m_[1]) ? mask : temp / m_[1]);
+
+ ::boost::uintmax_t rem =
+ ((temp - ::boost::uintmax_t(q) * m_[1]) << bits) + product[i - 2];
+
+ ::boost::uintmax_t diff = m_[0] * ::boost::uintmax_t(q);
+
+ int error = 0;
+ if(diff > rem) {
+ if(diff - rem > m) {
+ error = 2;
+ } else {
+ error = 1;
+ }
+ }
+ q -= error;
+ rem = rem + error * m - diff;
+
+ quotient[i - 2] = q;
+ product[i] = 0;
+ product[i-1] = (rem >> bits) & mask;
+ product[i-2] = rem & mask;
+ }
+
+ div_t result = {
+ ((::boost::uintmax_t(quotient[1]) << bits) | quotient[0]),
+ ((::boost::uintmax_t(product[1]) << bits) | product[0]) >> shift,
+ };
+ return result;
+}
+
+boost::uintmax_t muldiv(boost::uintmax_t a, boost::uintmax_t b, boost::uintmax_t m)
+{ return detail::muldivmod(a, b, m).quotient; }
+
+boost::uintmax_t mulmod(boost::uintmax_t a, boost::uintmax_t b, boost::uintmax_t m)
+{ return detail::muldivmod(a, b, m).remainder; }
+
+} // namespace detail
+} // namespace random
+} // namespace boost
+
+#include <boost/random/detail/enable_warnings.hpp>
+
+#endif // BOOST_RANDOM_DETAIL_LARGE_ARITHMETIC_HPP

Modified: trunk/boost/random/shuffle_order.hpp
==============================================================================
--- trunk/boost/random/shuffle_order.hpp (original)
+++ trunk/boost/random/shuffle_order.hpp 2011-03-16 23:45:39 EDT (Wed, 16 Mar 2011)
@@ -143,36 +143,15 @@
             // try to do it in the native type if we know that it won't
             // overflow
             j = k * off / (brange + 1);
- }
-#if !defined(BOOST_NO_INT64_T)
- else if(brange < (std::numeric_limits<uint64_t>::max)() / k) {
+ } else if(brange < (std::numeric_limits<uintmax_t>::max)() / k) {
             // Otherwise try to use uint64_t
             j = static_cast<base_unsigned>(
- static_cast<uint64_t>(off) * k /
- (static_cast<uint64_t>(brange) + 1));
- }
-#endif
- else {
- // If all else fails, fall back to a general algorithm that
- // never overflows.
-
- const base_unsigned r_mod_k = ((brange % k) + 1) % k;
- const base_unsigned bucket_size = (brange - k + 1)/k + 1;
- // if the candidate from the first round is zero, we're safe.
- base_unsigned candidate = 0;
- base_unsigned old_candidate;
- base_unsigned error = 0;
- do {
- old_candidate = candidate;
- candidate = (off - error) / bucket_size;
- base_unsigned possible = (off - error + 1) / bucket_size;
- error = possible - possible * (k - r_mod_k) / k;
- } while(old_candidate != candidate);
-
- j = candidate;
-
- // Would cause overflow
- // assert(j == uint64_t(off)*k/(uint64_t(brange)+1));
+ static_cast<uintmax_t>(off) * k /
+ (static_cast<uintmax_t>(brange) + 1));
+ } else {
+ boost::uintmax_t divisor =
+ static_cast<boost::uintmax_t>(brange) + 1;
+ j = static_cast<base_unsigned>(detail::muldiv(off, k, divisor));
         }
         // assert(0 <= j && j < k);
         y = v[j];


Boost-Commit list run by bdawes at acm.org, david.abrahams at rcn.com, gregod at cs.rpi.edu, cpdaniel at pacbell.net, john at johnmaddock.co.uk