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