re PR libstdc++/51795 (linear_congruential_engine doesn't work correctly)

2012-04-29  Marc Glisse  <marc.glisse@inria.fr>
	    Paolo Carlini  <paolo.carlini@oracle.com>

	PR libstdc++/51795
	* include/bits/stl_algobase.h (__lg<>(_Size)): Remove.
	(__lg(int), __lg(unsigned), __lg(long), __lg(unsigned long),
	__lg(long long), __lg(unsigned long long)): Define constexpr.
	* include/bits/random.h (_Mod<>): Overcome Schrage's algorithm
	limitations.
	(__mod): Adjust.
	(linear_congruential): Remove FIXME static_assert.
	* include/bits/random.tcc (_Mod<>): Adjust.
	* testsuite/26_numerics/random/linear_congruential_engine/operators/
	51795.cc: New.

Co-Authored-By: Paolo Carlini <paolo.carlini@oracle.com>

From-SVN: r186948

libstdc++-v3/ChangeLog

+2012-04-29  Marc Glisse  <marc.glisse@inria.fr>
+	    Paolo Carlini  <paolo.carlini@oracle.com>
+
+	PR libstdc++/51795
+	* include/bits/stl_algobase.h (__lg<>(_Size)): Remove.
+	(__lg(int), __lg(unsigned), __lg(long), __lg(unsigned long),
+	__lg(long long), __lg(unsigned long long)): Define constexpr.
+	* include/bits/random.h (_Mod<>): Overcome Schrage's algorithm
+	limitations.
+	(__mod): Adjust.
+	(linear_congruential): Remove FIXME static_assert.
+	* include/bits/random.tcc (_Mod<>): Adjust.
+	* testsuite/26_numerics/random/linear_congruential_engine/operators/
+	51795.cc: New.
+
 2012-04-29  Jonathan Wakely  <jwakely.gcc@gmail.com>
 
 	* include/std/functional (function::function(F)): LWG 2132: Disable

libstdc++-v3/include/bits/random.h

@@ -76,15 +76,78 @@ _GLIBCXX_END_NAMESPACE_VERSION
       struct _Shift<_UIntType, __w, true>
       { static const _UIntType __value = _UIntType(1) << __w; };
 
-    template<typename _Tp, _Tp __m, _Tp __a, _Tp __c, bool>
-      struct _Mod;
+    template<int __s,
+	     int __which = ((__s <= __CHAR_BIT__ * sizeof (int))
+			    + (__s <= __CHAR_BIT__ * sizeof (long))
+			    + (__s <= __CHAR_BIT__ * sizeof (long long))
+			    /* assume long long no bigger than __int128 */
+			    + (__s <= 128))>
+      struct _Select_uint_least_t
+      {
+	static_assert(__which < 0, /* needs to be dependent */
+		      "sorry, would be too much trouble for a slow result");
+      };
+
+    template<int __s>
+      struct _Select_uint_least_t<__s, 4>
+      { typedef unsigned int type; };
+
+    template<int __s>
+      struct _Select_uint_least_t<__s, 3>
+      { typedef unsigned long type; };
+
+    template<int __s>
+      struct _Select_uint_least_t<__s, 2>
+      { typedef unsigned long long type; };
+
+#ifdef _GLIBCXX_USE_INT128
+    template<int __s>
+      struct _Select_uint_least_t<__s, 1>
+      { typedef unsigned __int128 type; };
+#endif
+
+    // Assume a != 0, a < m, c < m, x < m.
+    template<typename _Tp, _Tp __m, _Tp __a, _Tp __c,
+	     bool __big_enough = (!(__m & (__m - 1))
+				  || (_Tp(-1) - __c) / __a >= __m - 1),
+             bool __schrage_ok = __m % __a < __m / __a>
+      struct _Mod
+      {
+	typedef typename _Select_uint_least_t<std::__lg(__a)
+					      + std::__lg(__m) + 2>::type _Tp2;
+	static _Tp
+	__calc(_Tp __x)
+	{ return static_cast<_Tp>((_Tp2(__a) * __x + __c) % __m); }
+      };
+
+    // Schrage.
+    template<typename _Tp, _Tp __m, _Tp __a, _Tp __c>
+      struct _Mod<_Tp, __m, __a, __c, false, true>
+      {
+	static _Tp
+	__calc(_Tp __x);
+      };
+
+    // Special cases:
+    // - for m == 2^n or m == 0, unsigned integer overflow is safe.
+    // - a * (m - 1) + c fits in _Tp, there is no overflow.
+    template<typename _Tp, _Tp __m, _Tp __a, _Tp __c, bool __s>
+      struct _Mod<_Tp, __m, __a, __c, true, __s>
+      {
+	static _Tp
+	__calc(_Tp __x)
+	{
+	  _Tp __res = __a * __x + __c;
+	  if (__m)
+	    __res %= __m;
+	  return __res;
+	}
+      };
 
-    // Dispatch based on modulus value to prevent divide-by-zero compile-time
-    // errors when m == 0.
     template<typename _Tp, _Tp __m, _Tp __a = 1, _Tp __c = 0>
       inline _Tp
       __mod(_Tp __x)
-      { return _Mod<_Tp, __m, __a, __c, __m == 0>::__calc(__x); }
+      { return _Mod<_Tp, __m, __a, __c>::__calc(__x); }
 
     /*
      * An adaptor class for converting the output of any Generator into
@@ -174,11 +237,6 @@ _GLIBCXX_BEGIN_NAMESPACE_VERSION
       static_assert(__m == 0u || (__a < __m && __c < __m),
 		    "template argument substituting __m out of bounds");
 
-      // XXX FIXME:
-      // _Mod::__calc should handle correctly __m % __a >= __m / __a too.
-      static_assert(__m % __a < __m / __a,
-		    "sorry, not implemented yet: try a smaller 'a' constant");
-
     public:
       /** The type of the generated random value. */
       typedef _UIntType result_type;

libstdc++-v3/include/bits/random.tcc

@@ -41,58 +41,42 @@ namespace std _GLIBCXX_VISIBILITY(default)
   {
   _GLIBCXX_BEGIN_NAMESPACE_VERSION
 
-    // General case for x = (ax + c) mod m -- use Schrage's algorithm to
-    // avoid integer overflow.
-    //
-    // Because a and c are compile-time integral constants the compiler
-    // kindly elides any unreachable paths.
+    // General case for x = (ax + c) mod m -- use Schrage's algorithm
+    // to avoid integer overflow.
     //
     // Preconditions:  a > 0, m > 0.
     //
-    // XXX FIXME: as-is, only works correctly for __m % __a < __m / __a. 
-    //
-    template<typename _Tp, _Tp __m, _Tp __a, _Tp __c, bool>
-      struct _Mod
-      {
-	static _Tp
-	__calc(_Tp __x)
-	{
-	  if (__a == 1)
-	    __x %= __m;
-	  else
-	    {
-	      static const _Tp __q = __m / __a;
-	      static const _Tp __r = __m % __a;
-
-	      _Tp __t1 = __a * (__x % __q);
-	      _Tp __t2 = __r * (__x / __q);
-	      if (__t1 >= __t2)
-		__x = __t1 - __t2;
-	      else
-		__x = __m - __t2 + __t1;
-	    }
-
-	  if (__c != 0)
-	    {
-	      const _Tp __d = __m - __x;
-	      if (__d > __c)
-		__x += __c;
-	      else
-		__x = __c - __d;
-	    }
-	  return __x;
-	}
-      };
-
-    // Special case for m == 0 -- use unsigned integer overflow as modulo
-    // operator.
+    // Note: only works correctly for __m % __a < __m / __a.
     template<typename _Tp, _Tp __m, _Tp __a, _Tp __c>
-      struct _Mod<_Tp, __m, __a, __c, true>
+      _Tp
+      _Mod<_Tp, __m, __a, __c, false, true>::
+      __calc(_Tp __x)
       {
-	static _Tp
-	__calc(_Tp __x)
-	{ return __a * __x + __c; }
-      };
+	if (__a == 1)
+	  __x %= __m;
+	else
+	  {
+	    static const _Tp __q = __m / __a;
+	    static const _Tp __r = __m % __a;
+
+	    _Tp __t1 = __a * (__x % __q);
+	    _Tp __t2 = __r * (__x / __q);
+	    if (__t1 >= __t2)
+	      __x = __t1 - __t2;
+	    else
+	      __x = __m - __t2 + __t1;
+	  }
+
+	if (__c != 0)
+	  {
+	    const _Tp __d = __m - __x;
+	    if (__d > __c)
+	      __x += __c;
+	    else
+	      __x = __c - __d;
+	  }
+	return __x;
+      }
 
     template<typename _InputIterator, typename _OutputIterator,
 	     typename _UnaryOperation>

libstdc++-v3/include/bits/stl_algobase.h

@@ -975,37 +975,27 @@ _GLIBCXX_BEGIN_NAMESPACE_VERSION
 
   /// This is a helper function for the sort routines and for random.tcc.
   //  Precondition: __n > 0.
-  template<typename _Size>
-    inline _Size
-    __lg(_Size __n)
-    {
-      _Size __k;
-      for (__k = 0; __n != 0; __n >>= 1)
-	++__k;
-      return __k - 1;
-    }
-
-  inline int
+  inline _GLIBCXX_CONSTEXPR int
   __lg(int __n)
   { return sizeof(int) * __CHAR_BIT__  - 1 - __builtin_clz(__n); }
 
-  inline unsigned
+  inline _GLIBCXX_CONSTEXPR unsigned
   __lg(unsigned __n)
   { return sizeof(int) * __CHAR_BIT__  - 1 - __builtin_clz(__n); }
 
-  inline long
+  inline _GLIBCXX_CONSTEXPR long
   __lg(long __n)
   { return sizeof(long) * __CHAR_BIT__ - 1 - __builtin_clzl(__n); }
 
-  inline unsigned long
+  inline _GLIBCXX_CONSTEXPR unsigned long
   __lg(unsigned long __n)
   { return sizeof(long) * __CHAR_BIT__ - 1 - __builtin_clzl(__n); }
 
-  inline long long
+  inline _GLIBCXX_CONSTEXPR long long
   __lg(long long __n)
   { return sizeof(long long) * __CHAR_BIT__ - 1 - __builtin_clzll(__n); }
 
-  inline unsigned long long
+  inline _GLIBCXX_CONSTEXPR unsigned long long
   __lg(unsigned long long __n)
   { return sizeof(long long) * __CHAR_BIT__ - 1 - __builtin_clzll(__n); }
 

libstdc++-v3/testsuite/26_numerics/random/linear_congruential_engine/operators/51795.cc

+// { dg-options "-std=gnu++11" }
+// { dg-require-cstdint "" }
+//
+// Copyright (C) 2012 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/>.
+
+// 26.5.3.1 class template linear_congruential_engine [rand.eng.lcong]
+
+#include <random>
+#include <testsuite_hooks.h>
+
+void test01()
+{
+  bool test __attribute__((unused)) = true;
+
+  typedef std::linear_congruential_engine<std::uint64_t, 1103515245ULL,
+					  12345, 2147483648ULL> engine;
+  engine eng(1103527590ULL);
+
+  for (unsigned it = 0; it < 1000; ++it)
+    {
+      std::uint64_t num = eng();
+      VERIFY( (num >= eng.min() && num <= eng.max()) );
+    }
+}
+
+int main()
+{
+  test01();
+  return 0;
+}