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. mersenne_twister_engine(command_queue & queue,result_type value=default_seed)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. mersenne_twister_engine(const mersenne_twister_engine<T> & 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. operator =(const mersenne_twister_engine<T> & other)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. ~mersenne_twister_engine()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. seed(result_type value,command_queue & queue)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 seed(command_queue & queue)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> generate(OutputIterator first,OutputIterator last,command_queue & queue)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_ generate(discard_iterator first,discard_iterator last,command_queue & queue)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> generate(OutputIterator first,OutputIterator last,Function op,command_queue & queue)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. discard(size_t z,command_queue & queue)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> fill(OutputIterator first,OutputIterator last,command_queue & queue)168 void fill(OutputIterator first, OutputIterator last, command_queue &queue) 169 { 170 generate(first, last, queue); 171 } 172 173 private: 174 /// \internal_ generate_state(command_queue & queue)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_ load_program()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 255