git.sesse.net Git - casparcg/blob - dependencies64/boost/boost/random/uniform_on_sphere.hpp

   1 /* boost random/uniform_on_sphere.hpp header file
   2  *
   3  * Copyright Jens Maurer 2000-2001
   4  * Copyright Steven Watanabe 2011
   5  * Distributed under the Boost Software License, Version 1.0. (See
   6  * accompanying file LICENSE_1_0.txt or copy at
   7  * http://www.boost.org/LICENSE_1_0.txt)
   8  *
   9  * See http://www.boost.org for most recent version including documentation.
  10  *
  11  * $Id$
  12  *
  13  * Revision history
  14  *  2001-02-18  moved to individual header files
  15  */
  16
  17 #ifndef BOOST_RANDOM_UNIFORM_ON_SPHERE_HPP
  18 #define BOOST_RANDOM_UNIFORM_ON_SPHERE_HPP
  19
  20 #include <vector>
  21 #include <algorithm>     // std::transform
  22 #include <functional>    // std::bind2nd, std::divides
  23 #include <boost/assert.hpp>
  24 #include <boost/random/detail/config.hpp>
  25 #include <boost/random/detail/operators.hpp>
  26 #include <boost/random/normal_distribution.hpp>
  27
  28 namespace boost {
  29 namespace random {
  30
  31 /**
  32  * Instantiations of class template uniform_on_sphere model a
  33  * \random_distribution. Such a distribution produces random
  34  * numbers uniformly distributed on the unit sphere of arbitrary
  35  * dimension @c dim. The @c Cont template parameter must be a STL-like
  36  * container type with begin and end operations returning non-const
  37  * ForwardIterators of type @c Cont::iterator.
  38  */
  39 template<class RealType = double, class Cont = std::vector<RealType> >
  40 class uniform_on_sphere
  41 {
  42 public:
  43     typedef RealType input_type;
  44     typedef Cont result_type;
  45
  46     class param_type
  47     {
  48     public:
  49
  50         typedef uniform_on_sphere distribution_type;
  51
  52         /**
  53          * Constructs the parameters of a uniform_on_sphere
  54          * distribution, given the dimension of the sphere.
  55          */
  56         explicit param_type(int dim_arg = 2) : _dim(dim_arg)
  57         {
  58             BOOST_ASSERT(_dim >= 0);
  59         }
  60
  61         /** Returns the dimension of the sphere. */
  62         int dim() const { return _dim; }
  63
  64         /** Writes the parameters to a @c std::ostream. */
  65         BOOST_RANDOM_DETAIL_OSTREAM_OPERATOR(os, param_type, parm)
  66         {
  67             os << parm._dim;
  68             return os;
  69         }
  70
  71         /** Reads the parameters from a @c std::istream. */
  72         BOOST_RANDOM_DETAIL_ISTREAM_OPERATOR(is, param_type, parm)
  73         {
  74             is >> parm._dim;
  75             return is;
  76         }
  77
  78         /** Returns true if the two sets of parameters are equal. */
  79         BOOST_RANDOM_DETAIL_EQUALITY_OPERATOR(param_type, lhs, rhs)
  80         { return lhs._dim == rhs._dim; }
  81
  82         /** Returns true if the two sets of parameters are different. */
  83         BOOST_RANDOM_DETAIL_INEQUALITY_OPERATOR(param_type)
  84
  85     private:
  86         int _dim;
  87     };
  88
  89     /**
  90      * Constructs a @c uniform_on_sphere distribution.
  91      * @c dim is the dimension of the sphere.
  92      *
  93      * Requires: dim >= 0
  94      */
  95     explicit uniform_on_sphere(int dim_arg = 2)
  96       : _container(dim_arg), _dim(dim_arg) { }
  97
  98     /**
  99      * Constructs a @c uniform_on_sphere distribution from its parameters.
 100      */
 101     explicit uniform_on_sphere(const param_type& parm)
 102       : _container(parm.dim()), _dim(parm.dim()) { }
 103
 104     // compiler-generated copy ctor and assignment operator are fine
 105
 106     /** Returns the dimension of the sphere. */
 107     int dim() const { return _dim; }
 108
 109     /** Returns the parameters of the distribution. */
 110     param_type param() const { return param_type(_dim); }
 111     /** Sets the parameters of the distribution. */
 112     void param(const param_type& parm)
 113     {
 114         _dim = parm.dim();
 115         _container.resize(_dim);
 116     }
 117
 118     /**
 119      * Returns the smallest value that the distribution can produce.
 120      * Note that this is required to approximate the standard library's
 121      * requirements.  The behavior is defined according to lexicographical
 122      * comparison so that for a container type of std::vector,
 123      * dist.min() <= x <= dist.max() where x is any value produced
 124      * by the distribution.
 125      */
 126     result_type min BOOST_PREVENT_MACRO_SUBSTITUTION () const
 127     {
 128         result_type result(_dim);
 129         if(_dim != 0) {
 130             result.front() = RealType(-1.0);
 131         }
 132         return result;
 133     }
 134     /**
 135      * Returns the largest value that the distribution can produce.
 136      * Note that this is required to approximate the standard library's
 137      * requirements.  The behavior is defined according to lexicographical
 138      * comparison so that for a container type of std::vector,
 139      * dist.min() <= x <= dist.max() where x is any value produced
 140      * by the distribution.
 141      */
 142     result_type max BOOST_PREVENT_MACRO_SUBSTITUTION () const
 143     {
 144         result_type result(_dim);
 145         if(_dim != 0) {
 146             result.front() = RealType(1.0);
 147         }
 148         return result;
 149     }
 150
 151     /**
 152      * Effects: Subsequent uses of the distribution do not depend
 153      * on values produced by any engine prior to invoking reset.
 154      */
 155     void reset() {}
 156
 157     /**
 158      * Returns a point uniformly distributed over the surface of
 159      * a sphere of dimension dim().
 160      */
 161     template<class Engine>
 162     const result_type & operator()(Engine& eng)
 163     {
 164         using std::sqrt;
 165         switch(_dim)
 166         {
 167         case 0: break;
 168         case 1:
 169             {
 170                 if(uniform_01<RealType>()(eng) < 0.5) {
 171                     *_container.begin() = -1;
 172                 } else {
 173                     *_container.begin() = 1;
 174                 }
 175             }
 176         case 2:
 177             {
 178                 uniform_01<RealType> uniform;
 179                 RealType sqsum;
 180                 RealType x, y;
 181                 do {
 182                     x = uniform(eng) * 2 - 1;
 183                     y = uniform(eng) * 2 - 1;
 184                     sqsum = x*x + y*y;
 185                 } while(sqsum == 0 || sqsum > 1);
 186                 RealType mult = 1/sqrt(sqsum);
 187                 typename Cont::iterator iter = _container.begin();
 188                 *iter = x * mult;
 189                 iter++;
 190                 *iter = y * mult;
 191                 break;
 192             }
 193         case 3:
 194             {
 195                 uniform_01<RealType> uniform;
 196                 RealType sqsum;
 197                 RealType x, y;
 198                 do {
 199                     x = uniform(eng) * 2 - 1;
 200                     y = uniform(eng) * 2 - 1;
 201                     sqsum = x*x + y*y;
 202                 } while(sqsum > 1);
 203                 RealType mult = 2 * sqrt(1 - sqsum);
 204                 typename Cont::iterator iter = _container.begin();
 205                 *iter = x * mult;
 206                 ++iter;
 207                 *iter = y * mult;
 208                 ++iter;
 209                 *iter = 2 * sqsum - 1;
 210                 break;
 211             }
 212         default:
 213             {
 214                 detail::unit_normal_distribution<RealType> normal;
 215                 RealType sqsum;
 216                 do {
 217                     sqsum = 0;
 218                     for(typename Cont::iterator it = _container.begin();
 219                         it != _container.end();
 220                         ++it) {
 221                         RealType val = normal(eng);
 222                         *it = val;
 223                         sqsum += val * val;
 224                     }
 225                 } while(sqsum == 0);
 226                 // for all i: result[i] /= sqrt(sqsum)
 227                 std::transform(_container.begin(), _container.end(), _container.begin(),
 228                                std::bind2nd(std::multiplies<RealType>(), 1/sqrt(sqsum)));
 229             }
 230         }
 231         return _container;
 232     }
 233
 234     /**
 235      * Returns a point uniformly distributed over the surface of
 236      * a sphere of dimension param.dim().
 237      */
 238     template<class Engine>
 239     result_type operator()(Engine& eng, const param_type& parm) const
 240     {
 241         return uniform_on_sphere(parm)(eng);
 242     }
 243
 244     /** Writes the distribution to a @c std::ostream. */
 245     BOOST_RANDOM_DETAIL_OSTREAM_OPERATOR(os, uniform_on_sphere, sd)
 246     {
 247         os << sd._dim;
 248         return os;
 249     }
 250
 251     /** Reads the distribution from a @c std::istream. */
 252     BOOST_RANDOM_DETAIL_ISTREAM_OPERATOR(is, uniform_on_sphere, sd)
 253     {
 254         is >> sd._dim;
 255         sd._container.resize(sd._dim);
 256         return is;
 257     }
 258
 259     /**
 260      * Returns true if the two distributions will produce identical
 261      * sequences of values, given equal generators.
 262      */
 263     BOOST_RANDOM_DETAIL_EQUALITY_OPERATOR(uniform_on_sphere, lhs, rhs)
 264     { return lhs._dim == rhs._dim; }
 265
 266     /**
 267      * Returns true if the two distributions may produce different
 268      * sequences of values, given equal generators.
 269      */
 270     BOOST_RANDOM_DETAIL_INEQUALITY_OPERATOR(uniform_on_sphere)
 271
 272 private:
 273     result_type _container;
 274     int _dim;
 275 };
 276
 277 } // namespace random
 278
 279 using random::uniform_on_sphere;
 280
 281 } // namespace boost
 282
 283 #endif // BOOST_RANDOM_UNIFORM_ON_SPHERE_HPP