1 //---------------------------------------------------------------------------// 2 // Copyright (c) 2014 Roshan <thisisroshansmail@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_LINEAR_CONGRUENTIAL_ENGINE_HPP 12 #define BOOST_COMPUTE_RANDOM_LINEAR_CONGRUENTIAL_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 /// 32 /// \class linear_congruential_engine 33 /// \brief 'Quick and Dirty' linear congruential engine 34 /// 35 /// Quick and dirty linear congruential engine to generate low quality 36 /// random numbers very quickly. For uses in which good quality of random 37 /// numbers is required(Monte-Carlo Simulations), use other engines like 38 /// Mersenne Twister instead. 39 /// 40 template<class T = uint_> 41 class linear_congruential_engine 42 { 43 public: 44 typedef T result_type; 45 static const T default_seed = 1; 46 static const T a = 1099087573; 47 static const size_t threads = 1024; 48 49 /// Creates a new linear_congruential_engine and seeds it with \p value. linear_congruential_engine(command_queue & queue,result_type value=default_seed)50 explicit linear_congruential_engine(command_queue &queue, 51 result_type value = default_seed) 52 : m_context(queue.get_context()), 53 m_multiplicands(m_context, threads * sizeof(result_type)) 54 { 55 // setup program 56 load_program(); 57 58 // seed state 59 seed(value, queue); 60 61 // generate multiplicands 62 generate_multiplicands(queue); 63 } 64 65 /// Creates a new linear_congruential_engine object as a copy of \p other. linear_congruential_engine(const linear_congruential_engine<T> & other)66 linear_congruential_engine(const linear_congruential_engine<T> &other) 67 : m_context(other.m_context), 68 m_program(other.m_program), 69 m_seed(other.m_seed), 70 m_multiplicands(other.m_multiplicands) 71 { 72 } 73 74 /// Copies \p other to \c *this. 75 linear_congruential_engine<T>& operator =(const linear_congruential_engine<T> & other)76 operator=(const linear_congruential_engine<T> &other) 77 { 78 if(this != &other){ 79 m_context = other.m_context; 80 m_program = other.m_program; 81 m_seed = other.m_seed; 82 m_multiplicands = other.m_multiplicands; 83 } 84 85 return *this; 86 } 87 88 /// Destroys the linear_congruential_engine object. ~linear_congruential_engine()89 ~linear_congruential_engine() 90 { 91 } 92 93 /// Seeds the random number generator with \p value. 94 /// 95 /// \param value seed value for the random-number generator 96 /// \param queue command queue to perform the operation 97 /// 98 /// If no seed value is provided, \c default_seed is used. seed(result_type value,command_queue & queue)99 void seed(result_type value, command_queue &queue) 100 { 101 (void) queue; 102 103 m_seed = value; 104 } 105 106 /// \overload seed(command_queue & queue)107 void seed(command_queue &queue) 108 { 109 seed(default_seed, queue); 110 } 111 112 /// Generates random numbers and stores them to the range [\p first, \p last). 113 template<class OutputIterator> generate(OutputIterator first,OutputIterator last,command_queue & queue)114 void generate(OutputIterator first, OutputIterator last, command_queue &queue) 115 { 116 size_t size = detail::iterator_range_size(first, last); 117 118 kernel fill_kernel(m_program, "fill"); 119 fill_kernel.set_arg(1, m_multiplicands); 120 fill_kernel.set_arg(2, first.get_buffer()); 121 122 size_t offset = 0; 123 124 for(;;){ 125 size_t count = 0; 126 if(size > threads){ 127 count = (std::min)(static_cast<size_t>(threads), size - offset); 128 } 129 else { 130 count = size; 131 } 132 fill_kernel.set_arg(0, static_cast<const uint_>(m_seed)); 133 fill_kernel.set_arg(3, static_cast<const uint_>(offset)); 134 queue.enqueue_1d_range_kernel(fill_kernel, 0, count, 0); 135 136 offset += count; 137 138 if(offset >= size){ 139 break; 140 } 141 142 update_seed(queue); 143 } 144 } 145 146 /// \internal_ generate(discard_iterator first,discard_iterator last,command_queue & queue)147 void generate(discard_iterator first, discard_iterator last, command_queue &queue) 148 { 149 (void) queue; 150 151 size_t size = detail::iterator_range_size(first, last); 152 uint_ max_mult = 153 detail::read_single_value<T>(m_multiplicands, threads-1, queue); 154 while(size >= threads) { 155 m_seed *= max_mult; 156 size -= threads; 157 } 158 m_seed *= 159 detail::read_single_value<T>(m_multiplicands, size-1, queue); 160 } 161 162 /// Generates random numbers, transforms them with \p op, and then stores 163 /// them to the range [\p first, \p last). 164 template<class OutputIterator, class Function> generate(OutputIterator first,OutputIterator last,Function op,command_queue & queue)165 void generate(OutputIterator first, OutputIterator last, Function op, command_queue &queue) 166 { 167 vector<T> tmp(std::distance(first, last), queue.get_context()); 168 generate(tmp.begin(), tmp.end(), queue); 169 transform(tmp.begin(), tmp.end(), first, op, queue); 170 } 171 172 /// Generates \p z random numbers and discards them. discard(size_t z,command_queue & queue)173 void discard(size_t z, command_queue &queue) 174 { 175 generate(discard_iterator(0), discard_iterator(z), queue); 176 } 177 178 private: 179 /// \internal_ 180 /// Generates the multiplicands for each thread generate_multiplicands(command_queue & queue)181 void generate_multiplicands(command_queue &queue) 182 { 183 kernel multiplicand_kernel = 184 m_program.create_kernel("multiplicand"); 185 multiplicand_kernel.set_arg(0, m_multiplicands); 186 187 queue.enqueue_task(multiplicand_kernel); 188 } 189 190 /// \internal_ update_seed(command_queue & queue)191 void update_seed(command_queue &queue) 192 { 193 m_seed *= 194 detail::read_single_value<T>(m_multiplicands, threads-1, queue); 195 } 196 197 /// \internal_ load_program()198 void load_program() 199 { 200 boost::shared_ptr<program_cache> cache = 201 program_cache::get_global_cache(m_context); 202 203 std::string cache_key = 204 std::string("__boost_linear_congruential_engine_") + type_name<T>(); 205 206 const char source[] = 207 "__kernel void multiplicand(__global uint *multiplicands)\n" 208 "{\n" 209 " uint a = 1099087573;\n" 210 " multiplicands[0] = a;\n" 211 " for(uint i = 1; i < 1024; i++){\n" 212 " multiplicands[i] = a * multiplicands[i-1];\n" 213 " }\n" 214 "}\n" 215 216 "__kernel void fill(const uint seed,\n" 217 " __global uint *multiplicands,\n" 218 " __global uint *result," 219 " const uint offset)\n" 220 "{\n" 221 " const uint i = get_global_id(0);\n" 222 " result[offset+i] = seed * multiplicands[i];\n" 223 "}\n"; 224 225 m_program = cache->get_or_build(cache_key, std::string(), source, m_context); 226 } 227 228 private: 229 context m_context; 230 program m_program; 231 T m_seed; 232 buffer m_multiplicands; 233 }; 234 235 } // end compute namespace 236 } // end boost namespace 237 238 #endif // BOOST_COMPUTE_RANDOM_LINEAR_CONGRUENTIAL_ENGINE_HPP 239