/* 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.8 2004/07/27 03:43:32 dgregor 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)

    {

      x += (c-m);

      if(x < 0)

        x += m;

      return x;

    }

  };

  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