williamr@2: /* boost random/detail/const_mod.hpp header file williamr@2: * williamr@2: * Copyright Jens Maurer 2000-2001 williamr@2: * Distributed under the Boost Software License, Version 1.0. (See williamr@2: * accompanying file LICENSE_1_0.txt or copy at williamr@2: * http://www.boost.org/LICENSE_1_0.txt) williamr@2: * williamr@2: * See http://www.boost.org for most recent version including documentation. williamr@2: * williamr@2: * $Id: const_mod.hpp,v 1.8 2004/07/27 03:43:32 dgregor Exp $ williamr@2: * williamr@2: * Revision history williamr@2: * 2001-02-18 moved to individual header files williamr@2: */ williamr@2: williamr@2: #ifndef BOOST_RANDOM_CONST_MOD_HPP williamr@2: #define BOOST_RANDOM_CONST_MOD_HPP williamr@2: williamr@2: #include williamr@2: #include williamr@2: #include williamr@2: #include williamr@2: #include williamr@2: williamr@2: namespace boost { williamr@2: namespace random { williamr@2: williamr@2: /* williamr@2: * Some random number generators require modular arithmetic. Put williamr@2: * everything we need here. williamr@2: * IntType must be an integral type. williamr@2: */ williamr@2: williamr@2: namespace detail { williamr@2: williamr@2: template williamr@2: struct do_add williamr@2: { }; williamr@2: williamr@2: template<> williamr@2: struct do_add williamr@2: { williamr@2: template williamr@2: static IntType add(IntType m, IntType x, IntType c) williamr@2: { williamr@2: x += (c-m); williamr@2: if(x < 0) williamr@2: x += m; williamr@2: return x; williamr@2: } williamr@2: }; williamr@2: williamr@2: template<> williamr@2: struct do_add williamr@2: { williamr@2: template williamr@2: static IntType add(IntType, IntType, IntType) williamr@2: { williamr@2: // difficult williamr@2: assert(!"const_mod::add with c too large"); williamr@2: return 0; williamr@2: } williamr@2: }; williamr@2: } // namespace detail williamr@2: williamr@2: #if !(defined(__BORLANDC__) && (__BORLANDC__ == 0x560)) williamr@2: williamr@2: template williamr@2: class const_mod williamr@2: { williamr@2: public: williamr@2: static IntType add(IntType x, IntType c) williamr@2: { williamr@2: if(c == 0) williamr@2: return x; williamr@2: else if(c <= traits::const_max - m) // i.e. m+c < max williamr@2: return add_small(x, c); williamr@2: else williamr@2: return detail::do_add::add(m, x, c); williamr@2: } williamr@2: williamr@2: static IntType mult(IntType a, IntType x) williamr@2: { williamr@2: if(a == 1) williamr@2: return x; williamr@2: else if(m <= traits::const_max/a) // i.e. a*m <= max williamr@2: return mult_small(a, x); williamr@2: else if(traits::is_signed && (m%a < m/a)) williamr@2: return mult_schrage(a, x); williamr@2: else { williamr@2: // difficult williamr@2: assert(!"const_mod::mult with a too large"); williamr@2: return 0; williamr@2: } williamr@2: } williamr@2: williamr@2: static IntType mult_add(IntType a, IntType x, IntType c) williamr@2: { williamr@2: if(m <= (traits::const_max-c)/a) // i.e. a*m+c <= max williamr@2: return (a*x+c) % m; williamr@2: else williamr@2: return add(mult(a, x), c); williamr@2: } williamr@2: williamr@2: static IntType invert(IntType x) williamr@2: { return x == 0 ? 0 : invert_euclidian(x); } williamr@2: williamr@2: private: williamr@2: typedef integer_traits traits; williamr@2: williamr@2: const_mod(); // don't instantiate williamr@2: williamr@2: static IntType add_small(IntType x, IntType c) williamr@2: { williamr@2: x += c; williamr@2: if(x >= m) williamr@2: x -= m; williamr@2: return x; williamr@2: } williamr@2: williamr@2: static IntType mult_small(IntType a, IntType x) williamr@2: { williamr@2: return a*x % m; williamr@2: } williamr@2: williamr@2: static IntType mult_schrage(IntType a, IntType value) williamr@2: { williamr@2: const IntType q = m / a; williamr@2: const IntType r = m % a; williamr@2: williamr@2: assert(r < q); // check that overflow cannot happen williamr@2: williamr@2: value = a*(value%q) - r*(value/q); williamr@2: // An optimizer bug in the SGI MIPSpro 7.3.1.x compiler requires this williamr@2: // convoluted formulation of the loop (Synge Todo) williamr@2: for(;;) { williamr@2: if (value > 0) williamr@2: break; williamr@2: value += m; williamr@2: } williamr@2: return value; williamr@2: } williamr@2: williamr@2: // invert c in the finite field (mod m) (m must be prime) williamr@2: static IntType invert_euclidian(IntType c) williamr@2: { williamr@2: // we are interested in the gcd factor for c, because this is our inverse williamr@2: BOOST_STATIC_ASSERT(m > 0); williamr@2: #if BOOST_WORKAROUND(__MWERKS__, BOOST_TESTED_AT(0x3003)) williamr@2: assert(boost::integer_traits::is_signed); williamr@2: #elif !defined(BOOST_NO_LIMITS_COMPILE_TIME_CONSTANTS) williamr@2: BOOST_STATIC_ASSERT(boost::integer_traits::is_signed); williamr@2: #endif williamr@2: assert(c > 0); williamr@2: IntType l1 = 0; williamr@2: IntType l2 = 1; williamr@2: IntType n = c; williamr@2: IntType p = m; williamr@2: for(;;) { williamr@2: IntType q = p / n; williamr@2: l1 -= q * l2; // this requires a signed IntType! williamr@2: p -= q * n; williamr@2: if(p == 0) williamr@2: return (l2 < 1 ? l2 + m : l2); williamr@2: IntType q2 = n / p; williamr@2: l2 -= q2 * l1; williamr@2: n -= q2 * p; williamr@2: if(n == 0) williamr@2: return (l1 < 1 ? l1 + m : l1); williamr@2: } williamr@2: } williamr@2: }; williamr@2: williamr@2: // The modulus is exactly the word size: rely on machine overflow handling. williamr@2: // Due to a GCC bug, we cannot partially specialize in the presence of williamr@2: // template value parameters. williamr@2: template<> williamr@2: class const_mod williamr@2: { williamr@2: typedef unsigned int IntType; williamr@2: public: williamr@2: static IntType add(IntType x, IntType c) { return x+c; } williamr@2: static IntType mult(IntType a, IntType x) { return a*x; } williamr@2: static IntType mult_add(IntType a, IntType x, IntType c) { return a*x+c; } williamr@2: williamr@2: // m is not prime, thus invert is not useful williamr@2: private: // don't instantiate williamr@2: const_mod(); williamr@2: }; williamr@2: williamr@2: template<> williamr@2: class const_mod williamr@2: { williamr@2: typedef unsigned long IntType; williamr@2: public: williamr@2: static IntType add(IntType x, IntType c) { return x+c; } williamr@2: static IntType mult(IntType a, IntType x) { return a*x; } williamr@2: static IntType mult_add(IntType a, IntType x, IntType c) { return a*x+c; } williamr@2: williamr@2: // m is not prime, thus invert is not useful williamr@2: private: // don't instantiate williamr@2: const_mod(); williamr@2: }; williamr@2: williamr@2: // the modulus is some power of 2: rely partly on machine overflow handling williamr@2: // we only specialize for rand48 at the moment williamr@2: #ifndef BOOST_NO_INT64_T williamr@2: template<> williamr@2: class const_mod williamr@2: { williamr@2: typedef uint64_t IntType; williamr@2: public: williamr@2: static IntType add(IntType x, IntType c) { return c == 0 ? x : mod(x+c); } williamr@2: static IntType mult(IntType a, IntType x) { return mod(a*x); } williamr@2: static IntType mult_add(IntType a, IntType x, IntType c) williamr@2: { return mod(a*x+c); } williamr@2: static IntType mod(IntType x) { return x &= ((uint64_t(1) << 48)-1); } williamr@2: williamr@2: // m is not prime, thus invert is not useful williamr@2: private: // don't instantiate williamr@2: const_mod(); williamr@2: }; williamr@2: #endif /* !BOOST_NO_INT64_T */ williamr@2: williamr@2: #else williamr@2: williamr@2: // williamr@2: // for some reason Borland C++ Builder 6 has problems with williamr@2: // the full specialisations of const_mod, define a generic version williamr@2: // instead, the compiler will optimise away the const-if statements: williamr@2: // williamr@2: williamr@2: template williamr@2: class const_mod williamr@2: { williamr@2: public: williamr@2: static IntType add(IntType x, IntType c) williamr@2: { williamr@2: if(0 == m) williamr@2: { williamr@2: return x+c; williamr@2: } williamr@2: else williamr@2: { williamr@2: if(c == 0) williamr@2: return x; williamr@2: else if(c <= traits::const_max - m) // i.e. m+c < max williamr@2: return add_small(x, c); williamr@2: else williamr@2: return detail::do_add::add(m, x, c); williamr@2: } williamr@2: } williamr@2: williamr@2: static IntType mult(IntType a, IntType x) williamr@2: { williamr@2: if(x == 0) williamr@2: { williamr@2: return a*x; williamr@2: } williamr@2: else williamr@2: { williamr@2: if(a == 1) williamr@2: return x; williamr@2: else if(m <= traits::const_max/a) // i.e. a*m <= max williamr@2: return mult_small(a, x); williamr@2: else if(traits::is_signed && (m%a < m/a)) williamr@2: return mult_schrage(a, x); williamr@2: else { williamr@2: // difficult williamr@2: assert(!"const_mod::mult with a too large"); williamr@2: return 0; williamr@2: } williamr@2: } williamr@2: } williamr@2: williamr@2: static IntType mult_add(IntType a, IntType x, IntType c) williamr@2: { williamr@2: if(m == 0) williamr@2: { williamr@2: return a*x+c; williamr@2: } williamr@2: else williamr@2: { williamr@2: if(m <= (traits::const_max-c)/a) // i.e. a*m+c <= max williamr@2: return (a*x+c) % m; williamr@2: else williamr@2: return add(mult(a, x), c); williamr@2: } williamr@2: } williamr@2: williamr@2: static IntType invert(IntType x) williamr@2: { return x == 0 ? 0 : invert_euclidian(x); } williamr@2: williamr@2: private: williamr@2: typedef integer_traits traits; williamr@2: williamr@2: const_mod(); // don't instantiate williamr@2: williamr@2: static IntType add_small(IntType x, IntType c) williamr@2: { williamr@2: x += c; williamr@2: if(x >= m) williamr@2: x -= m; williamr@2: return x; williamr@2: } williamr@2: williamr@2: static IntType mult_small(IntType a, IntType x) williamr@2: { williamr@2: return a*x % m; williamr@2: } williamr@2: williamr@2: static IntType mult_schrage(IntType a, IntType value) williamr@2: { williamr@2: const IntType q = m / a; williamr@2: const IntType r = m % a; williamr@2: williamr@2: assert(r < q); // check that overflow cannot happen williamr@2: williamr@2: value = a*(value%q) - r*(value/q); williamr@2: while(value <= 0) williamr@2: value += m; williamr@2: return value; williamr@2: } williamr@2: williamr@2: // invert c in the finite field (mod m) (m must be prime) williamr@2: static IntType invert_euclidian(IntType c) williamr@2: { williamr@2: // we are interested in the gcd factor for c, because this is our inverse williamr@2: BOOST_STATIC_ASSERT(m > 0); williamr@2: #ifndef BOOST_NO_LIMITS_COMPILE_TIME_CONSTANTS williamr@2: BOOST_STATIC_ASSERT(boost::integer_traits::is_signed); williamr@2: #endif williamr@2: assert(c > 0); williamr@2: IntType l1 = 0; williamr@2: IntType l2 = 1; williamr@2: IntType n = c; williamr@2: IntType p = m; williamr@2: for(;;) { williamr@2: IntType q = p / n; williamr@2: l1 -= q * l2; // this requires a signed IntType! williamr@2: p -= q * n; williamr@2: if(p == 0) williamr@2: return (l2 < 1 ? l2 + m : l2); williamr@2: IntType q2 = n / p; williamr@2: l2 -= q2 * l1; williamr@2: n -= q2 * p; williamr@2: if(n == 0) williamr@2: return (l1 < 1 ? l1 + m : l1); williamr@2: } williamr@2: } williamr@2: }; williamr@2: williamr@2: williamr@2: #endif williamr@2: williamr@2: } // namespace random williamr@2: } // namespace boost williamr@2: williamr@2: #endif // BOOST_RANDOM_CONST_MOD_HPP