williamr@2: /* boost random/detail/const_mod.hpp header file
williamr@2:  *
williamr@2:  * Copyright Jens Maurer 2000-2001
williamr@2:  * Distributed under the Boost Software License, Version 1.0. (See
williamr@2:  * accompanying file LICENSE_1_0.txt or copy at
williamr@2:  * http://www.boost.org/LICENSE_1_0.txt)
williamr@2:  *
williamr@2:  * See http://www.boost.org for most recent version including documentation.
williamr@2:  *
williamr@2:  * $Id: const_mod.hpp,v 1.8 2004/07/27 03:43:32 dgregor Exp $
williamr@2:  *
williamr@2:  * Revision history
williamr@2:  *  2001-02-18  moved to individual header files
williamr@2:  */
williamr@2: 
williamr@2: #ifndef BOOST_RANDOM_CONST_MOD_HPP
williamr@2: #define BOOST_RANDOM_CONST_MOD_HPP
williamr@2: 
williamr@2: #include <cassert>
williamr@2: #include <boost/static_assert.hpp>
williamr@2: #include <boost/cstdint.hpp>
williamr@2: #include <boost/integer_traits.hpp>
williamr@2: #include <boost/detail/workaround.hpp>
williamr@2: 
williamr@2: namespace boost {
williamr@2: namespace random {
williamr@2: 
williamr@2: /*
williamr@2:  * Some random number generators require modular arithmetic.  Put
williamr@2:  * everything we need here.
williamr@2:  * IntType must be an integral type.
williamr@2:  */
williamr@2: 
williamr@2: namespace detail {
williamr@2: 
williamr@2:   template<bool is_signed>
williamr@2:   struct do_add
williamr@2:   { };
williamr@2: 
williamr@2:   template<>
williamr@2:   struct do_add<true>
williamr@2:   {
williamr@2:     template<class IntType>
williamr@2:     static IntType add(IntType m, IntType x, IntType c)
williamr@2:     {
williamr@2:       x += (c-m);
williamr@2:       if(x < 0)
williamr@2:         x += m;
williamr@2:       return x;
williamr@2:     }
williamr@2:   };
williamr@2: 
williamr@2:   template<>
williamr@2:   struct do_add<false>
williamr@2:   {
williamr@2:     template<class IntType>
williamr@2:     static IntType add(IntType, IntType, IntType)
williamr@2:     {
williamr@2:       // difficult
williamr@2:       assert(!"const_mod::add with c too large");
williamr@2:       return 0;
williamr@2:     }
williamr@2:   };
williamr@2: } // namespace detail
williamr@2: 
williamr@2: #if !(defined(__BORLANDC__) && (__BORLANDC__ == 0x560))
williamr@2: 
williamr@2: template<class IntType, IntType m>
williamr@2: class const_mod
williamr@2: {
williamr@2: public:
williamr@2:   static IntType add(IntType x, IntType c)
williamr@2:   {
williamr@2:     if(c == 0)
williamr@2:       return x;
williamr@2:     else if(c <= traits::const_max - m)    // i.e. m+c < max
williamr@2:       return add_small(x, c);
williamr@2:     else
williamr@2:       return detail::do_add<traits::is_signed>::add(m, x, c);
williamr@2:   }
williamr@2: 
williamr@2:   static IntType mult(IntType a, IntType x)
williamr@2:   {
williamr@2:     if(a == 1)
williamr@2:       return x;
williamr@2:     else if(m <= traits::const_max/a)      // i.e. a*m <= max
williamr@2:       return mult_small(a, x);
williamr@2:     else if(traits::is_signed && (m%a < m/a))
williamr@2:       return mult_schrage(a, x);
williamr@2:     else {
williamr@2:       // difficult
williamr@2:       assert(!"const_mod::mult with a too large");
williamr@2:       return 0;
williamr@2:     }
williamr@2:   }
williamr@2: 
williamr@2:   static IntType mult_add(IntType a, IntType x, IntType c)
williamr@2:   {
williamr@2:     if(m <= (traits::const_max-c)/a)   // i.e. a*m+c <= max
williamr@2:       return (a*x+c) % m;
williamr@2:     else
williamr@2:       return add(mult(a, x), c);
williamr@2:   }
williamr@2: 
williamr@2:   static IntType invert(IntType x)
williamr@2:   { return x == 0 ? 0 : invert_euclidian(x); }
williamr@2: 
williamr@2: private:
williamr@2:   typedef integer_traits<IntType> traits;
williamr@2: 
williamr@2:   const_mod();      // don't instantiate
williamr@2: 
williamr@2:   static IntType add_small(IntType x, IntType c)
williamr@2:   {
williamr@2:     x += c;
williamr@2:     if(x >= m)
williamr@2:       x -= m;
williamr@2:     return x;
williamr@2:   }
williamr@2: 
williamr@2:   static IntType mult_small(IntType a, IntType x)
williamr@2:   {
williamr@2:     return a*x % m;
williamr@2:   }
williamr@2: 
williamr@2:   static IntType mult_schrage(IntType a, IntType value)
williamr@2:   {
williamr@2:     const IntType q = m / a;
williamr@2:     const IntType r = m % a;
williamr@2: 
williamr@2:     assert(r < q);        // check that overflow cannot happen
williamr@2: 
williamr@2:     value = a*(value%q) - r*(value/q);
williamr@2:     // An optimizer bug in the SGI MIPSpro 7.3.1.x compiler requires this
williamr@2:     // convoluted formulation of the loop (Synge Todo)
williamr@2:     for(;;) {
williamr@2:       if (value > 0)
williamr@2:         break;
williamr@2:       value += m;
williamr@2:     }
williamr@2:     return value;
williamr@2:   }
williamr@2: 
williamr@2:   // invert c in the finite field (mod m) (m must be prime)
williamr@2:   static IntType invert_euclidian(IntType c)
williamr@2:   {
williamr@2:     // we are interested in the gcd factor for c, because this is our inverse
williamr@2:     BOOST_STATIC_ASSERT(m > 0);
williamr@2: #if BOOST_WORKAROUND(__MWERKS__, BOOST_TESTED_AT(0x3003))
williamr@2:     assert(boost::integer_traits<IntType>::is_signed);
williamr@2: #elif !defined(BOOST_NO_LIMITS_COMPILE_TIME_CONSTANTS)
williamr@2:     BOOST_STATIC_ASSERT(boost::integer_traits<IntType>::is_signed);
williamr@2: #endif
williamr@2:     assert(c > 0);
williamr@2:     IntType l1 = 0;
williamr@2:     IntType l2 = 1;
williamr@2:     IntType n = c;
williamr@2:     IntType p = m;
williamr@2:     for(;;) {
williamr@2:       IntType q = p / n;
williamr@2:       l1 -= q * l2;           // this requires a signed IntType!
williamr@2:       p -= q * n;
williamr@2:       if(p == 0)
williamr@2:         return (l2 < 1 ? l2 + m : l2);
williamr@2:       IntType q2 = n / p;
williamr@2:       l2 -= q2 * l1;
williamr@2:       n -= q2 * p;
williamr@2:       if(n == 0)
williamr@2:         return (l1 < 1 ? l1 + m : l1);
williamr@2:     }
williamr@2:   }
williamr@2: };
williamr@2: 
williamr@2: // The modulus is exactly the word size: rely on machine overflow handling.
williamr@2: // Due to a GCC bug, we cannot partially specialize in the presence of
williamr@2: // template value parameters.
williamr@2: template<>
williamr@2: class const_mod<unsigned int, 0>
williamr@2: {
williamr@2:   typedef unsigned int IntType;
williamr@2: public:
williamr@2:   static IntType add(IntType x, IntType c) { return x+c; }
williamr@2:   static IntType mult(IntType a, IntType x) { return a*x; }
williamr@2:   static IntType mult_add(IntType a, IntType x, IntType c) { return a*x+c; }
williamr@2: 
williamr@2:   // m is not prime, thus invert is not useful
williamr@2: private:                      // don't instantiate
williamr@2:   const_mod();
williamr@2: };
williamr@2: 
williamr@2: template<>
williamr@2: class const_mod<unsigned long, 0>
williamr@2: {
williamr@2:   typedef unsigned long IntType;
williamr@2: public:
williamr@2:   static IntType add(IntType x, IntType c) { return x+c; }
williamr@2:   static IntType mult(IntType a, IntType x) { return a*x; }
williamr@2:   static IntType mult_add(IntType a, IntType x, IntType c) { return a*x+c; }
williamr@2: 
williamr@2:   // m is not prime, thus invert is not useful
williamr@2: private:                      // don't instantiate
williamr@2:   const_mod();
williamr@2: };
williamr@2: 
williamr@2: // the modulus is some power of 2: rely partly on machine overflow handling
williamr@2: // we only specialize for rand48 at the moment
williamr@2: #ifndef BOOST_NO_INT64_T
williamr@2: template<>
williamr@2: class const_mod<uint64_t, uint64_t(1) << 48>
williamr@2: {
williamr@2:   typedef uint64_t IntType;
williamr@2: public:
williamr@2:   static IntType add(IntType x, IntType c) { return c == 0 ? x : mod(x+c); }
williamr@2:   static IntType mult(IntType a, IntType x) { return mod(a*x); }
williamr@2:   static IntType mult_add(IntType a, IntType x, IntType c)
williamr@2:     { return mod(a*x+c); }
williamr@2:   static IntType mod(IntType x) { return x &= ((uint64_t(1) << 48)-1); }
williamr@2: 
williamr@2:   // m is not prime, thus invert is not useful
williamr@2: private:                      // don't instantiate
williamr@2:   const_mod();
williamr@2: };
williamr@2: #endif /* !BOOST_NO_INT64_T */
williamr@2: 
williamr@2: #else
williamr@2: 
williamr@2: //
williamr@2: // for some reason Borland C++ Builder 6 has problems with
williamr@2: // the full specialisations of const_mod, define a generic version
williamr@2: // instead, the compiler will optimise away the const-if statements:
williamr@2: //
williamr@2: 
williamr@2: template<class IntType, IntType m>
williamr@2: class const_mod
williamr@2: {
williamr@2: public:
williamr@2:   static IntType add(IntType x, IntType c)
williamr@2:   {
williamr@2:     if(0 == m)
williamr@2:     {
williamr@2:        return x+c;
williamr@2:     }
williamr@2:     else
williamr@2:     {
williamr@2:        if(c == 0)
williamr@2:          return x;
williamr@2:        else if(c <= traits::const_max - m)    // i.e. m+c < max
williamr@2:          return add_small(x, c);
williamr@2:        else
williamr@2:          return detail::do_add<traits::is_signed>::add(m, x, c);
williamr@2:     }
williamr@2:   }
williamr@2: 
williamr@2:   static IntType mult(IntType a, IntType x)
williamr@2:   {
williamr@2:     if(x == 0)
williamr@2:     {
williamr@2:        return a*x;
williamr@2:     }
williamr@2:     else
williamr@2:     {
williamr@2:        if(a == 1)
williamr@2:          return x;
williamr@2:        else if(m <= traits::const_max/a)      // i.e. a*m <= max
williamr@2:          return mult_small(a, x);
williamr@2:        else if(traits::is_signed && (m%a < m/a))
williamr@2:          return mult_schrage(a, x);
williamr@2:        else {
williamr@2:          // difficult
williamr@2:          assert(!"const_mod::mult with a too large");
williamr@2:          return 0;
williamr@2:        }
williamr@2:     }
williamr@2:   }
williamr@2: 
williamr@2:   static IntType mult_add(IntType a, IntType x, IntType c)
williamr@2:   {
williamr@2:     if(m == 0)
williamr@2:     {
williamr@2:        return a*x+c;
williamr@2:     }
williamr@2:     else
williamr@2:     {
williamr@2:        if(m <= (traits::const_max-c)/a)   // i.e. a*m+c <= max
williamr@2:          return (a*x+c) % m;
williamr@2:        else
williamr@2:          return add(mult(a, x), c);
williamr@2:     }
williamr@2:   }
williamr@2: 
williamr@2:   static IntType invert(IntType x)
williamr@2:   { return x == 0 ? 0 : invert_euclidian(x); }
williamr@2: 
williamr@2: private:
williamr@2:   typedef integer_traits<IntType> traits;
williamr@2: 
williamr@2:   const_mod();      // don't instantiate
williamr@2: 
williamr@2:   static IntType add_small(IntType x, IntType c)
williamr@2:   {
williamr@2:     x += c;
williamr@2:     if(x >= m)
williamr@2:       x -= m;
williamr@2:     return x;
williamr@2:   }
williamr@2: 
williamr@2:   static IntType mult_small(IntType a, IntType x)
williamr@2:   {
williamr@2:     return a*x % m;
williamr@2:   }
williamr@2: 
williamr@2:   static IntType mult_schrage(IntType a, IntType value)
williamr@2:   {
williamr@2:     const IntType q = m / a;
williamr@2:     const IntType r = m % a;
williamr@2: 
williamr@2:     assert(r < q);        // check that overflow cannot happen
williamr@2: 
williamr@2:     value = a*(value%q) - r*(value/q);
williamr@2:     while(value <= 0)
williamr@2:       value += m;
williamr@2:     return value;
williamr@2:   }
williamr@2: 
williamr@2:   // invert c in the finite field (mod m) (m must be prime)
williamr@2:   static IntType invert_euclidian(IntType c)
williamr@2:   {
williamr@2:     // we are interested in the gcd factor for c, because this is our inverse
williamr@2:     BOOST_STATIC_ASSERT(m > 0);
williamr@2: #ifndef BOOST_NO_LIMITS_COMPILE_TIME_CONSTANTS
williamr@2:     BOOST_STATIC_ASSERT(boost::integer_traits<IntType>::is_signed);
williamr@2: #endif
williamr@2:     assert(c > 0);
williamr@2:     IntType l1 = 0;
williamr@2:     IntType l2 = 1;
williamr@2:     IntType n = c;
williamr@2:     IntType p = m;
williamr@2:     for(;;) {
williamr@2:       IntType q = p / n;
williamr@2:       l1 -= q * l2;           // this requires a signed IntType!
williamr@2:       p -= q * n;
williamr@2:       if(p == 0)
williamr@2:         return (l2 < 1 ? l2 + m : l2);
williamr@2:       IntType q2 = n / p;
williamr@2:       l2 -= q2 * l1;
williamr@2:       n -= q2 * p;
williamr@2:       if(n == 0)
williamr@2:         return (l1 < 1 ? l1 + m : l1);
williamr@2:     }
williamr@2:   }
williamr@2: };
williamr@2: 
williamr@2: 
williamr@2: #endif
williamr@2: 
williamr@2: } // namespace random
williamr@2: } // namespace boost
williamr@2: 
williamr@2: #endif // BOOST_RANDOM_CONST_MOD_HPP