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 58649 2010-01-02 21:23:17Z 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 <cassert>
  20 #include <boost/static_assert.hpp>
  21 #include <boost/cstdint.hpp>
  22 #include <boost/integer_traits.hpp>
  23 #include <boost/detail/workaround.hpp>
  24
  25 #include <boost/random/detail/disable_warnings.hpp>
  26
  27 namespace boost {
  28 namespace random {
  29
  30 /*
  31  * Some random number generators require modular arithmetic.  Put
  32  * everything we need here.
  33  * IntType must be an integral type.
  34  */
  35
  36 namespace detail {
  37
  38   template<bool is_signed>
  39   struct do_add
  40   { };
  41
  42   template<>
  43   struct do_add<true>
  44   {
  45     template<class IntType>
  46     static IntType add(IntType m, IntType x, IntType c)
  47     {
  48       if (x < m - c)
  49         return x + c;
  50       else
  51         return x - (m-c);
  52     }
  53   };
  54
  55   template<>
  56   struct do_add<false>
  57   {
  58     template<class IntType>
  59     static IntType add(IntType, IntType, IntType)
  60     {
  61       // difficult
  62       assert(!"const_mod::add with c too large");
  63       return 0;
  64     }
  65   };
  66 } // namespace detail
  67
  68 #if !(defined(__BORLANDC__) && (__BORLANDC__ == 0x560))
  69
  70 template<class IntType, IntType m>
  71 class const_mod
  72 {
  73 public:
  74   static IntType add(IntType x, IntType c)
  75   {
  76     if(c == 0)
  77       return x;
  78     else if(c <= traits::const_max - m)    // i.e. m+c < max
  79       return add_small(x, c);
  80     else
  81       return detail::do_add<traits::is_signed>::add(m, x, c);
  82   }
  83
  84   static IntType mult(IntType a, IntType x)
  85   {
  86     if(a == 1)
  87       return x;
  88     else if(m <= traits::const_max/a)      // i.e. a*m <= max
  89       return mult_small(a, x);
  90     else if(traits::is_signed && (m%a < m/a))
  91       return mult_schrage(a, x);
  92     else {
  93       // difficult
  94       assert(!"const_mod::mult with a too large");
  95       return 0;
  96     }
  97   }
  98
  99   static IntType mult_add(IntType a, IntType x, IntType c)
 100   {
 101     if(m <= (traits::const_max-c)/a)   // i.e. a*m+c <= max
 102       return (a*x+c) % m;
 103     else
 104       return add(mult(a, x), c);
 105   }
 106
 107   static IntType invert(IntType x)
 108   { return x == 0 ? 0 : invert_euclidian(x); }
 109
 110 private:
 111   typedef integer_traits<IntType> traits;
 112
 113   const_mod();      // don't instantiate
 114
 115   static IntType add_small(IntType x, IntType c)
 116   {
 117     x += c;
 118     if(x >= m)
 119       x -= m;
 120     return x;
 121   }
 122
 123   static IntType mult_small(IntType a, IntType x)
 124   {
 125     return a*x % m;
 126   }
 127
 128   static IntType mult_schrage(IntType a, IntType value)
 129   {
 130     const IntType q = m / a;
 131     const IntType r = m % a;
 132
 133     assert(r < q);        // check that overflow cannot happen
 134
 135     value = a*(value%q) - r*(value/q);
 136     // An optimizer bug in the SGI MIPSpro 7.3.1.x compiler requires this
 137     // convoluted formulation of the loop (Synge Todo)
 138     for(;;) {
 139       if (value > 0)
 140         break;
 141       value += m;
 142     }
 143     return value;
 144   }
 145
 146   // invert c in the finite field (mod m) (m must be prime)
 147   static IntType invert_euclidian(IntType c)
 148   {
 149     // we are interested in the gcd factor for c, because this is our inverse
 150     BOOST_STATIC_ASSERT(m > 0);
 151 #if BOOST_WORKAROUND(__MWERKS__, BOOST_TESTED_AT(0x3003))
 152     assert(boost::integer_traits<IntType>::is_signed);
 153 #elif !defined(BOOST_NO_LIMITS_COMPILE_TIME_CONSTANTS)
 154     BOOST_STATIC_ASSERT(boost::integer_traits<IntType>::is_signed);
 155 #endif
 156     assert(c > 0);
 157     IntType l1 = 0;
 158     IntType l2 = 1;
 159     IntType n = c;
 160     IntType p = m;
 161     for(;;) {
 162       IntType q = p / n;
 163       l1 -= q * l2;           // this requires a signed IntType!
 164       p -= q * n;
 165       if(p == 0)
 166         return (l2 < 1 ? l2 + m : l2);
 167       IntType q2 = n / p;
 168       l2 -= q2 * l1;
 169       n -= q2 * p;
 170       if(n == 0)
 171         return (l1 < 1 ? l1 + m : l1);
 172     }
 173   }
 174 };
 175
 176 // The modulus is exactly the word size: rely on machine overflow handling.
 177 // Due to a GCC bug, we cannot partially specialize in the presence of
 178 // template value parameters.
 179 template<>
 180 class const_mod<unsigned int, 0>
 181 {
 182   typedef unsigned int IntType;
 183 public:
 184   static IntType add(IntType x, IntType c) { return x+c; }
 185   static IntType mult(IntType a, IntType x) { return a*x; }
 186   static IntType mult_add(IntType a, IntType x, IntType c) { return a*x+c; }
 187
 188   // m is not prime, thus invert is not useful
 189 private:                      // don't instantiate
 190   const_mod();
 191 };
 192
 193 template<>
 194 class const_mod<unsigned long, 0>
 195 {
 196   typedef unsigned long IntType;
 197 public:
 198   static IntType add(IntType x, IntType c) { return x+c; }
 199   static IntType mult(IntType a, IntType x) { return a*x; }
 200   static IntType mult_add(IntType a, IntType x, IntType c) { return a*x+c; }
 201
 202   // m is not prime, thus invert is not useful
 203 private:                      // don't instantiate
 204   const_mod();
 205 };
 206
 207 // the modulus is some power of 2: rely partly on machine overflow handling
 208 // we only specialize for rand48 at the moment
 209 #ifndef BOOST_NO_INT64_T
 210 template<>
 211 class const_mod<uint64_t, uint64_t(1) << 48>
 212 {
 213   typedef uint64_t IntType;
 214 public:
 215   static IntType add(IntType x, IntType c) { return c == 0 ? x : mod(x+c); }
 216   static IntType mult(IntType a, IntType x) { return mod(a*x); }
 217   static IntType mult_add(IntType a, IntType x, IntType c)
 218     { return mod(a*x+c); }
 219   static IntType mod(IntType x) { return x &= ((uint64_t(1) << 48)-1); }
 220
 221   // m is not prime, thus invert is not useful
 222 private:                      // don't instantiate
 223   const_mod();
 224 };
 225 #endif /* !BOOST_NO_INT64_T */
 226
 227 #else
 228
 229 //
 230 // for some reason Borland C++ Builder 6 has problems with
 231 // the full specialisations of const_mod, define a generic version
 232 // instead, the compiler will optimise away the const-if statements:
 233 //
 234
 235 template<class IntType, IntType m>
 236 class const_mod
 237 {
 238 public:
 239   static IntType add(IntType x, IntType c)
 240   {
 241     if(0 == m)
 242     {
 243        return x+c;
 244     }
 245     else
 246     {
 247        if(c == 0)
 248          return x;
 249        else if(c <= traits::const_max - m)    // i.e. m+c < max
 250          return add_small(x, c);
 251        else
 252          return detail::do_add<traits::is_signed>::add(m, x, c);
 253     }
 254   }
 255
 256   static IntType mult(IntType a, IntType x)
 257   {
 258     if(x == 0)
 259     {
 260        return a*x;
 261     }
 262     else
 263     {
 264        if(a == 1)
 265          return x;
 266        else if(m <= traits::const_max/a)      // i.e. a*m <= max
 267          return mult_small(a, x);
 268        else if(traits::is_signed && (m%a < m/a))
 269          return mult_schrage(a, x);
 270        else {
 271          // difficult
 272          assert(!"const_mod::mult with a too large");
 273          return 0;
 274        }
 275     }
 276   }
 277
 278   static IntType mult_add(IntType a, IntType x, IntType c)
 279   {
 280     if(m == 0)
 281     {
 282        return a*x+c;
 283     }
 284     else
 285     {
 286        if(m <= (traits::const_max-c)/a)   // i.e. a*m+c <= max
 287          return (a*x+c) % m;
 288        else
 289          return add(mult(a, x), c);
 290     }
 291   }
 292
 293   static IntType invert(IntType x)
 294   { return x == 0 ? 0 : invert_euclidian(x); }
 295
 296 private:
 297   typedef integer_traits<IntType> traits;
 298
 299   const_mod();      // don't instantiate
 300
 301   static IntType add_small(IntType x, IntType c)
 302   {
 303     x += c;
 304     if(x >= m)
 305       x -= m;
 306     return x;
 307   }
 308
 309   static IntType mult_small(IntType a, IntType x)
 310   {
 311     return a*x % m;
 312   }
 313
 314   static IntType mult_schrage(IntType a, IntType value)
 315   {
 316     const IntType q = m / a;
 317     const IntType r = m % a;
 318
 319     assert(r < q);        // check that overflow cannot happen
 320
 321     value = a*(value%q) - r*(value/q);
 322     while(value <= 0)
 323       value += m;
 324     return value;
 325   }
 326
 327   // invert c in the finite field (mod m) (m must be prime)
 328   static IntType invert_euclidian(IntType c)
 329   {
 330     // we are interested in the gcd factor for c, because this is our inverse
 331     BOOST_STATIC_ASSERT(m > 0);
 332 #ifndef BOOST_NO_LIMITS_COMPILE_TIME_CONSTANTS
 333     BOOST_STATIC_ASSERT(boost::integer_traits<IntType>::is_signed);
 334 #endif
 335     assert(c > 0);
 336     IntType l1 = 0;
 337     IntType l2 = 1;
 338     IntType n = c;
 339     IntType p = m;
 340     for(;;) {
 341       IntType q = p / n;
 342       l1 -= q * l2;           // this requires a signed IntType!
 343       p -= q * n;
 344       if(p == 0)
 345         return (l2 < 1 ? l2 + m : l2);
 346       IntType q2 = n / p;
 347       l2 -= q2 * l1;
 348       n -= q2 * p;
 349       if(n == 0)
 350         return (l1 < 1 ? l1 + m : l1);
 351     }
 352   }
 353 };
 354
 355
 356 #endif
 357
 358 } // namespace random
 359 } // namespace boost
 360
 361 #include <boost/random/detail/enable_warnings.hpp>
 362
 363 #endif // BOOST_RANDOM_CONST_MOD_HPP