[Boost-users] Random permutations using boost::random_number_generator and std::random_shuffle.

19 Jan 2011

      Hi,

We need the ability to randomly permute some range, and we already
have this functionality available in the computer algebra system
Maxima. We use Maxima for numerous things, often implementing basic
algorithms in it, then later writing more efficient versions in C++.

The problem we have is that we'd like to be able to make sure that any
function we write to randomly shuffle a range, matches up with the
Maxima "random_permutation" function.

We believe it uses the Mersenne twister (boost::mt19937) and Knuth
shuffle algorithms, and so we thought to implement these in C++ using
boost::uniform_distribution, boost:variate_generator etc.

However, we have two questions (open problems):

1) How does boost::variate_generator interact with mt19937, and the
distribution? Is this documented somewhere?

In the attached code, it seems we can construct
boost::uniform_distribution with any maximum value, and still the
number generator returns values outside of that range? Is this really
right?

2) What algorithm does std::random_shuffle implement? (Not boost
related specifically but hopefully someone will know)

Is this defined by the C++ standard? It doesn't seem to be the simple
Knuth shuffle algorithm (starting from top of range to bottom or vice
versa) as we implement in the attached code. If anyone could point me
in the right direction, I'd very much appreciate it.

Thanks

Matthew Gwynne
http://cs.swan.ac.uk/~csmg/

Code:

#include <vector>
#include <iostream>
#include <algorithm>
#include <iterator>
#include <cassert>

#include <boost/lexical_cast.hpp>
#include <boost/random/mersenne_twister.hpp>
#include <boost/random/uniform_int.hpp>
#include <boost/random/variate_generator.hpp>
#include <boost/random/random_number_generator.hpp>

namespace {

  enum { errcode_parameter = 1 };

  const std::string program = "RandomShuffle";
  const std::string err = "ERROR[" + program + "]: ";

  const std::string version = "0.0.3";

  const int default_seed = 1;

  typedef boost::mt19937 base_generator_type;
  base_generator_type base_rand_gen;
  inline void set_random(const int seed) {
    assert(seed >= 1);
    base_rand_gen.seed(seed);
  }

  const unsigned int default_N = 10;

  typedef std::vector<int> vec_t;
  inline void initialise(vec_t& a) {
    for (unsigned int i = 0; i < a.size(); ++i) a[i] = i+1;
  }
  inline void output(const vec_t& a) {
    for (unsigned int i = 0; i < a.size(); ++i) std::cout << a[i] << " ";
    std::cout << "\n";
  }

  template <class RandomAccessIterator, class RandomNumberGenerator>
  inline void random_shuffle(const RandomAccessIterator first, const
RandomAccessIterator last,
    RandomNumberGenerator& rand) {
    typedef typename
std::iterator_traits<RandomAccessIterator>::difference_type
difference_type;
    difference_type n = last - first;
    assert(n >= 0);
    if (n <= 1) return;
    for (RandomAccessIterator i = first; i != last-1; ++i,--n) {
      const difference_type r = rand(n);
      assert(r >= 0);
      assert(r < n);
      std::iter_swap(i, i+r);
    }
    assert(n == 1);
  }

}

int main(const int argc, const char* const argv[]) {
  if (argc > 3) {
    std::cerr << err << "\n At most two arguments are allowed "
      "(the seed for the random-number generator and the number of
elements).\n";
    return errcode_parameter;
  }

  const int seed = (argc == 1) ? default_seed :
boost::lexical_cast<int>(argv[1]);
  const unsigned int N = (argc < 3) ? default_N :
boost::lexical_cast<unsigned int>(argv[2]);

  vec_t V(N);

  typedef boost::uniform_int<> uniform_distribution_type;
  uniform_distribution_type
uniform_distribution(0,std::numeric_limits<int>::max()); // is this
correct???
  typedef boost::variate_generator<base_generator_type&,
uniform_distribution_type> generator_type;
  generator_type rand_gen(base_rand_gen, uniform_distribution);
  typedef boost::random_number_generator<generator_type> RandomNumberGenerator;
  RandomNumberGenerator rg(rand_gen);

  set_random(seed);
  std::cout << "Underlying random sequence:\n";
  for (long n = V.size(); n > 1; --n) std::cout << rg(n) << " ";
  std::cout << "\n\n";

  initialise(V);
  set_random(seed);
  std::random_shuffle(V.begin(), V.end(), rg);
  std::cout << "std::random_shuffle:\n";
  output(V);

  initialise(V);
  set_random(seed);
  ::random_shuffle(V.begin(), V.end(), rg);
  std::cout << "::random_shuffle:\n";
  output(V);
}