Optimize uniform_on_sphere for dims <= 3. Fixes #1059.

include/boost/random/uniform_on_sphere.hpp

@@ -152,7 +152,7 @@ class uniform_on_sphere
      * Effects: Subsequent uses of the distribution do not depend
      * on values produced by any engine prior to invoking reset.
      */
-    void reset() { _normal.reset(); }
+    void reset() {}
 
     /**
      * Returns a point uniformly distributed over the surface of
@@ -161,21 +161,72 @@ class uniform_on_sphere
     template<class Engine>
     const result_type & operator()(Engine& eng)
     {
-        if (!_container.empty()) {
-            RealType sqsum = 0;
-            do {
-                for(typename Cont::iterator it = _container.begin();
-                    it != _container.end();
-                    ++it) {
-                    RealType val = _normal(eng);
-                    *it = val;
-                    sqsum += val * val;
+        using std::sqrt;
+        switch(_dim)
+        {
+        case 0: break;
+        case 1:
+            {
+                if(uniform_01<RealType>()(eng) < 0.5) {
+                    *_container.begin() = -1;
+                } else {
+                    *_container.begin() = 1;
                 }
-            } while(sqsum == 0);
-            using std::sqrt;
-            // for all i: result[i] /= sqrt(sqsum)
-            std::transform(_container.begin(), _container.end(), _container.begin(),
-                           std::bind2nd(std::divides<RealType>(), sqrt(sqsum)));
+            }
+        case 2:
+            {
+                uniform_01<RealType> uniform;
+                RealType sqsum;
+                RealType x, y;
+                do {
+                    x = uniform(eng) * 2 - 1;
+                    y = uniform(eng) * 2 - 1;
+                    sqsum = x*x + y*y;
+                } while(sqsum == 0 || sqsum > 1);
+                RealType mult = 1/sqrt(sqsum);
+                typename Cont::iterator iter = _container.begin();
+                *iter = x * mult;
+                iter++;
+                *iter = y * mult;
+                break;
+            }
+        case 3:
+            {
+                uniform_01<RealType> uniform;
+                RealType sqsum;
+                RealType x, y, z;
+                do {
+                    x = uniform(eng) * 2 - 1;
+                    y = uniform(eng) * 2 - 1;
+                    sqsum = x*x + y*y;
+                } while(sqsum > 1);
+                RealType mult = 2 * sqrt(1 - sqsum);
+                typename Cont::iterator iter = _container.begin();
+                *iter = x * mult;
+                ++iter;
+                *iter = y * mult;
+                ++iter;
+                *iter = 2 * sqsum - 1;
+                break;
+            }
+        default:
+            {
+                detail::unit_normal_distribution<RealType> normal;
+                RealType sqsum;
+                do {
+                    sqsum = 0;
+                    for(typename Cont::iterator it = _container.begin();
+                        it != _container.end();
+                        ++it) {
+                        RealType val = normal(eng);
+                        *it = val;
+                        sqsum += val * val;
+                    }
+                } while(sqsum == 0);
+                // for all i: result[i] /= sqrt(sqsum)
+                std::transform(_container.begin(), _container.end(), _container.begin(),
+                               std::bind2nd(std::multiplies<RealType>(), 1/sqrt(sqsum)));
+            }
         }
         return _container;
     }
@@ -210,7 +261,7 @@ class uniform_on_sphere
      * sequences of values, given equal generators.
      */
     BOOST_RANDOM_DETAIL_EQUALITY_OPERATOR(uniform_on_sphere, lhs, rhs)
-    { return lhs._dim == rhs._dim && lhs._normal == rhs._normal; }
+    { return lhs._dim == rhs._dim; }
 
     /**
      * Returns true if the two distributions may produce different
@@ -219,7 +270,6 @@ class uniform_on_sphere
     BOOST_RANDOM_DETAIL_INEQUALITY_OPERATOR(uniform_on_sphere)
 
 private:
-    normal_distribution<RealType> _normal;
     result_type _container;
     int _dim;
 };

test/Jamfile.v2

@@ -107,6 +107,7 @@ run test_uniform_int.cpp ;
 run test_uniform_int_distribution.cpp /boost//unit_test_framework ;
 run test_uniform_real.cpp ;
 run test_uniform_real_distribution.cpp /boost//unit_test_framework ;
+run test_uniform_on_sphere.cpp ;
 run test_uniform_on_sphere_distribution.cpp /boost//unit_test_framework ;
 run test_uniform_smallint.cpp ;
 run test_uniform_smallint_distribution.cpp /boost//unit_test_framework ;

test/test_uniform_on_sphere.cpp

@@ -0,0 +1,45 @@
+/* test_uniform_on_sphere.cpp
+ *
+ * Copyright Steven Watanabe 2011
+ * Distributed under the Boost Software License, Version 1.0. (See
+ * accompanying file LICENSE_1_0.txt or copy at
+ * http://www.boost.org/LICENSE_1_0.txt)
+ *
+ * $Id$
+ *
+ */
+
+#include <boost/random/uniform_on_sphere.hpp>
+#include <boost/random/uniform_int.hpp>
+#include <boost/math/distributions/uniform.hpp>
+#include <cmath>
+
+class uniform_on_sphere_test {
+public:
+    typedef double result_type;
+    uniform_on_sphere_test(int dims, int x, int y)
+        : impl(dims), idx1(x), idx2(y) {}
+    template<class Engine>
+    result_type operator()(Engine& rng) {
+        const boost::random::uniform_on_sphere<>::result_type& tmp = impl(rng);
+        // This should be uniformly distributed in [-pi,pi)
+        return std::atan2(tmp[idx1], tmp[idx2]);
+    }
+private:
+    boost::random::uniform_on_sphere<> impl;
+    int idx1, idx2;
+};
+
+static const double pi = 3.14159265358979323846;
+
+#define BOOST_RANDOM_DISTRIBUTION uniform_on_sphere_test
+#define BOOST_RANDOM_DISTRIBUTION_NAME uniform_on_sphere
+#define BOOST_MATH_DISTRIBUTION boost::math::uniform
+#define BOOST_RANDOM_ARG1_TYPE double
+#define BOOST_RANDOM_ARG1_NAME n
+#define BOOST_RANDOM_ARG1_DEFAULT 6
+#define BOOST_RANDOM_ARG1_DISTRIBUTION(n) boost::uniform_int<>(2, n)
+#define BOOST_RANDOM_DISTRIBUTION_INIT (n, 0, n-1)
+#define BOOST_MATH_DISTRIBUTION_INIT (-pi, pi)
+
+#include "test_real_distribution.ipp"