| Message ID | 878umxtxei.fsf@x240.local.i-did-not-set--mail-host-address--so-tickle-me |
|---|---|
| State | New |
| Headers | show |
On Sat, 9 Aug 2014, Ulrich Drepper wrote: > How about the patch below? Looks good, with two details: > + 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(); I must be missing something obvious, but normal_distribution uses: __x = result_type(2.0) * __aurng() - 1.0; to get a number between -1 and 1, and I don't see where you do this rescaling. Does __aurng() already return a number in the right interval in this context? If so we may want to update our naming of variables to make that clearer. > + __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]); Assuming the 2 coordinates are obtained through a rescaling x->2*x-1, if __sq is not exactly 0, it must be between 2^-103 and 1 (for ieee double), so I am not sure hypot gains that much (at least in my mind hypot was mostly a gain close to 0 or infinity, but maybe it has more advantages). It can only hurt speed though, so not a big issue.
On Sat, Aug 9, 2014 at 1:40 PM, Marc Glisse <marc.glisse@inria.fr> wrote: > __x = result_type(2.0) * __aurng() - 1.0; You're right, we of course need the negatives as well. > Assuming the 2 coordinates are obtained through a rescaling x->2*x-1, if > __sq is not exactly 0, it must be between 2^-103 and 1 (for ieee > double), so I am not sure hypot gains that much (at least in my mind > hypot was mostly a gain close to 0 or infinity, but maybe it has more > advantages). It can only hurt speed though, so not a big issue. Depending on how similar in size the two values are, not using hypot() can drop quite a few bits. Especially with the scaling through division this error can be noticeable. Better be sure. Maybe at some point I have time to investigate the worst case scenario for the numbers in question but until this shows hypot isn't needed it's best to leave it in. I've committed the patch.
On Sun, Aug 10, 2014 at 2:00 AM, Ulrich Drepper <drepper@gmail.com> wrote: > On Sat, Aug 9, 2014 at 1:40 PM, Marc Glisse <marc.glisse@inria.fr> wrote: >> __x = result_type(2.0) * __aurng() - 1.0; > > You're right, we of course need the negatives as well. > >> Assuming the 2 coordinates are obtained through a rescaling x->2*x-1, if >> __sq is not exactly 0, it must be between 2^-103 and 1 (for ieee >> double), so I am not sure hypot gains that much (at least in my mind >> hypot was mostly a gain close to 0 or infinity, but maybe it has more >> advantages). It can only hurt speed though, so not a big issue. > > Depending on how similar in size the two values are, not using hypot() > can drop quite a few bits. Especially with the scaling through > division this error can be noticeable. Better be sure. Maybe at some > point I have time to investigate the worst case scenario for the > numbers in question but until this shows hypot isn't needed it's best > to leave it in. > > I've committed the patch. Hi, This causes lots of tests under libstdc++-v3/testsuite/ext/ failed on aarch64-none-elf and arm-none-eabi with below log message. In file included from /home/binche01/work/oban-work/build-aarch64-none-elf/obj/gcc2/aarch64-none-elf/libstdc++-v3/include/ext/random:3494:0, from /home/binche01/work/oban-work/src/gcc/libstdc++-v3/testsuite/ext/random/arcsine_distribution/cons/default.cc:23: /home/binche01/work/oban-work/src/gcc/libstdc++-v3/include/ext/random.tcc: In member function '__gnu_cxx::{anonymous}::uniform_on_sphere_helper<2ul, _RealType>::result_type __gnu_cxx::{anonymous}::uniform_on_sphere_helper<2ul, _RealType>::operator()(_NormalDistribution&, _UniformRandomNumberGenerator&)': /home/binche01/work/oban-work/src/gcc/libstdc++-v3/include/ext/random.tcc:1609:18: error: 'hypot' is not a member of 'std' auto __norm = std::hypot(__ret[0], __ret[1]); ^ /home/binche01/work/oban-work/src/gcc/libstdc++-v3/include/ext/random.tcc:1609:18: note: suggested alternative: In file included from /home/binche01/work/oban-work/build-aarch64-none-elf/obj/gcc2/aarch64-none-elf/libstdc++-v3/include/cmath:44:0, from /home/binche01/work/oban-work/build-aarch64-none-elf/obj/gcc2/aarch64-none-elf/libstdc++-v3/include/random:38, from /home/binche01/work/oban-work/build-aarch64-none-elf/obj/gcc2/aarch64-none-elf/libstdc++-v3/include/ext/random:38, from /home/binche01/work/oban-work/src/gcc/libstdc++-v3/testsuite/ext/random/arcsine_distribution/cons/default.cc:23: /home/binche01/work/oban-work/target-aarch64-none-elf/aarch64-none-elf/include/math.h:306:15: note: 'hypot' extern double hypot _PARAMS((double, double)); ^ FAIL: ext/random/arcsine_distribution/cons/default.cc (test for excess errors) Excess errors: /home/binche01/work/oban-work/src/gcc/libstdc++-v3/include/ext/random.tcc:1609:18: error: 'hypot' is not a member of 'std' Thanks, bin
On Wed, Aug 13, 2014 at 3:36 PM, Bin.Cheng <amker.cheng@gmail.com> wrote: > On Sun, Aug 10, 2014 at 2:00 AM, Ulrich Drepper <drepper@gmail.com> wrote: >> On Sat, Aug 9, 2014 at 1:40 PM, Marc Glisse <marc.glisse@inria.fr> wrote: >>> __x = result_type(2.0) * __aurng() - 1.0; >> >> You're right, we of course need the negatives as well. >> >>> Assuming the 2 coordinates are obtained through a rescaling x->2*x-1, if >>> __sq is not exactly 0, it must be between 2^-103 and 1 (for ieee >>> double), so I am not sure hypot gains that much (at least in my mind >>> hypot was mostly a gain close to 0 or infinity, but maybe it has more >>> advantages). It can only hurt speed though, so not a big issue. >> >> Depending on how similar in size the two values are, not using hypot() >> can drop quite a few bits. Especially with the scaling through >> division this error can be noticeable. Better be sure. Maybe at some >> point I have time to investigate the worst case scenario for the >> numbers in question but until this shows hypot isn't needed it's best >> to leave it in. >> >> I've committed the patch. > > Hi, > > This causes lots of tests under libstdc++-v3/testsuite/ext/ failed on > aarch64-none-elf and arm-none-eabi with below log message. > > In file included from > /home/binche01/work/oban-work/build-aarch64-none-elf/obj/gcc2/aarch64-none-elf/libstdc++-v3/include/ext/random:3494:0, > from > /home/binche01/work/oban-work/src/gcc/libstdc++-v3/testsuite/ext/random/arcsine_distribution/cons/default.cc:23: > /home/binche01/work/oban-work/src/gcc/libstdc++-v3/include/ext/random.tcc: > In member function > '__gnu_cxx::{anonymous}::uniform_on_sphere_helper<2ul, > _RealType>::result_type > __gnu_cxx::{anonymous}::uniform_on_sphere_helper<2ul, > _RealType>::operator()(_NormalDistribution&, > _UniformRandomNumberGenerator&)': > /home/binche01/work/oban-work/src/gcc/libstdc++-v3/include/ext/random.tcc:1609:18: > error: 'hypot' is not a member of 'std' > auto __norm = std::hypot(__ret[0], __ret[1]); > ^ > /home/binche01/work/oban-work/src/gcc/libstdc++-v3/include/ext/random.tcc:1609:18: > note: suggested alternative: > In file included from > /home/binche01/work/oban-work/build-aarch64-none-elf/obj/gcc2/aarch64-none-elf/libstdc++-v3/include/cmath:44:0, > from > /home/binche01/work/oban-work/build-aarch64-none-elf/obj/gcc2/aarch64-none-elf/libstdc++-v3/include/random:38, > from > /home/binche01/work/oban-work/build-aarch64-none-elf/obj/gcc2/aarch64-none-elf/libstdc++-v3/include/ext/random:38, > from > /home/binche01/work/oban-work/src/gcc/libstdc++-v3/testsuite/ext/random/arcsine_distribution/cons/default.cc:23: > /home/binche01/work/oban-work/target-aarch64-none-elf/aarch64-none-elf/include/math.h:306:15: > note: 'hypot' > extern double hypot _PARAMS((double, double)); > ^ > > FAIL: ext/random/arcsine_distribution/cons/default.cc (test for excess errors) > Excess errors: > /home/binche01/work/oban-work/src/gcc/libstdc++-v3/include/ext/random.tcc:1609:18: > error: 'hypot' is not a member of 'std' > https://gcc.gnu.org/bugzilla/show_bug.cgi?id=62118 filed. Thanks, bin
Ulrich-- On 08/13/2014 09:36 AM, Bin.Cheng wrote: > On Sun, Aug 10, 2014 at 2:00 AM, Ulrich Drepper <drepper@gmail.com> wrote: >> On Sat, Aug 9, 2014 at 1:40 PM, Marc Glisse <marc.glisse@inria.fr> wrote: >>> __x = result_type(2.0) * __aurng() - 1.0; >> You're right, we of course need the negatives as well. >> >>> Assuming the 2 coordinates are obtained through a rescaling x->2*x-1, if >>> __sq is not exactly 0, it must be between 2^-103 and 1 (for ieee >>> double), so I am not sure hypot gains that much (at least in my mind >>> hypot was mostly a gain close to 0 or infinity, but maybe it has more >>> advantages). It can only hurt speed though, so not a big issue. >> Depending on how similar in size the two values are, not using hypot() >> can drop quite a few bits. Especially with the scaling through >> division this error can be noticeable. Better be sure. Maybe at some >> point I have time to investigate the worst case scenario for the >> numbers in question but until this shows hypot isn't needed it's best >> to leave it in. >> >> I've committed the patch. > Hi, > > This causes lots of tests under libstdc++-v3/testsuite/ext/ failed on > aarch64-none-elf and arm-none-eabi with below log message. Please follow our usual rules, don't unconditionally use C99 functions like std::hypot, use _GLIBCXX_USE_C99_MATH_TR1 (many examples, eg in random.tcc). Thanks! Paolo.
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; +}