stc/exp/ml_stc_linux_make_v1.0/include/boost/random/detail/const_mod.hpp

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