diff options
Diffstat (limited to 'boost/compute/random/linear_congruential_engine.hpp')
| -rw-r--r-- | boost/compute/random/linear_congruential_engine.hpp | 238 |
1 files changed, 238 insertions, 0 deletions
diff --git a/boost/compute/random/linear_congruential_engine.hpp b/boost/compute/random/linear_congruential_engine.hpp new file mode 100644 index 0000000000..173a6c2881 --- /dev/null +++ b/boost/compute/random/linear_congruential_engine.hpp @@ -0,0 +1,238 @@ +//---------------------------------------------------------------------------// +// Copyright (c) 2014 Roshan <thisisroshansmail@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_LINEAR_CONGRUENTIAL_ENGINE_HPP +#define BOOST_COMPUTE_RANDOM_LINEAR_CONGRUENTIAL_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 linear_congruential_engine +/// \brief 'Quick and Dirty' linear congruential engine +/// +/// Quick and dirty linear congruential engine to generate low quality +/// random numbers very quickly. For uses in which good quality of random +/// numbers is required(Monte-Carlo Simulations), use other engines like +/// Mersenne Twister instead. +/// +template<class T = uint_> +class linear_congruential_engine +{ +public: + typedef T result_type; + static const T default_seed = 1; + static const T a = 1099087573; + static const size_t threads = 1024; + + /// Creates a new linear_congruential_engine and seeds it with \p value. + explicit linear_congruential_engine(command_queue &queue, + result_type value = default_seed) + : m_context(queue.get_context()), + m_multiplicands(m_context, threads * sizeof(result_type)) + { + // setup program + load_program(); + + // seed state + seed(value, queue); + + // generate multiplicands + generate_multiplicands(queue); + } + + /// Creates a new linear_congruential_engine object as a copy of \p other. + linear_congruential_engine(const linear_congruential_engine<T> &other) + : m_context(other.m_context), + m_program(other.m_program), + m_seed(other.m_seed), + m_multiplicands(other.m_multiplicands) + { + } + + /// Copies \p other to \c *this. + linear_congruential_engine<T>& + operator=(const linear_congruential_engine<T> &other) + { + if(this != &other){ + m_context = other.m_context; + m_program = other.m_program; + m_seed = other.m_seed; + m_multiplicands = other.m_multiplicands; + } + + return *this; + } + + /// Destroys the linear_congruential_engine object. + ~linear_congruential_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) + { + (void) queue; + + m_seed = value; + } + + /// \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) + { + size_t size = detail::iterator_range_size(first, last); + + kernel fill_kernel(m_program, "fill"); + fill_kernel.set_arg(1, m_multiplicands); + fill_kernel.set_arg(2, first.get_buffer()); + + size_t offset = 0; + + for(;;){ + size_t count = 0; + if(size > threads){ + count = (std::min)(static_cast<size_t>(threads), size - offset); + } + else { + count = size; + } + fill_kernel.set_arg(0, static_cast<const uint_>(m_seed)); + fill_kernel.set_arg(3, static_cast<const uint_>(offset)); + queue.enqueue_1d_range_kernel(fill_kernel, 0, count, 0); + + offset += count; + + if(offset >= size){ + break; + } + + update_seed(queue); + } + } + + /// \internal_ + void generate(discard_iterator first, discard_iterator last, command_queue &queue) + { + (void) queue; + + size_t size = detail::iterator_range_size(first, last); + uint_ max_mult = + detail::read_single_value<T>(m_multiplicands, threads-1, queue); + while(size >= threads) { + m_seed *= max_mult; + size -= threads; + } + m_seed *= + detail::read_single_value<T>(m_multiplicands, size-1, queue); + } + + /// 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); + } + +private: + /// \internal_ + /// Generates the multiplicands for each thread + void generate_multiplicands(command_queue &queue) + { + kernel multiplicand_kernel = + m_program.create_kernel("multiplicand"); + multiplicand_kernel.set_arg(0, m_multiplicands); + + queue.enqueue_task(multiplicand_kernel); + } + + /// \internal_ + void update_seed(command_queue &queue) + { + m_seed *= + detail::read_single_value<T>(m_multiplicands, threads-1, queue); + } + + /// \internal_ + void load_program() + { + boost::shared_ptr<program_cache> cache = + program_cache::get_global_cache(m_context); + + std::string cache_key = + std::string("__boost_linear_congruential_engine_") + type_name<T>(); + + const char source[] = + "__kernel void multiplicand(__global uint *multiplicands)\n" + "{\n" + " uint a = 1099087573;\n" + " multiplicands[0] = a;\n" + " for(uint i = 1; i < 1024; i++){\n" + " multiplicands[i] = a * multiplicands[i-1];\n" + " }\n" + "}\n" + + "__kernel void fill(const uint seed,\n" + " __global uint *multiplicands,\n" + " __global uint *result," + " const uint offset)\n" + "{\n" + " const uint i = get_global_id(0);\n" + " result[offset+i] = seed * multiplicands[i];\n" + "}\n"; + + m_program = cache->get_or_build(cache_key, std::string(), source, m_context); + } + +private: + context m_context; + program m_program; + T m_seed; + buffer m_multiplicands; +}; + +} // end compute namespace +} // end boost namespace + +#endif // BOOST_COMPUTE_RANDOM_LINEAR_CONGRUENTIAL_ENGINE_HPP |