ceph/src/boost/boost/compute/random/mersenne_twister_engine.hpp

   1 //---------------------------------------------------------------------------//
   2 // Copyright (c) 2013 Kyle Lutz <kyle.r.lutz@gmail.com>
   3 //
   4 // Distributed under the Boost Software License, Version 1.0
   5 // See accompanying file LICENSE_1_0.txt or copy at
   6 // http://www.boost.org/LICENSE_1_0.txt
   7 //
   8 // See http://boostorg.github.com/compute for more information.
   9 //---------------------------------------------------------------------------//
  10
  11 #ifndef BOOST_COMPUTE_RANDOM_MERSENNE_TWISTER_ENGINE_HPP
  12 #define BOOST_COMPUTE_RANDOM_MERSENNE_TWISTER_ENGINE_HPP
  13
  14 #include <algorithm>
  15
  16 #include <boost/compute/types.hpp>
  17 #include <boost/compute/buffer.hpp>
  18 #include <boost/compute/kernel.hpp>
  19 #include <boost/compute/context.hpp>
  20 #include <boost/compute/program.hpp>
  21 #include <boost/compute/command_queue.hpp>
  22 #include <boost/compute/algorithm/transform.hpp>
  23 #include <boost/compute/container/vector.hpp>
  24 #include <boost/compute/detail/iterator_range_size.hpp>
  25 #include <boost/compute/iterator/discard_iterator.hpp>
  26 #include <boost/compute/utility/program_cache.hpp>
  27
  28 namespace boost {
  29 namespace compute {
  30
  31 /// \class mersenne_twister_engine
  32 /// \brief Mersenne twister pseudorandom number generator.
  33 template<class T>
  34 class mersenne_twister_engine
  35 {
  36 public:
  37     typedef T result_type;
  38     static const T default_seed = 5489U;
  39     static const T n = 624;
  40     static const T m = 397;
  41
  42     /// Creates a new mersenne_twister_engine and seeds it with \p value.
  43     explicit mersenne_twister_engine(command_queue &queue,
  44                                      result_type value = default_seed)
  45         : m_context(queue.get_context()),
  46           m_state_buffer(m_context, n * sizeof(result_type))
  47     {
  48         // setup program
  49         load_program();
  50
  51         // seed state
  52         seed(value, queue);
  53     }
  54
  55     /// Creates a new mersenne_twister_engine object as a copy of \p other.
  56     mersenne_twister_engine(const mersenne_twister_engine<T> &other)
  57         : m_context(other.m_context),
  58           m_state_index(other.m_state_index),
  59           m_program(other.m_program),
  60           m_state_buffer(other.m_state_buffer)
  61     {
  62     }
  63
  64     /// Copies \p other to \c *this.
  65     mersenne_twister_engine<T>& operator=(const mersenne_twister_engine<T> &other)
  66     {
  67         if(this != &other){
  68             m_context = other.m_context;
  69             m_state_index = other.m_state_index;
  70             m_program = other.m_program;
  71             m_state_buffer = other.m_state_buffer;
  72         }
  73
  74         return *this;
  75     }
  76
  77     /// Destroys the mersenne_twister_engine object.
  78     ~mersenne_twister_engine()
  79     {
  80     }
  81
  82     /// Seeds the random number generator with \p value.
  83     ///
  84     /// \param value seed value for the random-number generator
  85     /// \param queue command queue to perform the operation
  86     ///
  87     /// If no seed value is provided, \c default_seed is used.
  88     void seed(result_type value, command_queue &queue)
  89     {
  90         kernel seed_kernel = m_program.create_kernel("seed");
  91         seed_kernel.set_arg(0, value);
  92         seed_kernel.set_arg(1, m_state_buffer);
  93
  94         queue.enqueue_task(seed_kernel);
  95
  96         m_state_index = 0;
  97     }
  98
  99     /// \overload
 100     void seed(command_queue &queue)
 101     {
 102         seed(default_seed, queue);
 103     }
 104
 105     /// Generates random numbers and stores them to the range [\p first, \p last).
 106     template<class OutputIterator>
 107     void generate(OutputIterator first, OutputIterator last, command_queue &queue)
 108     {
 109         const size_t size = detail::iterator_range_size(first, last);
 110
 111         kernel fill_kernel(m_program, "fill");
 112         fill_kernel.set_arg(0, m_state_buffer);
 113         fill_kernel.set_arg(2, first.get_buffer());
 114
 115         size_t offset = 0;
 116         size_t &p = m_state_index;
 117
 118         for(;;){
 119             size_t count = 0;
 120             if(size > n){
 121                 count = (std::min)(static_cast<size_t>(n), size - offset);
 122             }
 123             else {
 124                 count = size;
 125             }
 126             fill_kernel.set_arg(1, static_cast<const uint_>(p));
 127             fill_kernel.set_arg(3, static_cast<const uint_>(offset));
 128             queue.enqueue_1d_range_kernel(fill_kernel, 0, count, 0);
 129
 130             p += count;
 131             offset += count;
 132
 133             if(offset >= size){
 134                 break;
 135             }
 136
 137             generate_state(queue);
 138             p = 0;
 139         }
 140     }
 141
 142     /// \internal_
 143     void generate(discard_iterator first, discard_iterator last, command_queue &queue)
 144     {
 145         (void) queue;
 146
 147         m_state_index += std::distance(first, last);
 148     }
 149
 150     /// Generates random numbers, transforms them with \p op, and then stores
 151     /// them to the range [\p first, \p last).
 152     template<class OutputIterator, class Function>
 153     void generate(OutputIterator first, OutputIterator last, Function op, command_queue &queue)
 154     {
 155         vector<T> tmp(std::distance(first, last), queue.get_context());
 156         generate(tmp.begin(), tmp.end(), queue);
 157         transform(tmp.begin(), tmp.end(), first, op, queue);
 158     }
 159
 160     /// Generates \p z random numbers and discards them.
 161     void discard(size_t z, command_queue &queue)
 162     {
 163         generate(discard_iterator(0), discard_iterator(z), queue);
 164     }
 165
 166     /// \internal_ (deprecated)
 167     template<class OutputIterator>
 168     void fill(OutputIterator first, OutputIterator last, command_queue &queue)
 169     {
 170         generate(first, last, queue);
 171     }
 172
 173 private:
 174     /// \internal_
 175     void generate_state(command_queue &queue)
 176     {
 177         kernel generate_state_kernel =
 178             m_program.create_kernel("generate_state");
 179         generate_state_kernel.set_arg(0, m_state_buffer);
 180         queue.enqueue_task(generate_state_kernel);
 181     }
 182
 183     /// \internal_
 184     void load_program()
 185     {
 186         boost::shared_ptr<program_cache> cache =
 187             program_cache::get_global_cache(m_context);
 188
 189         std::string cache_key =
 190             std::string("__boost_mersenne_twister_engine_") + type_name<T>();
 191
 192         const char source[] =
 193             "static uint twiddle(uint u, uint v)\n"
 194             "{\n"
 195             "    return (((u & 0x80000000U) | (v & 0x7FFFFFFFU)) >> 1) ^\n"
 196             "           ((v & 1U) ? 0x9908B0DFU : 0x0U);\n"
 197             "}\n"
 198
 199             "__kernel void generate_state(__global uint *state)\n"
 200             "{\n"
 201             "    const uint n = 624;\n"
 202             "    const uint m = 397;\n"
 203             "    for(uint i = 0; i < (n - m); i++)\n"
 204             "        state[i] = state[i+m] ^ twiddle(state[i], state[i+1]);\n"
 205             "    for(uint i = n - m; i < (n - 1); i++)\n"
 206             "        state[i] = state[i+m-n] ^ twiddle(state[i], state[i+1]);\n"
 207             "    state[n-1] = state[m-1] ^ twiddle(state[n-1], state[0]);\n"
 208             "}\n"
 209
 210             "__kernel void seed(const uint s, __global uint *state)\n"
 211             "{\n"
 212             "    const uint n = 624;\n"
 213             "    state[0] = s & 0xFFFFFFFFU;\n"
 214             "    for(uint i = 1; i < n; i++){\n"
 215             "        state[i] = 1812433253U * (state[i-1] ^ (state[i-1] >> 30)) + i;\n"
 216             "        state[i] &= 0xFFFFFFFFU;\n"
 217             "    }\n"
 218             "    generate_state(state);\n"
 219             "}\n"
 220
 221             "static uint random_number(__global uint *state, const uint p)\n"
 222             "{\n"
 223             "    uint x = state[p];\n"
 224             "    x ^= (x >> 11);\n"
 225             "    x ^= (x << 7) & 0x9D2C5680U;\n"
 226             "    x ^= (x << 15) & 0xEFC60000U;\n"
 227             "    return x ^ (x >> 18);\n"
 228             "}\n"
 229
 230             "__kernel void fill(__global uint *state,\n"
 231             "                   const uint state_index,\n"
 232             "                   __global uint *vector,\n"
 233             "                   const uint offset)\n"
 234             "{\n"
 235             "    const uint i = get_global_id(0);\n"
 236             "    vector[offset+i] = random_number(state, state_index + i);\n"
 237             "}\n";
 238
 239         m_program = cache->get_or_build(cache_key, std::string(), source, m_context);
 240     }
 241
 242 private:
 243     context m_context;
 244     size_t m_state_index;
 245     program m_program;
 246     buffer m_state_buffer;
 247 };
 248
 249 typedef mersenne_twister_engine<uint_> mt19937;
 250
 251 } // end compute namespace
 252 } // end boost namespace
 253
 254 #endif // BOOST_COMPUTE_RANDOM_MERSENNE_TWISTER_ENGINE_HPP