summaryrefslogtreecommitdiffstats
path: root/contrib/src/boost/random/detail/const_mod.hpp
diff options
context:
space:
mode:
Diffstat (limited to 'contrib/src/boost/random/detail/const_mod.hpp')
-rw-r--r--contrib/src/boost/random/detail/const_mod.hpp432
1 files changed, 216 insertions, 216 deletions
diff --git a/contrib/src/boost/random/detail/const_mod.hpp b/contrib/src/boost/random/detail/const_mod.hpp
index 9f0d0e5..e0a43ab 100644
--- a/contrib/src/boost/random/detail/const_mod.hpp
+++ b/contrib/src/boost/random/detail/const_mod.hpp
@@ -1,216 +1,216 @@
-/* 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$
- *
- * Revision history
- * 2001-02-18 moved to individual header files
- */
-
-#ifndef BOOST_RANDOM_CONST_MOD_HPP
-#define BOOST_RANDOM_CONST_MOD_HPP
-
-#include <boost/assert.hpp>
-#include <boost/static_assert.hpp>
-#include <boost/integer_traits.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 {
-
-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 suppress_warnings = (m == 0);
- BOOST_ASSERT(suppress_warnings == 0);
- return x % (m + suppress_warnings);
- }
- }
-
- static IntType add(IntType x, IntType c)
- {
- 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(x < m - c)
- return x + c;
- else
- return x - (m - c);
- }
-
- static IntType mult(IntType a, IntType x)
- {
- 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
- return mult_general(a, x);
- }
-
- static IntType mult_add(IntType a, IntType x, IntType c)
- {
- 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 suppress_warnings = (m == 0);
- BOOST_ASSERT(suppress_warnings == 0);
- return (a*x+c) % (m + suppress_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 : (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 mult_small(IntType a, IntType x)
- {
- IntType suppress_warnings = (m == 0);
- BOOST_ASSERT(suppress_warnings == 0);
- return a*x % (m + suppress_warnings);
- }
-
- static IntType mult_schrage(IntType a, IntType value)
- {
- const IntType q = m / a;
- const IntType r = m % a;
-
- BOOST_ASSERT(r < q); // check that overflow cannot happen
-
- 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));
- }
- }
-
- 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_ASSERT(c > 0);
- IntType l1 = 0;
- IntType l2 = 1;
- IntType n = c;
- IntType p = m;
- for(;;) {
- IntType q = p / n;
- l1 += q * l2;
- p -= q * n;
- if(p == 0)
- return l2;
- IntType q2 = n / p;
- l2 += q2 * l1;
- n -= q2 * p;
- if(n == 0)
- return m - l1;
- }
- }
-
- // 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_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(;;) {
- if(p == 0)
- return l2;
- IntType q2 = n / p;
- l2 += q2 * l1;
- n -= q2 * p;
- if(n == 0)
- return m - l1;
- q = p / n;
- l1 += q * l2;
- p -= q * n;
- }
- }
-};
-
-} // namespace random
-} // namespace boost
-
-#include <boost/random/detail/enable_warnings.hpp>
-
-#endif // BOOST_RANDOM_CONST_MOD_HPP
+/* 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$
+ *
+ * Revision history
+ * 2001-02-18 moved to individual header files
+ */
+
+#ifndef BOOST_RANDOM_CONST_MOD_HPP
+#define BOOST_RANDOM_CONST_MOD_HPP
+
+#include <boost/assert.hpp>
+#include <boost/static_assert.hpp>
+#include <boost/integer_traits.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 {
+
+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 suppress_warnings = (m == 0);
+ BOOST_ASSERT(suppress_warnings == 0);
+ return x % (m + suppress_warnings);
+ }
+ }
+
+ static IntType add(IntType x, IntType c)
+ {
+ 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(x < m - c)
+ return x + c;
+ else
+ return x - (m - c);
+ }
+
+ static IntType mult(IntType a, IntType x)
+ {
+ 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
+ return mult_general(a, x);
+ }
+
+ static IntType mult_add(IntType a, IntType x, IntType c)
+ {
+ 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 suppress_warnings = (m == 0);
+ BOOST_ASSERT(suppress_warnings == 0);
+ return (a*x+c) % (m + suppress_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 : (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 mult_small(IntType a, IntType x)
+ {
+ IntType suppress_warnings = (m == 0);
+ BOOST_ASSERT(suppress_warnings == 0);
+ return a*x % (m + suppress_warnings);
+ }
+
+ static IntType mult_schrage(IntType a, IntType value)
+ {
+ const IntType q = m / a;
+ const IntType r = m % a;
+
+ BOOST_ASSERT(r < q); // check that overflow cannot happen
+
+ 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));
+ }
+ }
+
+ 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_ASSERT(c > 0);
+ IntType l1 = 0;
+ IntType l2 = 1;
+ IntType n = c;
+ IntType p = m;
+ for(;;) {
+ IntType q = p / n;
+ l1 += q * l2;
+ p -= q * n;
+ if(p == 0)
+ return l2;
+ IntType q2 = n / p;
+ l2 += q2 * l1;
+ n -= q2 * p;
+ if(n == 0)
+ return m - l1;
+ }
+ }
+
+ // 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_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(;;) {
+ if(p == 0)
+ return l2;
+ IntType q2 = n / p;
+ l2 += q2 * l1;
+ n -= q2 * p;
+ if(n == 0)
+ return m - l1;
+ q = p / n;
+ l1 += q * l2;
+ p -= q * n;
+ }
+ }
+};
+
+} // namespace random
+} // namespace boost
+
+#include <boost/random/detail/enable_warnings.hpp>
+
+#endif // BOOST_RANDOM_CONST_MOD_HPP