boost/compute/random/mersenne_twister_engine.hpp
//---------------------------------------------------------------------------//
// Copyright (c) 2013 Kyle Lutz <kyle.r.lutz@gmail.com>
//
// 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://boostorg.github.com/compute for more information.
//---------------------------------------------------------------------------//
#ifndef BOOST_COMPUTE_RANDOM_MERSENNE_TWISTER_ENGINE_HPP
#define BOOST_COMPUTE_RANDOM_MERSENNE_TWISTER_ENGINE_HPP
#include <algorithm>
#include <boost/compute/types.hpp>
#include <boost/compute/buffer.hpp>
#include <boost/compute/kernel.hpp>
#include <boost/compute/context.hpp>
#include <boost/compute/program.hpp>
#include <boost/compute/command_queue.hpp>
#include <boost/compute/algorithm/transform.hpp>
#include <boost/compute/container/vector.hpp>
#include <boost/compute/detail/iterator_range_size.hpp>
#include <boost/compute/iterator/discard_iterator.hpp>
#include <boost/compute/utility/program_cache.hpp>
namespace boost {
namespace compute {
/// \class mersenne_twister_engine
/// \brief Mersenne twister pseudorandom number generator.
template<class T>
class mersenne_twister_engine
{
public:
typedef T result_type;
static const T default_seed = 5489U;
static const T n = 624;
static const T m = 397;
/// Creates a new mersenne_twister_engine and seeds it with \p value.
explicit mersenne_twister_engine(command_queue &queue,
result_type value = default_seed)
: m_context(queue.get_context()),
m_state_buffer(m_context, n * sizeof(result_type))
{
// setup program
load_program();
// seed state
seed(value, queue);
}
/// Creates a new mersenne_twister_engine object as a copy of \p other.
mersenne_twister_engine(const mersenne_twister_engine<T> &other)
: m_context(other.m_context),
m_state_index(other.m_state_index),
m_program(other.m_program),
m_state_buffer(other.m_state_buffer)
{
}
/// Copies \p other to \c *this.
mersenne_twister_engine<T>& operator=(const mersenne_twister_engine<T> &other)
{
if(this != &other){
m_context = other.m_context;
m_state_index = other.m_state_index;
m_program = other.m_program;
m_state_buffer = other.m_state_buffer;
}
return *this;
}
/// Destroys the mersenne_twister_engine object.
~mersenne_twister_engine()
{
}
/// Seeds the random number generator with \p value.
///
/// \param value seed value for the random-number generator
/// \param queue command queue to perform the operation
///
/// If no seed value is provided, \c default_seed is used.
void seed(result_type value, command_queue &queue)
{
kernel seed_kernel = m_program.create_kernel("seed");
seed_kernel.set_arg(0, value);
seed_kernel.set_arg(1, m_state_buffer);
queue.enqueue_task(seed_kernel);
m_state_index = 0;
}
/// \overload
void seed(command_queue &queue)
{
seed(default_seed, queue);
}
/// Generates random numbers and stores them to the range [\p first, \p last).
template<class OutputIterator>
void generate(OutputIterator first, OutputIterator last, command_queue &queue)
{
const size_t size = detail::iterator_range_size(first, last);
kernel fill_kernel(m_program, "fill");
fill_kernel.set_arg(0, m_state_buffer);
fill_kernel.set_arg(2, first.get_buffer());
size_t offset = 0;
size_t &p = m_state_index;
for(;;){
size_t count = 0;
if(size > n){
count = (std::min)(static_cast<size_t>(n), size - offset);
}
else {
count = size;
}
fill_kernel.set_arg(1, static_cast<const uint_>(p));
fill_kernel.set_arg(3, static_cast<const uint_>(offset));
queue.enqueue_1d_range_kernel(fill_kernel, 0, count, 0);
p += count;
offset += count;
if(offset >= size){
break;
}
generate_state(queue);
p = 0;
}
}
/// \internal_
void generate(discard_iterator first, discard_iterator last, command_queue &queue)
{
(void) queue;
m_state_index += std::distance(first, last);
}
/// Generates random numbers, transforms them with \p op, and then stores
/// them to the range [\p first, \p last).
template<class OutputIterator, class Function>
void generate(OutputIterator first, OutputIterator last, Function op, command_queue &queue)
{
vector<T> tmp(std::distance(first, last), queue.get_context());
generate(tmp.begin(), tmp.end(), queue);
transform(tmp.begin(), tmp.end(), first, op, queue);
}
/// Generates \p z random numbers and discards them.
void discard(size_t z, command_queue &queue)
{
generate(discard_iterator(0), discard_iterator(z), queue);
}
/// \internal_ (deprecated)
template<class OutputIterator>
void fill(OutputIterator first, OutputIterator last, command_queue &queue)
{
generate(first, last, queue);
}
private:
/// \internal_
void generate_state(command_queue &queue)
{
kernel generate_state_kernel =
m_program.create_kernel("generate_state");
generate_state_kernel.set_arg(0, m_state_buffer);
queue.enqueue_task(generate_state_kernel);
}
/// \internal_
void load_program()
{
boost::shared_ptr<program_cache> cache =
program_cache::get_global_cache(m_context);
std::string cache_key =
std::string("__boost_mersenne_twister_engine_") + type_name<T>();
const char source[] =
"static uint twiddle(uint u, uint v)\n"
"{\n"
" return (((u & 0x80000000U) | (v & 0x7FFFFFFFU)) >> 1) ^\n"
" ((v & 1U) ? 0x9908B0DFU : 0x0U);\n"
"}\n"
"__kernel void generate_state(__global uint *state)\n"
"{\n"
" const uint n = 624;\n"
" const uint m = 397;\n"
" for(uint i = 0; i < (n - m); i++)\n"
" state[i] = state[i+m] ^ twiddle(state[i], state[i+1]);\n"
" for(uint i = n - m; i < (n - 1); i++)\n"
" state[i] = state[i+m-n] ^ twiddle(state[i], state[i+1]);\n"
" state[n-1] = state[m-1] ^ twiddle(state[n-1], state[0]);\n"
"}\n"
"__kernel void seed(const uint s, __global uint *state)\n"
"{\n"
" const uint n = 624;\n"
" state[0] = s & 0xFFFFFFFFU;\n"
" for(uint i = 1; i < n; i++){\n"
" state[i] = 1812433253U * (state[i-1] ^ (state[i-1] >> 30)) + i;\n"
" state[i] &= 0xFFFFFFFFU;\n"
" }\n"
" generate_state(state);\n"
"}\n"
"static uint random_number(__global uint *state, const uint p)\n"
"{\n"
" uint x = state[p];\n"
" x ^= (x >> 11);\n"
" x ^= (x << 7) & 0x9D2C5680U;\n"
" x ^= (x << 15) & 0xEFC60000U;\n"
" return x ^ (x >> 18);\n"
"}\n"
"__kernel void fill(__global uint *state,\n"
" const uint state_index,\n"
" __global uint *vector,\n"
" const uint offset)\n"
"{\n"
" const uint i = get_global_id(0);\n"
" vector[offset+i] = random_number(state, state_index + i);\n"
"}\n";
m_program = cache->get_or_build(cache_key, std::string(), source, m_context);
}
private:
context m_context;
size_t m_state_index;
program m_program;
buffer m_state_buffer;
};
typedef mersenne_twister_engine<uint_> mt19937;
} // end compute namespace
} // end boost namespace
#endif // BOOST_COMPUTE_RANDOM_MERSENNE_TWISTER_ENGINE_HPP