| /* boost random/detail/const_mod.hpp header file |
| * |
| * Copyright Jens Maurer 2000-2001 |
| * 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: const_mod.hpp 58649 2010-01-02 21:23:17Z 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/static_assert.hpp> |
| #include <boost/cstdint.hpp> |
| #include <boost/integer_traits.hpp> |
| #include <boost/detail/workaround.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 add(IntType x, IntType c) |
| { |
| 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(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 <= (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); |
| // 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 value; |
| } |
| |
| // 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); |
| 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! |
| p -= q * n; |
| if(p == 0) |
| return (l2 < 1 ? l2 + m : l2); |
| IntType q2 = n / p; |
| l2 -= q2 * l1; |
| n -= q2 * p; |
| if(n == 0) |
| return (l1 < 1 ? l1 + 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) |
| { |
| // 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); |
| 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! |
| p -= q * n; |
| if(p == 0) |
| return (l2 < 1 ? l2 + m : l2); |
| IntType q2 = n / p; |
| l2 -= q2 * l1; |
| n -= q2 * p; |
| if(n == 0) |
| return (l1 < 1 ? l1 + m : l1); |
| } |
| } |
| }; |
| |
| |
| #endif |
| |
| } // namespace random |
| } // namespace boost |
| |
| #include <boost/random/detail/enable_warnings.hpp> |
| |
| #endif // BOOST_RANDOM_CONST_MOD_HPP |