Re: [PATCH] Fix Bug 83237 - Values returned by std::poisson_distribution are not distributed correctly

[Paolo Carlini]

Hi,

On 23/12/2017 23:10, Michele Pezzutti wrote:

I got confirmation from Luc.

He also added it to the errata file---the entries regarding p. 511, page 6 of http://luc.devroye.org/errors.pdf

Nice. Then I'm going to commit the below, tested x86_64-linux. Thanks again!

Cheers,
Paolo.

//////////////////

2017-12-24  Michele Pezzutti <mp...@tiscali.it>

        PR libstdc++/83237
        * include/bits/random.tcc (poisson_distribution<>::operator()):
        Fix __x = 1 case - see updated Errata of Devroye's treatise.
        * testsuite/26_numerics/random/poisson_distribution/operators/
        values.cc: Add test.
        * testsuite/26_numerics/random/pr60037-neg.cc: Adjust dg-error
        line number.

Index: include/bits/random.tcc
===================================================================
--- include/bits/random.tcc     (revision 255991)
+++ include/bits/random.tcc     (working copy)
@@ -1301,6 +1301,8 @@ _GLIBCXX_BEGIN_NAMESPACE_VERSION
            const double __c2 = __param._M_c2b + __c1;
            const double __c3 = __c2 + 1;
            const double __c4 = __c3 + 1;
+           // 1 / 78
+           const double __178 = 0.0128205128205128205128205128205128L;
            // e^(1 / 78)
            const double __e178 = 1.0129030479320018583185514777512983L;
            const double __c5 = __c4 + __e178;
@@ -1340,7 +1342,11 @@ _GLIBCXX_BEGIN_NAMESPACE_VERSION
                else if (__u <= __c4)
                  __x = 0;
                else if (__u <= __c5)
-                 __x = 1;
+                 {
+                   __x = 1;
+                   // Only in the Errata, see libstdc++/83237.
+                   __w = __178;
+                 }
                else
                  {
                    const double __v = -std::log(1.0 - __aurng());
Index: testsuite/26_numerics/random/poisson_distribution/operators/values.cc
===================================================================
--- testsuite/26_numerics/random/poisson_distribution/operators/values.cc       
(revision 255991)
+++ testsuite/26_numerics/random/poisson_distribution/operators/values.cc       
(working copy)
@@ -42,6 +42,12 @@ void test01()
   std::poisson_distribution<> pd3(30.0);
   auto bpd3 = std::bind(pd3, eng);
   testDiscreteDist(bpd3, [](int n) { return poisson_pdf(n, 30.0); } );
+
+  // libstdc++/83237
+  std::poisson_distribution<> pd4(37.17);
+  auto bpd4 = std::bind(pd4, eng);
+  testDiscreteDist<100, 2000000>(bpd4, [](int n)
+                                { return poisson_pdf(n, 37.17); } );
 }
 
 int main()
Index: testsuite/26_numerics/random/pr60037-neg.cc
===================================================================
--- testsuite/26_numerics/random/pr60037-neg.cc (revision 255991)
+++ testsuite/26_numerics/random/pr60037-neg.cc (working copy)
@@ -11,4 +11,4 @@ auto x = std::generate_canonical<std::size_t,
 
 // { dg-error "static assertion failed: template argument must be a floating 
point type" "" { target *-*-* } 156 }
 
-// { dg-error "static assertion failed: template argument must be a floating 
point type" "" { target *-*-* } 3311 }
+// { dg-error "static assertion failed: template argument must be a floating 
point type" "" { target *-*-* } 3317 }