Boost-Graph
view release on metacpan or search on metacpan
include/boost/random/detail/const_mod.hpp view on Meta::CPAN
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:
( run in 1.168 second using v1.01-cache-2.11-cpan-39bf76dae61 )