4

I am trying to reimplement c++11 uniform_int_distribution without templates and specifically for mt19937, because I wanted to port functionality in this typical use case to other languages without template facility, and it would be nice to have a more readable version. The difficulty from reading gcc implementation far exceeds what I expected from this mathematically simple conversion from one uniform distribution to another. (Yes, I know this is specializing the general, and I wouldn't gain new functionality from this practice)

I looked into gcc 4.8.1 headers. The most used mt19937 class is a typedef:

  typedef mersenne_twister_engine<
    uint_fast32_t,
    32, 624, 397, 31,
    0x9908b0dfUL, 11,
    0xffffffffUL, 7,
    0x9d2c5680UL, 15,
    0xefc60000UL, 18, 1812433253UL> mt19937;

The code for uniform_int_distribution in gcc is heavily templated, and not very readable to me. I wonder how I can simplify/specialize that code to non-template code, just for the mt19937 case.

The most relevant piece of code I found from 4.8.1/include/c++/bits/random.tcc is attached in the end (with double underscore __ removed for clarity).

I tried to specialize the code but wasn't very successful. For starter, I tried to figure out the range of mt19937: the min is 0, the max is (in random.h):

  static constexpr result_type
  max()
  { return __detail::_Shift<_UIntType, __w>::__value - 1; }

, which involves complicated template programming that's not easily readable. I figured maybe it's better to ask than reverse engineer the templates.

So my questions are:

  1. What' the range (max) for the specific type mt19937?
  2. How to specialize the rest of the code with types and values specific to uint32 and mt19973.

Thanks in advance.

--- gcc 4.8.1 code for sampling from mt19937 to uniform_int_distribution ---

  template<typename _IntType>
    template<typename _ForwardIterator,
         typename _UniformRandomNumberGenerator>
      void
      uniform_int_distribution<_IntType>::
      generate_impl(_ForwardIterator f, _ForwardIterator t,
              _UniformRandomNumberGenerator& urng,
              const param_type& param)
      {
    glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>)
    typedef typename _UniformRandomNumberGenerator::result_type
      _Gresult_type;
    typedef typename std::make_unsigned<result_type>::type utype;
    typedef typename std::common_type<_Gresult_type, utype>::type
      uctype;

    const uctype urngmin = urng.min();
    const uctype urngmax = urng.max();
    const uctype urngrange = urngmax - urngmin;
    const uctype urange
      = uctype(param.b()) - uctype(param.a());

    uctype ret;

    if (urngrange > urange)
      {
        if (detail::_Power_of_2(urngrange + 1)
        && detail::_Power_of_2(urange + 1))
          {
        while (f != t)
          {
            ret = uctype(urng()) - urngmin;
            *f++ = (ret & urange) + param.a();
          }
          }
        else
          {
        // downscaling
        const uctype uerange = urange + 1; // urange can be zero
        const uctype scaling = urngrange / uerange;
        const uctype past = uerange * scaling;
        while (f != t)
          {
            do
              ret = uctype(urng()) - urngmin;
            while (ret >= past);
            *f++ = ret / scaling + param.a();
          }
          }
      }
    else if (urngrange < urange)
      {
        // upscaling
        /*
          Note that every value in [0, urange]
          can be written uniquely as

          (urngrange + 1) * high + low

          where

          high in [0, urange / (urngrange + 1)]

          and

          low in [0, urngrange].
        */
        uctype tmp; // wraparound control
        while (f != t)
          {
        do
          {
            const uctype uerngrange = urngrange + 1;
            tmp = (uerngrange * operator()
                 (urng, param_type(0, urange / uerngrange)));
            ret = tmp + (uctype(urng()) - urngmin);
          }
        while (ret > urange || ret < tmp);
        *f++ = ret;
          }
      }
    else
      while (f != t)
        *f++ = uctype(urng()) - urngmin + param.a();
      }
thor
  • 21,418
  • 31
  • 87
  • 173

1 Answers1

3
  1. mt19937 generates integers in [0,2^32-1]:

    std::mt19937 mt_gen;
    std::cout << mt_gen.min() << '\n';
    std::cout << mt_gen.max() << '\n';
    

gives

    0
    4294967295

2. If I understand correctly, you want to specialise the template "by hand", but without figuring out exactly what uniform_int_distribution does? Once you have mt19937 (e.g., http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/MT2002/CODES/mt19937ar.c), generating uniformly distributed integers in a given range [low, high] is conceptually simple, but there are a few details that require attention (for example, meticulously checking for off-by-one errors). The second answer here Generating a uniform distribution of INTEGERS in C (after the accepted one) might be helpful.

Community
  • 1
  • 1
Drake
  • 857
  • 5
  • 10
  • +1 Thanks. It's straight forward to figure out the range this way. – thor Mar 20 '14 at 11:35
  • Is there actual code/pseudo-code for the process describe therein? – thor Mar 20 '14 at 14:29
  • I am using C++11 so I don't have a link handy that would be suitable for your needs. If you have access to the book "Numerical Recipes in C", then you may want to look at the 7th chapter. Especially regarding the extraction of bits from already generated random numbers, and how random the bits can be expected to be. – Drake Mar 21 '14 at 07:38