From patchwork Fri Mar 25 10:10:04 2011 Content-Type: text/plain; charset="utf-8" MIME-Version: 1.0 Content-Transfer-Encoding: 7bit X-Patchwork-Submitter: Paolo Carlini X-Patchwork-Id: 88358 Return-Path: X-Original-To: incoming@patchwork.ozlabs.org Delivered-To: patchwork-incoming@bilbo.ozlabs.org Received: from sourceware.org (server1.sourceware.org [209.132.180.131]) by ozlabs.org (Postfix) with SMTP id 5C6F81007D5 for ; Fri, 25 Mar 2011 21:10:21 +1100 (EST) Received: (qmail 6086 invoked by alias); 25 Mar 2011 10:10:18 -0000 Received: (qmail 6068 invoked by uid 22791); 25 Mar 2011 10:10:17 -0000 X-SWARE-Spam-Status: No, hits=-2.3 required=5.0 tests=AWL, BAYES_00, DKIM_SIGNED, DKIM_VALID, DKIM_VALID_AU, FREEMAIL_FROM, RCVD_IN_DNSWL_LOW X-Spam-Check-By: sourceware.org Received: from mail-fx0-f47.google.com (HELO mail-fx0-f47.google.com) (209.85.161.47) by sourceware.org (qpsmtpd/0.43rc1) with ESMTP; Fri, 25 Mar 2011 10:10:12 +0000 Received: by fxm19 with SMTP id 19so1096101fxm.20 for ; Fri, 25 Mar 2011 03:10:10 -0700 (PDT) Received: by 10.223.26.205 with SMTP id f13mr667847fac.49.1301047810515; Fri, 25 Mar 2011 03:10:10 -0700 (PDT) Received: from macbook-pro-di-paolo-carlini.local (47.Red-88-2-96.staticIP.rima-tde.net [88.2.96.47]) by mx.google.com with ESMTPS id 21sm326812fav.41.2011.03.25.03.10.06 (version=SSLv3 cipher=OTHER); Fri, 25 Mar 2011 03:10:08 -0700 (PDT) Message-ID: <4D8C69FC.30902@oracle.com> Date: Fri, 25 Mar 2011 11:10:04 +0100 From: Paolo Carlini User-Agent: Mozilla/5.0 (Macintosh; U; Intel Mac OS X 10.6; it; rv:1.9.2.15) Gecko/20110303 Thunderbird/3.1.9 MIME-Version: 1.0 To: Paolo Carlini CC: Gcc Patch List , libstdc++ Subject: Re: [v3] Fix negative_binomial_distribution References: <4D8B7878.70305@oracle.com> In-Reply-To: <4D8B7878.70305@oracle.com> Mailing-List: contact gcc-patches-help@gcc.gnu.org; run by ezmlm Precedence: bulk List-Id: List-Unsubscribe: List-Archive: List-Post: List-Help: Sender: gcc-patches-owner@gcc.gnu.org Delivered-To: mailing list gcc-patches@gcc.gnu.org ... tweaking the fix like this makes for slightly faster repeated calls (spares a division) and also makes clearer that we had a plain typo p / (1 - p) for (1 - p) / p, grrr.. Double checked Devroye in the meanwhile. Committed to mainline, will be in 4.6.1. Paolo. /////////////// 2011-03-25 Paolo Carlini * include/bits/random.h (negative_binomial_distribution<>:: negative_binomial_distribution(_IntType, double), negative_binomial_distribution<>:: negative_binomial_distribution(const param_type&)): Tweak construction of _M_gd. * include/bits/random.tcc (negative_binomial_distribution<>:: operator()): Adjust. Index: include/bits/random.tcc =================================================================== --- include/bits/random.tcc (revision 171411) +++ include/bits/random.tcc (working copy) @@ -1075,7 +1075,7 @@ return __is; } - // This is Leger's algorithm. + // This is Leger's algorithm, also in Devroye, Ch. X, Example 1.5. template template typename negative_binomial_distribution<_IntType>::result_type @@ -1085,8 +1085,7 @@ const double __y = _M_gd(__urng); // XXX Is the constructor too slow? - std::poisson_distribution __poisson(__y * (1.0 - p()) - / p()); + std::poisson_distribution __poisson(__y); return __poisson(__urng); } @@ -1100,10 +1099,10 @@ typedef typename std::gamma_distribution::param_type param_type; - const double __y = _M_gd(__urng, param_type(__p.k(), 1.0)); + const double __y = + _M_gd(__urng, param_type(__p.k(), (1.0 - __p.p()) / __p.p())); - std::poisson_distribution __poisson(__y * (1.0 - __p.p()) - / __p.p() ); + std::poisson_distribution __poisson(__y); return __poisson(__urng); } Index: include/bits/random.h =================================================================== --- include/bits/random.h (revision 171411) +++ include/bits/random.h (working copy) @@ -3804,12 +3804,12 @@ explicit negative_binomial_distribution(_IntType __k = 1, double __p = 0.5) - : _M_param(__k, __p), _M_gd(__k, 1.0) + : _M_param(__k, __p), _M_gd(__k, (1.0 - __p) / __p) { } explicit negative_binomial_distribution(const param_type& __p) - : _M_param(__p), _M_gd(__p.k(), 1.0) + : _M_param(__p), _M_gd(__p.k(), (1.0 - __p.p()) / __p.p()) { } /**