compute/random/mersenne_twister_engine.hpp

0001 //---------------------------------------------------------------------------//
0002 // Copyright (c) 2013 Kyle Lutz <kyle.r.lutz@gmail.com>
0003 //
0004 // Distributed under the Boost Software License, Version 1.0
0005 // See accompanying file LICENSE_1_0.txt or copy at
0006 // http://www.boost.org/LICENSE_1_0.txt
0007 //
0008 // See http://boostorg.github.com/compute for more information.
0009 //---------------------------------------------------------------------------//
0010
0011 #ifndef BOOST_COMPUTE_RANDOM_MERSENNE_TWISTER_ENGINE_HPP
0012 #define BOOST_COMPUTE_RANDOM_MERSENNE_TWISTER_ENGINE_HPP
0013
0014 #include <algorithm>
0015
0016 #include <boost/compute/types.hpp>
0017 #include <boost/compute/buffer.hpp>
0018 #include <boost/compute/kernel.hpp>
0019 #include <boost/compute/context.hpp>
0020 #include <boost/compute/program.hpp>
0021 #include <boost/compute/command_queue.hpp>
0022 #include <boost/compute/algorithm/transform.hpp>
0023 #include <boost/compute/container/vector.hpp>
0024 #include <boost/compute/detail/iterator_range_size.hpp>
0025 #include <boost/compute/iterator/discard_iterator.hpp>
0026 #include <boost/compute/utility/program_cache.hpp>
0027
0028 namespace boost {
0029 namespace compute {
0030
0031 /// \class mersenne_twister_engine
0032 /// \brief Mersenne twister pseudorandom number generator.
0033 template<class T>
0034 class mersenne_twister_engine
0035 {
0036 public:
0037     typedef T result_type;
0038     static const T default_seed = 5489U;
0039     static const T n = 624;
0040     static const T m = 397;
0041
0042     /// Creates a new mersenne_twister_engine and seeds it with \p value.
0043     explicit mersenne_twister_engine(command_queue &queue,
0044                                      result_type value = default_seed)
0045         : m_context(queue.get_context()),
0046           m_state_buffer(m_context, n * sizeof(result_type))
0047     {
0048         // setup program
0049         load_program();
0050
0051         // seed state
0052         seed(value, queue);
0053     }
0054
0055     /// Creates a new mersenne_twister_engine object as a copy of \p other.
0056     mersenne_twister_engine(const mersenne_twister_engine<T> &other)
0057         : m_context(other.m_context),
0058           m_state_index(other.m_state_index),
0059           m_program(other.m_program),
0060           m_state_buffer(other.m_state_buffer)
0061     {
0062     }
0063
0064     /// Copies \p other to \c *this.
0065     mersenne_twister_engine<T>& operator=(const mersenne_twister_engine<T> &other)
0066     {
0067         if(this != &other){
0068             m_context = other.m_context;
0069             m_state_index = other.m_state_index;
0070             m_program = other.m_program;
0071             m_state_buffer = other.m_state_buffer;
0072         }
0073
0074         return *this;
0075     }
0076
0077     /// Destroys the mersenne_twister_engine object.
0078     ~mersenne_twister_engine()
0079     {
0080     }
0081
0082     /// Seeds the random number generator with \p value.
0083     ///
0084     /// \param value seed value for the random-number generator
0085     /// \param queue command queue to perform the operation
0086     ///
0087     /// If no seed value is provided, \c default_seed is used.
0088     void seed(result_type value, command_queue &queue)
0089     {
0090         kernel seed_kernel = m_program.create_kernel("seed");
0091         seed_kernel.set_arg(0, value);
0092         seed_kernel.set_arg(1, m_state_buffer);
0093
0094         queue.enqueue_task(seed_kernel);
0095
0096         m_state_index = 0;
0097     }
0098
0099     /// \overload
0100     void seed(command_queue &queue)
0101     {
0102         seed(default_seed, queue);
0103     }
0104
0105     /// Generates random numbers and stores them to the range [\p first, \p last).
0106     template<class OutputIterator>
0107     void generate(OutputIterator first, OutputIterator last, command_queue &queue)
0108     {
0109         const size_t size = detail::iterator_range_size(first, last);
0110
0111         kernel fill_kernel(m_program, "fill");
0112         fill_kernel.set_arg(0, m_state_buffer);
0113         fill_kernel.set_arg(2, first.get_buffer());
0114
0115         size_t offset = 0;
0116         size_t &p = m_state_index;
0117
0118         for(;;){
0119             size_t count = 0;
0120             if(size > n){
0121                 count = (std::min)(static_cast<size_t>(n), size - offset);
0122             }
0123             else {
0124                 count = size;
0125             }
0126             fill_kernel.set_arg(1, static_cast<const uint_>(p));
0127             fill_kernel.set_arg(3, static_cast<const uint_>(offset));
0128             queue.enqueue_1d_range_kernel(fill_kernel, 0, count, 0);
0129
0130             p += count;
0131             offset += count;
0132
0133             if(offset >= size){
0134                 break;
0135             }
0136
0137             generate_state(queue);
0138             p = 0;
0139         }
0140     }
0141
0142     /// \internal_
0143     void generate(discard_iterator first, discard_iterator last, command_queue &queue)
0144     {
0145         (void) queue;
0146
0147         m_state_index += std::distance(first, last);
0148     }
0149
0150     /// Generates random numbers, transforms them with \p op, and then stores
0151     /// them to the range [\p first, \p last).
0152     template<class OutputIterator, class Function>
0153     void generate(OutputIterator first, OutputIterator last, Function op, command_queue &queue)
0154     {
0155         vector<T> tmp(std::distance(first, last), queue.get_context());
0156         generate(tmp.begin(), tmp.end(), queue);
0157         transform(tmp.begin(), tmp.end(), first, op, queue);
0158     }
0159
0160     /// Generates \p z random numbers and discards them.
0161     void discard(size_t z, command_queue &queue)
0162     {
0163         generate(discard_iterator(0), discard_iterator(z), queue);
0164     }
0165
0166     /// \internal_ (deprecated)
0167     template<class OutputIterator>
0168     void fill(OutputIterator first, OutputIterator last, command_queue &queue)
0169     {
0170         generate(first, last, queue);
0171     }
0172
0173 private:
0174     /// \internal_
0175     void generate_state(command_queue &queue)
0176     {
0177         kernel generate_state_kernel =
0178             m_program.create_kernel("generate_state");
0179         generate_state_kernel.set_arg(0, m_state_buffer);
0180         queue.enqueue_task(generate_state_kernel);
0181     }
0182
0183     /// \internal_
0184     void load_program()
0185     {
0186         boost::shared_ptr<program_cache> cache =
0187             program_cache::get_global_cache(m_context);
0188
0189         std::string cache_key =
0190             std::string("__boost_mersenne_twister_engine_") + type_name<T>();
0191
0192         const char source[] =
0193             "static uint twiddle(uint u, uint v)\n"
0194             "{\n"
0195             "    return (((u & 0x80000000U) | (v & 0x7FFFFFFFU)) >> 1) ^\n"
0196             "           ((v & 1U) ? 0x9908B0DFU : 0x0U);\n"
0197             "}\n"
0198
0199             "__kernel void generate_state(__global uint *state)\n"
0200             "{\n"
0201             "    const uint n = 624;\n"
0202             "    const uint m = 397;\n"
0203             "    for(uint i = 0; i < (n - m); i++)\n"
0204             "        state[i] = state[i+m] ^ twiddle(state[i], state[i+1]);\n"
0205             "    for(uint i = n - m; i < (n - 1); i++)\n"
0206             "        state[i] = state[i+m-n] ^ twiddle(state[i], state[i+1]);\n"
0207             "    state[n-1] = state[m-1] ^ twiddle(state[n-1], state[0]);\n"
0208             "}\n"
0209
0210             "__kernel void seed(const uint s, __global uint *state)\n"
0211             "{\n"
0212             "    const uint n = 624;\n"
0213             "    state[0] = s & 0xFFFFFFFFU;\n"
0214             "    for(uint i = 1; i < n; i++){\n"
0215             "        state[i] = 1812433253U * (state[i-1] ^ (state[i-1] >> 30)) + i;\n"
0216             "        state[i] &= 0xFFFFFFFFU;\n"
0217             "    }\n"
0218             "    generate_state(state);\n"
0219             "}\n"
0220
0221             "static uint random_number(__global uint *state, const uint p)\n"
0222             "{\n"
0223             "    uint x = state[p];\n"
0224             "    x ^= (x >> 11);\n"
0225             "    x ^= (x << 7) & 0x9D2C5680U;\n"
0226             "    x ^= (x << 15) & 0xEFC60000U;\n"
0227             "    return x ^ (x >> 18);\n"
0228             "}\n"
0229
0230             "__kernel void fill(__global uint *state,\n"
0231             "                   const uint state_index,\n"
0232             "                   __global uint *vector,\n"
0233             "                   const uint offset)\n"
0234             "{\n"
0235             "    const uint i = get_global_id(0);\n"
0236             "    vector[offset+i] = random_number(state, state_index + i);\n"
0237             "}\n";
0238
0239         m_program = cache->get_or_build(cache_key, std::string(), source, m_context);
0240     }
0241
0242 private:
0243     context m_context;
0244     size_t m_state_index;
0245     program m_program;
0246     buffer m_state_buffer;
0247 };
0248
0249 typedef mersenne_twister_engine<uint_> mt19937;
0250
0251 } // end compute namespace
0252 } // end boost namespace
0253
0254 #endif // BOOST_COMPUTE_RANDOM_MERSENNE_TWISTER_ENGINE_HPP