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
28namespace boost {
29namespace 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 */
39template<class RealType = double, class Cont = std::vector<RealType> >
40class uniform_on_sphere
41{
42public:
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 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 std::transform(_container.begin(), _container.end(), _container.begin(),
229 std::bind2nd(std::multiplies<RealType>(), 1/sqrt(sqsum)));
230 }
231 }
232 return _container;
233 }
234
235 /**
236 * Returns a point uniformly distributed over the surface of
237 * a sphere of dimension param.dim().
238 */
239 template<class Engine>
240 result_type operator()(Engine& eng, const param_type& parm) const
241 {
242 return uniform_on_sphere(parm)(eng);
243 }
244
245 /** Writes the distribution to a @c std::ostream. */
246 BOOST_RANDOM_DETAIL_OSTREAM_OPERATOR(os, uniform_on_sphere, sd)
247 {
248 os << sd._dim;
249 return os;
250 }
251
252 /** Reads the distribution from a @c std::istream. */
253 BOOST_RANDOM_DETAIL_ISTREAM_OPERATOR(is, uniform_on_sphere, sd)
254 {
255 is >> sd._dim;
256 sd._container.resize(sd._dim);
257 return is;
258 }
259
260 /**
261 * Returns true if the two distributions will produce identical
262 * sequences of values, given equal generators.
263 */
264 BOOST_RANDOM_DETAIL_EQUALITY_OPERATOR(uniform_on_sphere, lhs, rhs)
265 { return lhs._dim == rhs._dim; }
266
267 /**
268 * Returns true if the two distributions may produce different
269 * sequences of values, given equal generators.
270 */
271 BOOST_RANDOM_DETAIL_INEQUALITY_OPERATOR(uniform_on_sphere)
272
273private:
274 result_type _container;
275 int _dim;
276};
277
278} // namespace random
279
280using random::uniform_on_sphere;
281
282} // namespace boost
283
284#endif // BOOST_RANDOM_UNIFORM_ON_SPHERE_HPP
285