Chris@16: /* boost random/detail/const_mod.hpp header file
Chris@16:  *
Chris@16:  * Copyright Jens Maurer 2000-2001
Chris@16:  * Distributed under the Boost Software License, Version 1.0. (See
Chris@16:  * accompanying file LICENSE_1_0.txt or copy at
Chris@16:  * http://www.boost.org/LICENSE_1_0.txt)
Chris@16:  *
Chris@16:  * See http://www.boost.org for most recent version including documentation.
Chris@16:  *
Chris@101:  * $Id$
Chris@16:  *
Chris@16:  * Revision history
Chris@16:  *  2001-02-18  moved to individual header files
Chris@16:  */
Chris@16: 
Chris@16: #ifndef BOOST_RANDOM_CONST_MOD_HPP
Chris@16: #define BOOST_RANDOM_CONST_MOD_HPP
Chris@16: 
Chris@16: #include <boost/assert.hpp>
Chris@16: #include <boost/static_assert.hpp>
Chris@16: #include <boost/integer_traits.hpp>
Chris@16: #include <boost/type_traits/make_unsigned.hpp>
Chris@16: #include <boost/random/detail/large_arithmetic.hpp>
Chris@16: 
Chris@16: #include <boost/random/detail/disable_warnings.hpp>
Chris@16: 
Chris@16: namespace boost {
Chris@16: namespace random {
Chris@16: 
Chris@16: template<class IntType, IntType m>
Chris@16: class const_mod
Chris@16: {
Chris@16: public:
Chris@16:   static IntType apply(IntType x)
Chris@16:   {
Chris@16:     if(((unsigned_m() - 1) & unsigned_m()) == 0)
Chris@16:       return (unsigned_type(x)) & (unsigned_m() - 1);
Chris@16:     else {
Chris@101:       IntType suppress_warnings = (m == 0);
Chris@101:       BOOST_ASSERT(suppress_warnings == 0);
Chris@101:       return x % (m + suppress_warnings);
Chris@16:     }
Chris@16:   }
Chris@16: 
Chris@16:   static IntType add(IntType x, IntType c)
Chris@16:   {
Chris@16:     if(((unsigned_m() - 1) & unsigned_m()) == 0)
Chris@16:       return (unsigned_type(x) + unsigned_type(c)) & (unsigned_m() - 1);
Chris@16:     else if(c == 0)
Chris@16:       return x;
Chris@16:     else if(x < m - c)
Chris@16:       return x + c;
Chris@16:     else
Chris@16:       return x - (m - c);
Chris@16:   }
Chris@16: 
Chris@16:   static IntType mult(IntType a, IntType x)
Chris@16:   {
Chris@16:     if(((unsigned_m() - 1) & unsigned_m()) == 0)
Chris@16:       return unsigned_type(a) * unsigned_type(x) & (unsigned_m() - 1);
Chris@16:     else if(a == 0)
Chris@16:       return 0;
Chris@16:     else if(a == 1)
Chris@16:       return x;
Chris@16:     else if(m <= traits::const_max/a)      // i.e. a*m <= max
Chris@16:       return mult_small(a, x);
Chris@16:     else if(traits::is_signed && (m%a < m/a))
Chris@16:       return mult_schrage(a, x);
Chris@16:     else
Chris@16:       return mult_general(a, x);
Chris@16:   }
Chris@16: 
Chris@16:   static IntType mult_add(IntType a, IntType x, IntType c)
Chris@16:   {
Chris@16:     if(((unsigned_m() - 1) & unsigned_m()) == 0)
Chris@16:       return (unsigned_type(a) * unsigned_type(x) + unsigned_type(c)) & (unsigned_m() - 1);
Chris@16:     else if(a == 0)
Chris@16:       return c;
Chris@16:     else if(m <= (traits::const_max-c)/a) {  // i.e. a*m+c <= max
Chris@101:       IntType suppress_warnings = (m == 0);
Chris@101:       BOOST_ASSERT(suppress_warnings == 0);
Chris@101:       return (a*x+c) % (m + suppress_warnings);
Chris@16:     } else
Chris@16:       return add(mult(a, x), c);
Chris@16:   }
Chris@16: 
Chris@16:   static IntType pow(IntType a, boost::uintmax_t exponent)
Chris@16:   {
Chris@16:       IntType result = 1;
Chris@16:       while(exponent != 0) {
Chris@16:           if(exponent % 2 == 1) {
Chris@16:               result = mult(result, a);
Chris@16:           }
Chris@16:           a = mult(a, a);
Chris@16:           exponent /= 2;
Chris@16:       }
Chris@16:       return result;
Chris@16:   }
Chris@16: 
Chris@16:   static IntType invert(IntType x)
Chris@16:   { return x == 0 ? 0 : (m == 0? invert_euclidian0(x) : invert_euclidian(x)); }
Chris@16: 
Chris@16: private:
Chris@16:   typedef integer_traits<IntType> traits;
Chris@16:   typedef typename make_unsigned<IntType>::type unsigned_type;
Chris@16: 
Chris@16:   const_mod();      // don't instantiate
Chris@16: 
Chris@16:   static IntType mult_small(IntType a, IntType x)
Chris@16:   {
Chris@101:     IntType suppress_warnings = (m == 0);
Chris@101:     BOOST_ASSERT(suppress_warnings == 0);
Chris@101:     return a*x % (m + suppress_warnings);
Chris@16:   }
Chris@16: 
Chris@16:   static IntType mult_schrage(IntType a, IntType value)
Chris@16:   {
Chris@16:     const IntType q = m / a;
Chris@16:     const IntType r = m % a;
Chris@16: 
Chris@16:     BOOST_ASSERT(r < q);        // check that overflow cannot happen
Chris@16: 
Chris@16:     return sub(a*(value%q), r*(value/q));
Chris@16:   }
Chris@16: 
Chris@16:   static IntType mult_general(IntType a, IntType b)
Chris@16:   {
Chris@16:     IntType suppress_warnings = (m == 0);
Chris@16:     BOOST_ASSERT(suppress_warnings == 0);
Chris@16:     IntType modulus = m + suppress_warnings;
Chris@16:     BOOST_ASSERT(modulus == m);
Chris@16:     if(::boost::uintmax_t(modulus) <=
Chris@16:         (::std::numeric_limits< ::boost::uintmax_t>::max)() / modulus)
Chris@16:     {
Chris@16:       return static_cast<IntType>(boost::uintmax_t(a) * b % modulus);
Chris@16:     } else {
Chris@16:       return static_cast<IntType>(detail::mulmod(a, b, modulus));
Chris@16:     }
Chris@16:   }
Chris@16: 
Chris@16:   static IntType sub(IntType a, IntType b)
Chris@16:   {
Chris@16:     if(a < b)
Chris@16:       return m - (b - a);
Chris@16:     else
Chris@16:       return a - b;
Chris@16:   }
Chris@16: 
Chris@16:   static unsigned_type unsigned_m()
Chris@16:   {
Chris@16:       if(m == 0) {
Chris@16:           return unsigned_type((std::numeric_limits<IntType>::max)()) + 1;
Chris@16:       } else {
Chris@16:           return unsigned_type(m);
Chris@16:       }
Chris@16:   }
Chris@16: 
Chris@16:   // invert c in the finite field (mod m) (m must be prime)
Chris@16:   static IntType invert_euclidian(IntType c)
Chris@16:   {
Chris@16:     // we are interested in the gcd factor for c, because this is our inverse
Chris@16:     BOOST_ASSERT(c > 0);
Chris@16:     IntType l1 = 0;
Chris@16:     IntType l2 = 1;
Chris@16:     IntType n = c;
Chris@16:     IntType p = m;
Chris@16:     for(;;) {
Chris@16:       IntType q = p / n;
Chris@16:       l1 += q * l2;
Chris@16:       p -= q * n;
Chris@16:       if(p == 0)
Chris@16:         return l2;
Chris@16:       IntType q2 = n / p;
Chris@16:       l2 += q2 * l1;
Chris@16:       n -= q2 * p;
Chris@16:       if(n == 0)
Chris@16:         return m - l1;
Chris@16:     }
Chris@16:   }
Chris@16: 
Chris@16:   // invert c in the finite field (mod m) (c must be relatively prime to m)
Chris@16:   static IntType invert_euclidian0(IntType c)
Chris@16:   {
Chris@16:     // we are interested in the gcd factor for c, because this is our inverse
Chris@16:     BOOST_ASSERT(c > 0);
Chris@16:     if(c == 1) return 1;
Chris@16:     IntType l1 = 0;
Chris@16:     IntType l2 = 1;
Chris@16:     IntType n = c;
Chris@16:     IntType p = m;
Chris@16:     IntType max = (std::numeric_limits<IntType>::max)();
Chris@16:     IntType q = max / n;
Chris@16:     BOOST_ASSERT(max % n != n - 1 && "c must be relatively prime to m.");
Chris@16:     l1 += q * l2;
Chris@16:     p = max - q * n + 1;
Chris@16:     for(;;) {
Chris@16:       if(p == 0)
Chris@16:         return l2;
Chris@16:       IntType q2 = n / p;
Chris@16:       l2 += q2 * l1;
Chris@16:       n -= q2 * p;
Chris@16:       if(n == 0)
Chris@16:         return m - l1;
Chris@16:       q = p / n;
Chris@16:       l1 += q * l2;
Chris@16:       p -= q * n;
Chris@16:     }
Chris@16:   }
Chris@16: };
Chris@16: 
Chris@16: } // namespace random
Chris@16: } // namespace boost
Chris@16: 
Chris@16: #include <boost/random/detail/enable_warnings.hpp>
Chris@16: 
Chris@16: #endif // BOOST_RANDOM_CONST_MOD_HPP