ceph/src/boost/libs/compute/include/boost/compute/random/normal_distribution.hpp

   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.
  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.
  51     ~normal_distribution()
  52     {
  53     }
  54
  55     /// Returns the mean value of the distribution.
  56     result_type mean() const
  57     {
  58         return m_mean;
  59     }
  60
  61     /// Returns the standard-deviation of the distribution.
  62     result_type stddev() const
  63     {
  64         return m_stddev;
  65     }
  66
  67     /// Returns the minimum value of the distribution.
  68     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.
  74     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>
  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