From mboxrd@z Thu Jan 1 00:00:00 1970 Return-Path: Received: (qmail 3849 invoked by alias); 9 Aug 2014 16:56:11 -0000 Mailing-List: contact gcc-patches-help@gcc.gnu.org; run by ezmlm Precedence: bulk List-Id: List-Archive: List-Post: List-Help: Sender: gcc-patches-owner@gcc.gnu.org Received: (qmail 3830 invoked by uid 89); 9 Aug 2014 16:56:11 -0000 Authentication-Results: sourceware.org; auth=none X-Virus-Found: No X-Spam-SWARE-Status: No, score=-0.4 required=5.0 tests=AWL,BAYES_00,FREEMAIL_FROM,RCVD_IN_DNSWL_LOW,SPF_PASS autolearn=ham version=3.3.2 X-Spam-User: qpsmtpd, 2 recipients X-HELO: mail-qa0-f49.google.com Received: from mail-qa0-f49.google.com (HELO mail-qa0-f49.google.com) (209.85.216.49) by sourceware.org (qpsmtpd/0.93/v0.84-503-g423c35a) with (AES128-SHA encrypted) ESMTPS; Sat, 09 Aug 2014 16:56:09 +0000 Received: by mail-qa0-f49.google.com with SMTP id dc16so6615901qab.36 for ; Sat, 09 Aug 2014 09:56:06 -0700 (PDT) X-Received: by 10.224.95.6 with SMTP id b6mr49223248qan.17.1407603366867; Sat, 09 Aug 2014 09:56:06 -0700 (PDT) Received: from x240.local (cpe-68-173-15-30.nyc.res.rr.com. [68.173.15.30]) by mx.google.com with ESMTPSA id a8sm11689652qat.28.2014.08.09.09.56.06 for (version=TLSv1.2 cipher=ECDHE-RSA-AES128-GCM-SHA256 bits=128/128); Sat, 09 Aug 2014 09:56:06 -0700 (PDT) From: Ulrich Drepper X-Google-Original-From: Ulrich Drepper Received: from x240.local (localhost.localdomain [127.0.0.1]) by x240.local (8.14.8/8.14.8) with ESMTP id s79Gu5iG028145; Sat, 9 Aug 2014 12:56:05 -0400 Received: (from drepper@localhost) by x240.local (8.14.8/8.14.8/Submit) id s79Gu59G028144; Sat, 9 Aug 2014 12:56:05 -0400 To: libstdc++@gcc.gnu.org Cc: Jonathan Wakely , GCC Patches Subject: Re: [PATCH] libstdc++: add uniform on sphere distribution References: <87a98eow4m.fsf@x240.local.i-did-not-set--mail-host-address--so-tickle-me> <20140723102908.GM2361@redhat.com> <87d2catvi3.fsf@x240.local.i-did-not-set--mail-host-address--so-tickle-me> Date: Sat, 09 Aug 2014 16:56:00 -0000 In-Reply-To: (Marc Glisse's message of "Sat, 9 Aug 2014 17:33:43 +0200 (CEST)") Message-ID: <878umxtxei.fsf@x240.local.i-did-not-set--mail-host-address--so-tickle-me> User-Agent: Gnus/5.13 (Gnus v5.13) Emacs/24.3 (gnu/linux) MIME-Version: 1.0 Content-Type: text/plain X-SW-Source: 2014-08/txt/msg00956.txt.bz2 Marc Glisse 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 * 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 + class uniform_on_sphere_helper + { + typedef typename uniform_on_sphere_distribution<_Dimen, _RealType>::result_type result_type; + + public: + template + 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 + class uniform_on_sphere_helper<2, _RealType> + { + typedef typename uniform_on_sphere_distribution<2, _RealType>:: + result_type result_type; + + public: + template + 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 template 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 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 // . -// 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 #include 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 // . -// 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 #include 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 // . -// 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 #include @@ -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 +// +// 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 +// . + +// Class template uniform_on_sphere_distribution [rand.dist.ext.uniform_on_sphere] +// Concept RandomNumberDistribution [rand.concept.dist] + +#include +#include +#include + +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; +}