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