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

Paolo Carlini paolo.carlini@oracle.com
Sun Dec 24 15:54:00 GMT 2017


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.

//////////////////
-------------- next part --------------
2017-12-24  Michele Pezzutti <mpezz@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.
-------------- next part --------------
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 }


More information about the Libstdc++ mailing list