CODEWITHSUNDEEP
c++templatesgccc++11random
thor
thor

Reputation: 22510

c++11 random number generation: how to re-implement `uniform_int_distribution` with `mt19937` without templates

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();
      }

Upvotes: 4

Views: 2473

Answers (1)

Drake
Drake

Reputation: 857

  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.

Upvotes: 3

Related Questions