[PATCH] libstdc++: add uniform on sphere distribution
Ulrich Drepper
drepper@gmail.com
Sat Aug 9 16:56:00 GMT 2014
Marc Glisse <marc.glisse@inria.fr> writes:
> On Sat, 9 Aug 2014, Ulrich Drepper wrote:
> Yes, you still need the normalization step (divide by the norm).
I guess we can do this.
How about the patch below? Instead of specializing the entire class for
_Dimen==2 I've added a class at the implementation level.
I've also improved existing tests and add some new ones.
2014-08-09 Ulrich Drepper <drepper@gmail.com>
* include/ext/random.tcc (uniform_on_sphere_helper): Define.
(uniform_on_sphere_distribution::operator()): Use the new helper
class for the implementation.
* testsuite/ext/random/uniform_on_sphere_distribution/operators/
equal.cc: Remove bogus part of comment.
* testsuite/ext/random/uniform_on_sphere_distribution/operators/
inequal.cc: Likewise.
* testsuite/ext/random/uniform_on_sphere_distribution/operators/
serialize.cc: Add check to verify result of serialzation and
deserialization.
* testsuite/ext/random/uniform_on_sphere_distribution/operators/
generate.cc: New file.
diff --git a/libstdc++-v3/include/ext/random.tcc b/libstdc++-v3/include/ext/random.tcc
index 05361d8..d536ecb 100644
--- a/libstdc++-v3/include/ext/random.tcc
+++ b/libstdc++-v3/include/ext/random.tcc
@@ -1540,6 +1540,83 @@ _GLIBCXX_BEGIN_NAMESPACE_VERSION
}
+ namespace {
+
+ // Helper class for the uniform_on_sphere_distribution generation
+ // function.
+ template<std::size_t _Dimen, typename _RealType>
+ class uniform_on_sphere_helper
+ {
+ typedef typename uniform_on_sphere_distribution<_Dimen, _RealType>::result_type result_type;
+
+ public:
+ template<typename _NormalDistribution, typename _UniformRandomNumberGenerator>
+ result_type operator()(_NormalDistribution& __nd,
+ _UniformRandomNumberGenerator& __urng)
+ {
+ result_type __ret;
+ typename result_type::value_type __norm;
+
+ do
+ {
+ auto __sum = _RealType(0);
+
+ std::generate(__ret.begin(), __ret.end(),
+ [&__nd, &__urng, &__sum](){
+ _RealType __t = __nd(__urng);
+ __sum += __t * __t;
+ return __t; });
+ __norm = std::sqrt(__sum);
+ }
+ while (__norm == _RealType(0) || ! std::isfinite(__norm));
+
+ std::transform(__ret.begin(), __ret.end(), __ret.begin(),
+ [__norm](_RealType __val){ return __val / __norm; });
+
+ return __ret;
+ }
+ };
+
+
+ template<typename _RealType>
+ class uniform_on_sphere_helper<2, _RealType>
+ {
+ typedef typename uniform_on_sphere_distribution<2, _RealType>::
+ result_type result_type;
+
+ public:
+ template<typename _NormalDistribution,
+ typename _UniformRandomNumberGenerator>
+ result_type operator()(_NormalDistribution&,
+ _UniformRandomNumberGenerator& __urng)
+ {
+ result_type __ret;
+ _RealType __sq;
+ std::__detail::_Adaptor<_UniformRandomNumberGenerator,
+ _RealType> __aurng(__urng);
+
+ do
+ {
+ __ret[0] = __aurng();
+ __ret[1] = __aurng();
+
+ __sq = __ret[0] * __ret[0] + __ret[1] * __ret[1];
+ }
+ while (__sq == _RealType(0) || __sq > _RealType(1));
+
+ // Yes, we do not just use sqrt(__sq) because hypot() is more
+ // accurate.
+ auto __norm = std::hypot(__ret[0], __ret[1]);
+ __ret[0] /= __norm;
+ __ret[1] /= __norm;
+
+ return __ret;
+ }
+ };
+
+ }
+
+
template<std::size_t _Dimen, typename _RealType>
template<typename _UniformRandomNumberGenerator>
typename uniform_on_sphere_distribution<_Dimen, _RealType>::result_type
@@ -1547,18 +1624,8 @@ _GLIBCXX_BEGIN_NAMESPACE_VERSION
operator()(_UniformRandomNumberGenerator& __urng,
const param_type& __p)
{
- result_type __ret;
- _RealType __sum = _RealType(0);
-
- std::generate(__ret.begin(), __ret.end(),
- [&__urng, &__sum, this](){ _RealType __t = _M_nd(__urng);
- __sum += __t * __t;
- return __t; });
- auto __norm = std::sqrt(__sum);
- std::transform(__ret.begin(), __ret.end(), __ret.begin(),
- [__norm](_RealType __val){ return __val / __norm; });
-
- return __ret;
+ uniform_on_sphere_helper<_Dimen, _RealType> __helper;
+ return __helper(_M_nd, __urng);
}
template<std::size_t _Dimen, typename _RealType>
diff --git a/libstdc++-v3/testsuite/ext/random/uniform_on_sphere_distribution/operators/equal.cc b/libstdc++-v3/testsuite/ext/random/uniform_on_sphere_distribution/operators/equal.cc
index 35a024e..f5b8d17 100644
--- a/libstdc++-v3/testsuite/ext/random/uniform_on_sphere_distribution/operators/equal.cc
+++ b/libstdc++-v3/testsuite/ext/random/uniform_on_sphere_distribution/operators/equal.cc
@@ -20,7 +20,7 @@
// with this library; see the file COPYING3. If not see
// <http://www.gnu.org/licenses/>.
-// 26.5.8.4.5 Class template uniform_on_sphere_distribution [rand.dist.ext.uniform_on_sphere]
+// Class template uniform_on_sphere_distribution [rand.dist.ext.uniform_on_sphere]
#include <ext/random>
#include <testsuite_hooks.h>
diff --git a/libstdc++-v3/testsuite/ext/random/uniform_on_sphere_distribution/operators/inequal.cc b/libstdc++-v3/testsuite/ext/random/uniform_on_sphere_distribution/operators/inequal.cc
index 9f8e8c8..2675652 100644
--- a/libstdc++-v3/testsuite/ext/random/uniform_on_sphere_distribution/operators/inequal.cc
+++ b/libstdc++-v3/testsuite/ext/random/uniform_on_sphere_distribution/operators/inequal.cc
@@ -20,7 +20,7 @@
// with this library; see the file COPYING3. If not see
// <http://www.gnu.org/licenses/>.
-// 26.5.8.4.5 Class template uniform_on_sphere_distribution [rand.dist.ext.uniform_on_sphere]
+// Class template uniform_on_sphere_distribution [rand.dist.ext.uniform_on_sphere]
#include <ext/random>
#include <testsuite_hooks.h>
diff --git a/libstdc++-v3/testsuite/ext/random/uniform_on_sphere_distribution/operators/serialize.cc b/libstdc++-v3/testsuite/ext/random/uniform_on_sphere_distribution/operators/serialize.cc
index 80264ff..e9a758c 100644
--- a/libstdc++-v3/testsuite/ext/random/uniform_on_sphere_distribution/operators/serialize.cc
+++ b/libstdc++-v3/testsuite/ext/random/uniform_on_sphere_distribution/operators/serialize.cc
@@ -20,8 +20,8 @@
// with this library; see the file COPYING3. If not see
// <http://www.gnu.org/licenses/>.
-// 26.4.8.3.* Class template uniform_on_sphere_distribution [rand.dist.ext.uniform_on_sphere]
-// 26.4.2.4 Concept RandomNumberDistribution [rand.concept.dist]
+// Class template uniform_on_sphere_distribution [rand.dist.ext.uniform_on_sphere]
+// Concept RandomNumberDistribution [rand.concept.dist]
#include <ext/random>
#include <sstream>
@@ -40,6 +40,8 @@ test01()
str << u;
str >> v;
+
+ VERIFY( u == v );
}
int
--- /dev/null 2014-07-15 08:50:39.432896277 -0400
+++ b/libstdc++-v3/testsuite/ext/random/uniform_on_sphere_distribution/operators/generate.cc 2014-08-09 08:07:29.787244291 -0400
@@ -0,0 +1,60 @@
+// { dg-options "-std=gnu++11" }
+// { dg-require-cstdint "" }
+//
+// 2014-08-09 Ulrich Drepper <drepper@gmail.com>
+//
+// Copyright (C) 2014 Free Software Foundation, Inc.
+//
+// This file is part of the GNU ISO C++ Library. This library is free
+// software; you can redistribute it and/or modify it under the
+// terms of the GNU General Public License as published by the
+// Free Software Foundation; either version 3, or (at your option)
+// any later version.
+//
+// This library is distributed in the hope that it will be useful,
+// but WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+// GNU General Public License for more details.
+//
+// You should have received a copy of the GNU General Public License along
+// with this library; see the file COPYING3. If not see
+// <http://www.gnu.org/licenses/>.
+
+// Class template uniform_on_sphere_distribution [rand.dist.ext.uniform_on_sphere]
+// Concept RandomNumberDistribution [rand.concept.dist]
+
+#include <ext/random>
+#include <sstream>
+#include <testsuite_hooks.h>
+
+void
+test01()
+{
+ bool test [[gnu::unused]] = true;
+ std::minstd_rand0 rng;
+
+ __gnu_cxx::uniform_on_sphere_distribution<3> u3;
+
+ for (size_t n = 0; n < 1000; ++n)
+ {
+ auto r = u3(rng);
+
+ VERIFY (r[0] != 0.0 || r[1] != 0.0 || r[2] != 0.0);
+ }
+
+ __gnu_cxx::uniform_on_sphere_distribution<2> u2;
+
+ for (size_t n = 0; n < 1000; ++n)
+ {
+ auto r = u2(rng);
+
+ VERIFY (r[0] != 0.0 || r[1] != 0.0);
+ }
+}
+
+int
+main()
+{
+ test01();
+ return 0;
+}
More information about the Libstdc++
mailing list