1 //---------------------------------------------------------------------------// 2 // Copyright (c) 2013-2014 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_NORMAL_DISTRIBUTION_HPP 12 #define BOOST_COMPUTE_RANDOM_NORMAL_DISTRIBUTION_HPP 13 14 #include <limits> 15 16 #include <boost/assert.hpp> 17 #include <boost/type_traits.hpp> 18 19 #include <boost/compute/command_queue.hpp> 20 #include <boost/compute/function.hpp> 21 #include <boost/compute/types/fundamental.hpp> 22 #include <boost/compute/type_traits/make_vector_type.hpp> 23 24 namespace boost { 25 namespace compute { 26 27 /// \class normal_distribution 28 /// \brief Produces random, normally-distributed floating-point numbers. 29 /// 30 /// The following example shows how to setup a normal distribution to 31 /// produce random \c float values centered at \c 5: 32 /// 33 /// \snippet test/test_normal_distribution.cpp generate 34 /// 35 /// \see default_random_engine, uniform_real_distribution 36 template<class RealType = float> 37 class normal_distribution 38 { 39 public: 40 typedef RealType result_type; 41 42 /// Creates a new normal distribution producing numbers with the given 43 /// \p mean and \p stddev. normal_distribution(RealType mean=0.f,RealType stddev=1.f)44 normal_distribution(RealType mean = 0.f, RealType stddev = 1.f) 45 : m_mean(mean), 46 m_stddev(stddev) 47 { 48 } 49 50 /// Destroys the normal distribution object. ~normal_distribution()51 ~normal_distribution() 52 { 53 } 54 55 /// Returns the mean value of the distribution. mean() const56 result_type mean() const 57 { 58 return m_mean; 59 } 60 61 /// Returns the standard-deviation of the distribution. stddev() const62 result_type stddev() const 63 { 64 return m_stddev; 65 } 66 67 /// Returns the minimum value of the distribution. BOOST_PREVENT_MACRO_SUBSTITUTION() const68 result_type min BOOST_PREVENT_MACRO_SUBSTITUTION () const 69 { 70 return -std::numeric_limits<RealType>::infinity(); 71 } 72 73 /// Returns the maximum value of the distribution. BOOST_PREVENT_MACRO_SUBSTITUTION() const74 result_type max BOOST_PREVENT_MACRO_SUBSTITUTION () const 75 { 76 return std::numeric_limits<RealType>::infinity(); 77 } 78 79 /// Generates normally-distributed floating-point numbers and stores 80 /// them to the range [\p first, \p last). 81 template<class OutputIterator, class Generator> generate(OutputIterator first,OutputIterator last,Generator & generator,command_queue & queue)82 void generate(OutputIterator first, 83 OutputIterator last, 84 Generator &generator, 85 command_queue &queue) 86 { 87 typedef typename make_vector_type<RealType, 2>::type RealType2; 88 89 size_t count = detail::iterator_range_size(first, last); 90 91 vector<uint_> tmp(count, queue.get_context()); 92 generator.generate(tmp.begin(), tmp.end(), queue); 93 94 BOOST_COMPUTE_FUNCTION(RealType2, box_muller, (const uint2_ x), 95 { 96 const RealType one = 1; 97 const RealType two = 2; 98 99 // Use nextafter to push values down into [0,1) range; without this, floating point rounding can 100 // lead to have x1 = 1, but that would lead to taking the log of 0, which would result in negative 101 // infinities; by pushing the values off 1 towards 0, we ensure this won't happen. 102 const RealType x1 = nextafter(x.x / (RealType) UINT_MAX, (RealType) 0); 103 const RealType x2 = x.y / (RealType) UINT_MAX; 104 105 const RealType rho = sqrt(-two * log(one-x1)); 106 107 const RealType z1 = rho * cos(two * M_PI_F * x2); 108 const RealType z2 = rho * sin(two * M_PI_F * x2); 109 110 return (RealType2)(MEAN, MEAN) + (RealType2)(z1, z2) * (RealType2)(STDDEV, STDDEV); 111 }); 112 113 box_muller.define("MEAN", boost::lexical_cast<std::string>(m_mean)); 114 box_muller.define("STDDEV", boost::lexical_cast<std::string>(m_stddev)); 115 box_muller.define("RealType", type_name<RealType>()); 116 box_muller.define("RealType2", type_name<RealType2>()); 117 118 transform( 119 make_buffer_iterator<uint2_>(tmp.get_buffer(), 0), 120 make_buffer_iterator<uint2_>(tmp.get_buffer(), count / 2), 121 make_buffer_iterator<RealType2>(first.get_buffer(), 0), 122 box_muller, 123 queue 124 ); 125 } 126 127 private: 128 RealType m_mean; 129 RealType m_stddev; 130 131 BOOST_STATIC_ASSERT_MSG( 132 boost::is_floating_point<RealType>::value, 133 "Template argument must be a floating point type" 134 ); 135 }; 136 137 } // end compute namespace 138 } // end boost namespace 139 140 #endif // BOOST_COMPUTE_RANDOM_NORMAL_DISTRIBUTION_HPP 141