random/detail/const_mod.hpp

0001 /* boost random/detail/const_mod.hpp header file
0002  *
0003  * Copyright Jens Maurer 2000-2001
0004  * Distributed under the Boost Software License, Version 1.0. (See
0005  * accompanying file LICENSE_1_0.txt or copy at
0006  * http://www.boost.org/LICENSE_1_0.txt)
0007  *
0008  * See http://www.boost.org for most recent version including documentation.
0009  *
0010  * $Id$
0011  *
0012  * Revision history
0013  *  2001-02-18  moved to individual header files
0014  */
0015
0016 #ifndef BOOST_RANDOM_CONST_MOD_HPP
0017 #define BOOST_RANDOM_CONST_MOD_HPP
0018
0019 #include <boost/assert.hpp>
0020 #include <boost/static_assert.hpp>
0021 #include <boost/integer_traits.hpp>
0022 #include <boost/type_traits/make_unsigned.hpp>
0023 #include <boost/random/detail/large_arithmetic.hpp>
0024
0025 #include <boost/random/detail/disable_warnings.hpp>
0026
0027 namespace boost {
0028 namespace random {
0029
0030 template<class IntType, IntType m>
0031 class const_mod
0032 {
0033 public:
0034   static IntType apply(IntType x)
0035   {
0036     if(((unsigned_m() - 1) & unsigned_m()) == 0)
0037       return (unsigned_type(x)) & (unsigned_m() - 1);
0038     else {
0039       IntType suppress_warnings = (m == 0);
0040       BOOST_ASSERT(suppress_warnings == 0);
0041       return x % (m + suppress_warnings);
0042     }
0043   }
0044
0045   static IntType add(IntType x, IntType c)
0046   {
0047     if(((unsigned_m() - 1) & unsigned_m()) == 0)
0048       return (unsigned_type(x) + unsigned_type(c)) & (unsigned_m() - 1);
0049     else if(c == 0)
0050       return x;
0051     else if(x < m - c)
0052       return x + c;
0053     else
0054       return x - (m - c);
0055   }
0056
0057   static IntType mult(IntType a, IntType x)
0058   {
0059     if(((unsigned_m() - 1) & unsigned_m()) == 0)
0060       return unsigned_type(a) * unsigned_type(x) & (unsigned_m() - 1);
0061     else if(a == 0)
0062       return 0;
0063     else if(a == 1)
0064       return x;
0065     else if(m <= traits::const_max/a)      // i.e. a*m <= max
0066       return mult_small(a, x);
0067     else if(traits::is_signed && (m%a < m/a))
0068       return mult_schrage(a, x);
0069     else
0070       return mult_general(a, x);
0071   }
0072
0073   static IntType mult_add(IntType a, IntType x, IntType c)
0074   {
0075     if(((unsigned_m() - 1) & unsigned_m()) == 0)
0076       return (unsigned_type(a) * unsigned_type(x) + unsigned_type(c)) & (unsigned_m() - 1);
0077     else if(a == 0)
0078       return c;
0079     else if(m <= (traits::const_max-c)/a) {  // i.e. a*m+c <= max
0080       IntType suppress_warnings = (m == 0);
0081       BOOST_ASSERT(suppress_warnings == 0);
0082       return (a*x+c) % (m + suppress_warnings);
0083     } else
0084       return add(mult(a, x), c);
0085   }
0086
0087   static IntType pow(IntType a, boost::uintmax_t exponent)
0088   {
0089       IntType result = 1;
0090       while(exponent != 0) {
0091           if(exponent % 2 == 1) {
0092               result = mult(result, a);
0093           }
0094           a = mult(a, a);
0095           exponent /= 2;
0096       }
0097       return result;
0098   }
0099
0100   static IntType invert(IntType x)
0101   { return x == 0 ? 0 : (m == 0? invert_euclidian0(x) : invert_euclidian(x)); }
0102
0103 private:
0104   typedef integer_traits<IntType> traits;
0105   typedef typename make_unsigned<IntType>::type unsigned_type;
0106
0107   const_mod();      // don't instantiate
0108
0109   static IntType mult_small(IntType a, IntType x)
0110   {
0111     IntType suppress_warnings = (m == 0);
0112     BOOST_ASSERT(suppress_warnings == 0);
0113     return a*x % (m + suppress_warnings);
0114   }
0115
0116   static IntType mult_schrage(IntType a, IntType value)
0117   {
0118     const IntType q = m / a;
0119     const IntType r = m % a;
0120
0121     BOOST_ASSERT(r < q);        // check that overflow cannot happen
0122
0123     return sub(a*(value%q), r*(value/q));
0124   }
0125
0126   static IntType mult_general(IntType a, IntType b)
0127   {
0128     IntType suppress_warnings = (m == 0);
0129     BOOST_ASSERT(suppress_warnings == 0);
0130     IntType modulus = m + suppress_warnings;
0131     BOOST_ASSERT(modulus == m);
0132     if(::boost::uintmax_t(modulus) <=
0133         (::std::numeric_limits< ::boost::uintmax_t>::max)() / modulus)
0134     {
0135       return static_cast<IntType>(boost::uintmax_t(a) * b % modulus);
0136     } else {
0137       return static_cast<IntType>(detail::mulmod(a, b, modulus));
0138     }
0139   }
0140
0141   static IntType sub(IntType a, IntType b)
0142   {
0143     if(a < b)
0144       return m - (b - a);
0145     else
0146       return a - b;
0147   }
0148
0149   static unsigned_type unsigned_m()
0150   {
0151       if(m == 0) {
0152           return unsigned_type((std::numeric_limits<IntType>::max)()) + 1;
0153       } else {
0154           return unsigned_type(m);
0155       }
0156   }
0157
0158   // invert c in the finite field (mod m) (m must be prime)
0159   static IntType invert_euclidian(IntType c)
0160   {
0161     // we are interested in the gcd factor for c, because this is our inverse
0162     BOOST_ASSERT(c > 0);
0163     IntType l1 = 0;
0164     IntType l2 = 1;
0165     IntType n = c;
0166     IntType p = m;
0167     for(;;) {
0168       IntType q = p / n;
0169       l1 += q * l2;
0170       p -= q * n;
0171       if(p == 0)
0172         return l2;
0173       IntType q2 = n / p;
0174       l2 += q2 * l1;
0175       n -= q2 * p;
0176       if(n == 0)
0177         return m - l1;
0178     }
0179   }
0180
0181   // invert c in the finite field (mod m) (c must be relatively prime to m)
0182   static IntType invert_euclidian0(IntType c)
0183   {
0184     // we are interested in the gcd factor for c, because this is our inverse
0185     BOOST_ASSERT(c > 0);
0186     if(c == 1) return 1;
0187     IntType l1 = 0;
0188     IntType l2 = 1;
0189     IntType n = c;
0190     IntType p = m;
0191     IntType max = (std::numeric_limits<IntType>::max)();
0192     IntType q = max / n;
0193     BOOST_ASSERT(max % n != n - 1 && "c must be relatively prime to m.");
0194     l1 += q * l2;
0195     p = max - q * n + 1;
0196     for(;;) {
0197       if(p == 0)
0198         return l2;
0199       IntType q2 = n / p;
0200       l2 += q2 * l1;
0201       n -= q2 * p;
0202       if(n == 0)
0203         return m - l1;
0204       q = p / n;
0205       l1 += q * l2;
0206       p -= q * n;
0207     }
0208   }
0209 };
0210
0211 } // namespace random
0212 } // namespace boost
0213
0214 #include <boost/random/detail/enable_warnings.hpp>
0215
0216 #endif // BOOST_RANDOM_CONST_MOD_HPP