1 /* boost random/detail/const_mod.hpp header file 2 * 3 * Copyright Jens Maurer 2000-2001 4 * Distributed under the Boost Software License, Version 1.0. (See 5 * accompanying file LICENSE_1_0.txt or copy at 6 * http://www.boost.org/LICENSE_1_0.txt) 7 * 8 * See http://www.boost.org for most recent version including documentation. 9 * 10 * $Id$ 11 * 12 * Revision history 13 * 2001-02-18 moved to individual header files 14 */ 15 16 #ifndef BOOST_RANDOM_CONST_MOD_HPP 17 #define BOOST_RANDOM_CONST_MOD_HPP 18 19 #include <boost/assert.hpp> 20 #include <boost/static_assert.hpp> 21 #include <boost/integer_traits.hpp> 22 #include <boost/type_traits/make_unsigned.hpp> 23 #include <boost/random/detail/large_arithmetic.hpp> 24 25 #include <boost/random/detail/disable_warnings.hpp> 26 27 namespace boost { 28 namespace random { 29 30 template<class IntType, IntType m> 31 class const_mod 32 { 33 public: apply(IntType x)34 static IntType apply(IntType x) 35 { 36 if(((unsigned_m() - 1) & unsigned_m()) == 0) 37 return (unsigned_type(x)) & (unsigned_m() - 1); 38 else { 39 IntType suppress_warnings = (m == 0); 40 BOOST_ASSERT(suppress_warnings == 0); 41 return x % (m + suppress_warnings); 42 } 43 } 44 add(IntType x,IntType c)45 static IntType add(IntType x, IntType c) 46 { 47 if(((unsigned_m() - 1) & unsigned_m()) == 0) 48 return (unsigned_type(x) + unsigned_type(c)) & (unsigned_m() - 1); 49 else if(c == 0) 50 return x; 51 else if(x < m - c) 52 return x + c; 53 else 54 return x - (m - c); 55 } 56 mult(IntType a,IntType x)57 static IntType mult(IntType a, IntType x) 58 { 59 if(((unsigned_m() - 1) & unsigned_m()) == 0) 60 return unsigned_type(a) * unsigned_type(x) & (unsigned_m() - 1); 61 else if(a == 0) 62 return 0; 63 else if(a == 1) 64 return x; 65 else if(m <= traits::const_max/a) // i.e. a*m <= max 66 return mult_small(a, x); 67 else if(traits::is_signed && (m%a < m/a)) 68 return mult_schrage(a, x); 69 else 70 return mult_general(a, x); 71 } 72 mult_add(IntType a,IntType x,IntType c)73 static IntType mult_add(IntType a, IntType x, IntType c) 74 { 75 if(((unsigned_m() - 1) & unsigned_m()) == 0) 76 return (unsigned_type(a) * unsigned_type(x) + unsigned_type(c)) & (unsigned_m() - 1); 77 else if(a == 0) 78 return c; 79 else if(m <= (traits::const_max-c)/a) { // i.e. a*m+c <= max 80 IntType suppress_warnings = (m == 0); 81 BOOST_ASSERT(suppress_warnings == 0); 82 return (a*x+c) % (m + suppress_warnings); 83 } else 84 return add(mult(a, x), c); 85 } 86 pow(IntType a,boost::uintmax_t exponent)87 static IntType pow(IntType a, boost::uintmax_t exponent) 88 { 89 IntType result = 1; 90 while(exponent != 0) { 91 if(exponent % 2 == 1) { 92 result = mult(result, a); 93 } 94 a = mult(a, a); 95 exponent /= 2; 96 } 97 return result; 98 } 99 invert(IntType x)100 static IntType invert(IntType x) 101 { return x == 0 ? 0 : (m == 0? invert_euclidian0(x) : invert_euclidian(x)); } 102 103 private: 104 typedef integer_traits<IntType> traits; 105 typedef typename make_unsigned<IntType>::type unsigned_type; 106 107 const_mod(); // don't instantiate 108 mult_small(IntType a,IntType x)109 static IntType mult_small(IntType a, IntType x) 110 { 111 IntType suppress_warnings = (m == 0); 112 BOOST_ASSERT(suppress_warnings == 0); 113 return a*x % (m + suppress_warnings); 114 } 115 mult_schrage(IntType a,IntType value)116 static IntType mult_schrage(IntType a, IntType value) 117 { 118 const IntType q = m / a; 119 const IntType r = m % a; 120 121 BOOST_ASSERT(r < q); // check that overflow cannot happen 122 123 return sub(a*(value%q), r*(value/q)); 124 } 125 mult_general(IntType a,IntType b)126 static IntType mult_general(IntType a, IntType b) 127 { 128 IntType suppress_warnings = (m == 0); 129 BOOST_ASSERT(suppress_warnings == 0); 130 IntType modulus = m + suppress_warnings; 131 BOOST_ASSERT(modulus == m); 132 if(::boost::uintmax_t(modulus) <= 133 (::std::numeric_limits< ::boost::uintmax_t>::max)() / modulus) 134 { 135 return static_cast<IntType>(boost::uintmax_t(a) * b % modulus); 136 } else { 137 return static_cast<IntType>(detail::mulmod(a, b, modulus)); 138 } 139 } 140 sub(IntType a,IntType b)141 static IntType sub(IntType a, IntType b) 142 { 143 if(a < b) 144 return m - (b - a); 145 else 146 return a - b; 147 } 148 unsigned_m()149 static unsigned_type unsigned_m() 150 { 151 if(m == 0) { 152 return unsigned_type((std::numeric_limits<IntType>::max)()) + 1; 153 } else { 154 return unsigned_type(m); 155 } 156 } 157 158 // invert c in the finite field (mod m) (m must be prime) invert_euclidian(IntType c)159 static IntType invert_euclidian(IntType c) 160 { 161 // we are interested in the gcd factor for c, because this is our inverse 162 BOOST_ASSERT(c > 0); 163 IntType l1 = 0; 164 IntType l2 = 1; 165 IntType n = c; 166 IntType p = m; 167 for(;;) { 168 IntType q = p / n; 169 l1 += q * l2; 170 p -= q * n; 171 if(p == 0) 172 return l2; 173 IntType q2 = n / p; 174 l2 += q2 * l1; 175 n -= q2 * p; 176 if(n == 0) 177 return m - l1; 178 } 179 } 180 181 // invert c in the finite field (mod m) (c must be relatively prime to m) invert_euclidian0(IntType c)182 static IntType invert_euclidian0(IntType c) 183 { 184 // we are interested in the gcd factor for c, because this is our inverse 185 BOOST_ASSERT(c > 0); 186 if(c == 1) return 1; 187 IntType l1 = 0; 188 IntType l2 = 1; 189 IntType n = c; 190 IntType p = m; 191 IntType max = (std::numeric_limits<IntType>::max)(); 192 IntType q = max / n; 193 BOOST_ASSERT(max % n != n - 1 && "c must be relatively prime to m."); 194 l1 += q * l2; 195 p = max - q * n + 1; 196 for(;;) { 197 if(p == 0) 198 return l2; 199 IntType q2 = n / p; 200 l2 += q2 * l1; 201 n -= q2 * p; 202 if(n == 0) 203 return m - l1; 204 q = p / n; 205 l1 += q * l2; 206 p -= q * n; 207 } 208 } 209 }; 210 211 } // namespace random 212 } // namespace boost 213 214 #include <boost/random/detail/enable_warnings.hpp> 215 216 #endif // BOOST_RANDOM_CONST_MOD_HPP 217