*
* See http://www.boost.org for most recent version including documentation.
*
- * $Id: const_mod.hpp 58649 2010-01-02 21:23:17Z steven_watanabe $
+ * $Id: const_mod.hpp 71018 2011-04-05 21:27:52Z 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 <cassert>
+#include <boost/assert.hpp>
#include <boost/static_assert.hpp>
-#include <boost/cstdint.hpp>
#include <boost/integer_traits.hpp>
-#include <boost/detail/workaround.hpp>
+#include <boost/type_traits/make_unsigned.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.
- */
-
-namespace detail {
-
- template<bool is_signed>
- struct do_add
- { };
-
- template<>
- struct do_add<true>
- {
- template<class IntType>
- 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<false>
- {
- template<class IntType>
- 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 IntType, IntType m>
class const_mod
{
public:
+ static IntType apply(IntType x)
+ {
+ if(((unsigned_m() - 1) & unsigned_m()) == 0)
+ return (unsigned_type(x)) & (unsigned_m() - 1);
+ else {
+ IntType supress_warnings = (m == 0);
+ BOOST_ASSERT(supress_warnings == 0);
+ return x % (m + supress_warnings);
+ }
+ }
+
static IntType add(IntType x, IntType c)
{
- if(c == 0)
+ if(((unsigned_m() - 1) & unsigned_m()) == 0)
+ return (unsigned_type(x) + unsigned_type(c)) & (unsigned_m() - 1);
+ else if(c == 0)
return x;
- else if(c <= traits::const_max - m) // i.e. m+c < max
- return add_small(x, c);
+ else if(x < m - c)
+ return x + c;
else
- return detail::do_add<traits::is_signed>::add(m, x, c);
+ return x - (m - c);
}
static IntType mult(IntType a, IntType x)
{
- if(a == 1)
+ if(((unsigned_m() - 1) & unsigned_m()) == 0)
+ return unsigned_type(a) * unsigned_type(x) & (unsigned_m() - 1);
+ else if(a == 0)
+ return 0;
+ 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;
- }
+ else
+ return mult_general(a, x);
}
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
+ if(((unsigned_m() - 1) & unsigned_m()) == 0)
+ return (unsigned_type(a) * unsigned_type(x) + unsigned_type(c)) & (unsigned_m() - 1);
+ else if(a == 0)
+ return c;
+ else if(m <= (traits::const_max-c)/a) { // i.e. a*m+c <= max
+ IntType supress_warnings = (m == 0);
+ BOOST_ASSERT(supress_warnings == 0);
+ return (a*x+c) % (m + supress_warnings);
+ } else
return add(mult(a, x), c);
}
+ static IntType pow(IntType a, boost::uintmax_t exponent)
+ {
+ IntType result = 1;
+ while(exponent != 0) {
+ if(exponent % 2 == 1) {
+ result = mult(result, a);
+ }
+ a = mult(a, a);
+ exponent /= 2;
+ }
+ return result;
+ }
+
static IntType invert(IntType x)
- { return x == 0 ? 0 : invert_euclidian(x); }
+ { return x == 0 ? 0 : (m == 0? invert_euclidian0(x) : invert_euclidian(x)); }
private:
typedef integer_traits<IntType> traits;
+ typedef typename make_unsigned<IntType>::type unsigned_type;
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;
+ IntType supress_warnings = (m == 0);
+ BOOST_ASSERT(supress_warnings == 0);
+ return a*x % (m + supress_warnings);
}
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
+ BOOST_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 sub(a*(value%q), r*(value/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));
}
- return value;
+ }
+
+ static IntType sub(IntType a, IntType b)
+ {
+ if(a < b)
+ return m - (b - a);
+ else
+ return a - b;
+ }
+
+ static unsigned_type unsigned_m()
+ {
+ if(m == 0) {
+ return unsigned_type((std::numeric_limits<IntType>::max)()) + 1;
+ } else {
+ return unsigned_type(m);
+ }
}
// 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<IntType>::is_signed);
-#elif !defined(BOOST_NO_LIMITS_COMPILE_TIME_CONSTANTS)
- BOOST_STATIC_ASSERT(boost::integer_traits<IntType>::is_signed);
-#endif
- assert(c > 0);
+ BOOST_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!
+ l1 += q * l2;
p -= q * n;
if(p == 0)
- return (l2 < 1 ? l2 + m : l2);
+ return l2;
IntType q2 = n / p;
- l2 -= q2 * l1;
+ l2 += q2 * l1;
n -= q2 * p;
if(n == 0)
- return (l1 < 1 ? l1 + m : l1);
+ return 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<unsigned int, 0>
-{
- 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<unsigned long, 0>
-{
- 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<uint64_t, uint64_t(1) << 48>
-{
- 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 IntType, IntType m>
-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<traits::is_signed>::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<IntType> 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)
+ // invert c in the finite field (mod m) (c must be relatively prime to m)
+ static IntType invert_euclidian0(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<IntType>::is_signed);
-#endif
- assert(c > 0);
+ BOOST_ASSERT(c > 0);
+ if(c == 1) return 1;
IntType l1 = 0;
IntType l2 = 1;
IntType n = c;
IntType p = m;
+ IntType max = (std::numeric_limits<IntType>::max)();
+ IntType q = max / n;
+ BOOST_ASSERT(max % n != n - 1 && "c must be relatively prime to m.");
+ l1 += q * l2;
+ p = max - q * n + 1;
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);
+ return l2;
IntType q2 = n / p;
- l2 -= q2 * l1;
+ l2 += q2 * l1;
n -= q2 * p;
if(n == 0)
- return (l1 < 1 ? l1 + m : l1);
+ return m - l1;
+ q = p / n;
+ l1 += q * l2;
+ p -= q * n;
}
}
};
-
-#endif
-
} // namespace random
} // namespace boost
projects / rsem.git / blobdiff