boost/random/detail/const_mod.hpp

   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: const_mod.hpp 71018 2011-04-05 21:27:52Z steven_watanabe $
  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:
  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 supress_warnings = (m == 0);
  40       BOOST_ASSERT(supress_warnings == 0);
  41       return x % (m + supress_warnings);
  42     }
  43   }
  44
  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
  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
  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 supress_warnings = (m == 0);
  81       BOOST_ASSERT(supress_warnings == 0);
  82       return (a*x+c) % (m + supress_warnings);
  83     } else
  84       return add(mult(a, x), c);
  85   }
  86
  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
 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
 109   static IntType mult_small(IntType a, IntType x)
 110   {
 111     IntType supress_warnings = (m == 0);
 112     BOOST_ASSERT(supress_warnings == 0);
 113     return a*x % (m + supress_warnings);
 114   }
 115
 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
 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
 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
 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)
 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)
 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