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 */ param_type(int dim_arg=2)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. */ dim() const62 int dim() const { return _dim; } 63 64 /** Writes the parameters to a @c std::ostream. */ BOOST_RANDOM_DETAIL_OSTREAM_OPERATOR(os,param_type,parm)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. */ BOOST_RANDOM_DETAIL_ISTREAM_OPERATOR(is,param_type,parm)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. */ BOOST_RANDOM_DETAIL_EQUALITY_OPERATOR(param_type,lhs,rhs)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 */ uniform_on_sphere(int dim_arg=2)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 */ uniform_on_sphere(const param_type & parm)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. */ dim() const107 int dim() const { return _dim; } 108 109 /** Returns the parameters of the distribution. */ param() const110 param_type param() const { return param_type(_dim); } 111 /** Sets the parameters of the distribution. */ param(const param_type & parm)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 */ BOOST_PREVENT_MACRO_SUBSTITUTION() const126 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 */ BOOST_PREVENT_MACRO_SUBSTITUTION() const142 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 */ reset()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> operator ()(Engine & eng)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 break; 176 } 177 case 2: 178 { 179 uniform_01<RealType> uniform; 180 RealType sqsum; 181 RealType x, y; 182 do { 183 x = uniform(eng) * 2 - 1; 184 y = uniform(eng) * 2 - 1; 185 sqsum = x*x + y*y; 186 } while(sqsum == 0 || sqsum > 1); 187 RealType mult = 1/sqrt(sqsum); 188 typename Cont::iterator iter = _container.begin(); 189 *iter = x * mult; 190 iter++; 191 *iter = y * mult; 192 break; 193 } 194 case 3: 195 { 196 uniform_01<RealType> uniform; 197 RealType sqsum; 198 RealType x, y; 199 do { 200 x = uniform(eng) * 2 - 1; 201 y = uniform(eng) * 2 - 1; 202 sqsum = x*x + y*y; 203 } while(sqsum > 1); 204 RealType mult = 2 * sqrt(1 - sqsum); 205 typename Cont::iterator iter = _container.begin(); 206 *iter = x * mult; 207 ++iter; 208 *iter = y * mult; 209 ++iter; 210 *iter = 2 * sqsum - 1; 211 break; 212 } 213 default: 214 { 215 detail::unit_normal_distribution<RealType> normal; 216 RealType sqsum; 217 do { 218 sqsum = 0; 219 for(typename Cont::iterator it = _container.begin(); 220 it != _container.end(); 221 ++it) { 222 RealType val = normal(eng); 223 *it = val; 224 sqsum += val * val; 225 } 226 } while(sqsum == 0); 227 // for all i: result[i] /= sqrt(sqsum) 228 RealType inverse_distance = 1 / sqrt(sqsum); 229 for(typename Cont::iterator it = _container.begin(); 230 it != _container.end(); 231 ++it) { 232 *it *= inverse_distance; 233 } 234 } 235 } 236 return _container; 237 } 238 239 /** 240 * Returns a point uniformly distributed over the surface of 241 * a sphere of dimension param.dim(). 242 */ 243 template<class Engine> operator ()(Engine & eng,const param_type & parm) const244 result_type operator()(Engine& eng, const param_type& parm) const 245 { 246 return uniform_on_sphere(parm)(eng); 247 } 248 249 /** Writes the distribution to a @c std::ostream. */ BOOST_RANDOM_DETAIL_OSTREAM_OPERATOR(os,uniform_on_sphere,sd)250 BOOST_RANDOM_DETAIL_OSTREAM_OPERATOR(os, uniform_on_sphere, sd) 251 { 252 os << sd._dim; 253 return os; 254 } 255 256 /** Reads the distribution from a @c std::istream. */ BOOST_RANDOM_DETAIL_ISTREAM_OPERATOR(is,uniform_on_sphere,sd)257 BOOST_RANDOM_DETAIL_ISTREAM_OPERATOR(is, uniform_on_sphere, sd) 258 { 259 is >> sd._dim; 260 sd._container.resize(sd._dim); 261 return is; 262 } 263 264 /** 265 * Returns true if the two distributions will produce identical 266 * sequences of values, given equal generators. 267 */ BOOST_RANDOM_DETAIL_EQUALITY_OPERATOR(uniform_on_sphere,lhs,rhs)268 BOOST_RANDOM_DETAIL_EQUALITY_OPERATOR(uniform_on_sphere, lhs, rhs) 269 { return lhs._dim == rhs._dim; } 270 271 /** 272 * Returns true if the two distributions may produce different 273 * sequences of values, given equal generators. 274 */ 275 BOOST_RANDOM_DETAIL_INEQUALITY_OPERATOR(uniform_on_sphere) 276 277 private: 278 result_type _container; 279 int _dim; 280 }; 281 282 } // namespace random 283 284 using random::uniform_on_sphere; 285 286 } // namespace boost 287 288 #endif // BOOST_RANDOM_UNIFORM_ON_SPHERE_HPP 289