random.h (negative_binomial_distribution<>:: negative_binomial_distribution(_IntType...
[gcc.git] / libstdc++-v3 / include / bits / random.tcc
index c8b6590a28d64ae9f2834f49eb7111ea45fae580..01ee5824cd9f7df64c564e0910d285f749c53f73 100644 (file)
@@ -1,6 +1,6 @@
 // random number generation (out of line) -*- C++ -*-
 
-// Copyright (C) 2009 Free Software Foundation, Inc.
+// Copyright (C) 2009, 2010, 2011 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
 
 /** @file bits/random.tcc
  *  This is an internal header file, included by other library headers.
- *  You should not attempt to use it directly.
+ *  Do not attempt to use it directly. @headername{random}
  */
 
-#include <numeric>
-#include <algorithm>
+#ifndef _RANDOM_TCC
+#define _RANDOM_TCC 1
 
-namespace std
+#include <numeric> // std::accumulate and std::partial_sum
+
+namespace std _GLIBCXX_VISIBILITY(default)
 {
   /*
    * (Further) implementation-space details.
    */
   namespace __detail
   {
+  _GLIBCXX_BEGIN_NAMESPACE_VERSION
+
     // General case for x = (ax + c) mod m -- use Schrage's algorithm to
     // avoid integer overflow.
     //
@@ -87,8 +91,39 @@ namespace std
        __calc(_Tp __x)
        { return __a * __x + __c; }
       };
+
+    template<typename _InputIterator, typename _OutputIterator,
+            typename _UnaryOperation>
+      _OutputIterator
+      __transform(_InputIterator __first, _InputIterator __last,
+                 _OutputIterator __result, _UnaryOperation __unary_op)
+      {
+       for (; __first != __last; ++__first, ++__result)
+         *__result = __unary_op(*__first);
+       return __result;
+      }
+
+  _GLIBCXX_END_NAMESPACE_VERSION
   } // namespace __detail
 
+_GLIBCXX_BEGIN_NAMESPACE_VERSION
+
+  template<typename _UIntType, _UIntType __a, _UIntType __c, _UIntType __m>
+    constexpr _UIntType
+    linear_congruential_engine<_UIntType, __a, __c, __m>::multiplier;
+
+  template<typename _UIntType, _UIntType __a, _UIntType __c, _UIntType __m>
+    constexpr _UIntType
+    linear_congruential_engine<_UIntType, __a, __c, __m>::increment;
+
+  template<typename _UIntType, _UIntType __a, _UIntType __c, _UIntType __m>
+    constexpr _UIntType
+    linear_congruential_engine<_UIntType, __a, __c, __m>::modulus;
+
+  template<typename _UIntType, _UIntType __a, _UIntType __c, _UIntType __m>
+    constexpr _UIntType
+    linear_congruential_engine<_UIntType, __a, __c, __m>::default_seed;
+
   /**
    * Seeds the LCR with integral value @p __s, adjusted so that the
    * ring identity is never a member of the convergence set.
@@ -109,24 +144,25 @@ namespace std
    * Seeds the LCR engine with a value generated by @p __q.
    */
   template<typename _UIntType, _UIntType __a, _UIntType __c, _UIntType __m>
-    void
-    linear_congruential_engine<_UIntType, __a, __c, __m>::
-    seed(seed_seq& __q)
-    {
-      const _UIntType __k0 = __m == 0 ? std::numeric_limits<_UIntType>::digits
-                                     : std::__lg(__m);
-      const _UIntType __k = (__k0 + 31) / 32;
-      _UIntType __arr[__k + 3];
-      __q.generate(__arr + 0, __arr + __k + 3);
-      _UIntType __factor = 1u;
-      _UIntType __sum = 0u;
-      for (size_t __j = 0; __j < __k; ++__j)
-        {
-          __sum += __arr[__j + 3] * __factor;
-          __factor *= __detail::_Shift<_UIntType, 32>::__value;
-        }
-      seed(__sum);
-    }
+    template<typename _Sseq>
+      typename std::enable_if<std::is_class<_Sseq>::value>::type
+      linear_congruential_engine<_UIntType, __a, __c, __m>::
+      seed(_Sseq& __q)
+      {
+       const _UIntType __k0 = __m == 0 ? std::numeric_limits<_UIntType>::digits
+                                       : std::__lg(__m);
+       const _UIntType __k = (__k0 + 31) / 32;
+       uint_least32_t __arr[__k + 3];
+       __q.generate(__arr + 0, __arr + __k + 3);
+       _UIntType __factor = 1u;
+       _UIntType __sum = 0u;
+       for (size_t __j = 0; __j < __k; ++__j)
+         {
+           __sum += __arr[__j + 3] * __factor;
+           __factor *= __detail::_Shift<_UIntType, 32>::__value;
+         }
+       seed(__sum);
+      }
 
   template<typename _UIntType, _UIntType __a, _UIntType __c, _UIntType __m,
           typename _CharT, typename _Traits>
@@ -169,6 +205,133 @@ namespace std
     }
 
 
+  template<typename _UIntType,
+          size_t __w, size_t __n, size_t __m, size_t __r,
+          _UIntType __a, size_t __u, _UIntType __d, size_t __s,
+          _UIntType __b, size_t __t, _UIntType __c, size_t __l,
+          _UIntType __f>
+    constexpr size_t
+    mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d,
+                           __s, __b, __t, __c, __l, __f>::word_size;
+
+  template<typename _UIntType,
+          size_t __w, size_t __n, size_t __m, size_t __r,
+          _UIntType __a, size_t __u, _UIntType __d, size_t __s,
+          _UIntType __b, size_t __t, _UIntType __c, size_t __l,
+          _UIntType __f>
+    constexpr size_t
+    mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d,
+                           __s, __b, __t, __c, __l, __f>::state_size;
+
+  template<typename _UIntType,
+          size_t __w, size_t __n, size_t __m, size_t __r,
+          _UIntType __a, size_t __u, _UIntType __d, size_t __s,
+          _UIntType __b, size_t __t, _UIntType __c, size_t __l,
+          _UIntType __f>
+    constexpr size_t
+    mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d,
+                           __s, __b, __t, __c, __l, __f>::shift_size;
+
+  template<typename _UIntType,
+          size_t __w, size_t __n, size_t __m, size_t __r,
+          _UIntType __a, size_t __u, _UIntType __d, size_t __s,
+          _UIntType __b, size_t __t, _UIntType __c, size_t __l,
+          _UIntType __f>
+    constexpr size_t
+    mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d,
+                           __s, __b, __t, __c, __l, __f>::mask_bits;
+
+  template<typename _UIntType,
+          size_t __w, size_t __n, size_t __m, size_t __r,
+          _UIntType __a, size_t __u, _UIntType __d, size_t __s,
+          _UIntType __b, size_t __t, _UIntType __c, size_t __l,
+          _UIntType __f>
+    constexpr _UIntType
+    mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d,
+                           __s, __b, __t, __c, __l, __f>::xor_mask;
+
+  template<typename _UIntType,
+          size_t __w, size_t __n, size_t __m, size_t __r,
+          _UIntType __a, size_t __u, _UIntType __d, size_t __s,
+          _UIntType __b, size_t __t, _UIntType __c, size_t __l,
+          _UIntType __f>
+    constexpr size_t
+    mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d,
+                           __s, __b, __t, __c, __l, __f>::tempering_u;
+   
+  template<typename _UIntType,
+          size_t __w, size_t __n, size_t __m, size_t __r,
+          _UIntType __a, size_t __u, _UIntType __d, size_t __s,
+          _UIntType __b, size_t __t, _UIntType __c, size_t __l,
+          _UIntType __f>
+    constexpr _UIntType
+    mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d,
+                           __s, __b, __t, __c, __l, __f>::tempering_d;
+
+  template<typename _UIntType,
+          size_t __w, size_t __n, size_t __m, size_t __r,
+          _UIntType __a, size_t __u, _UIntType __d, size_t __s,
+          _UIntType __b, size_t __t, _UIntType __c, size_t __l,
+          _UIntType __f>
+    constexpr size_t
+    mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d,
+                           __s, __b, __t, __c, __l, __f>::tempering_s;
+
+  template<typename _UIntType,
+          size_t __w, size_t __n, size_t __m, size_t __r,
+          _UIntType __a, size_t __u, _UIntType __d, size_t __s,
+          _UIntType __b, size_t __t, _UIntType __c, size_t __l,
+          _UIntType __f>
+    constexpr _UIntType
+    mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d,
+                           __s, __b, __t, __c, __l, __f>::tempering_b;
+
+  template<typename _UIntType,
+          size_t __w, size_t __n, size_t __m, size_t __r,
+          _UIntType __a, size_t __u, _UIntType __d, size_t __s,
+          _UIntType __b, size_t __t, _UIntType __c, size_t __l,
+          _UIntType __f>
+    constexpr size_t
+    mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d,
+                           __s, __b, __t, __c, __l, __f>::tempering_t;
+
+  template<typename _UIntType,
+          size_t __w, size_t __n, size_t __m, size_t __r,
+          _UIntType __a, size_t __u, _UIntType __d, size_t __s,
+          _UIntType __b, size_t __t, _UIntType __c, size_t __l,
+          _UIntType __f>
+    constexpr _UIntType
+    mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d,
+                           __s, __b, __t, __c, __l, __f>::tempering_c;
+
+  template<typename _UIntType,
+          size_t __w, size_t __n, size_t __m, size_t __r,
+          _UIntType __a, size_t __u, _UIntType __d, size_t __s,
+          _UIntType __b, size_t __t, _UIntType __c, size_t __l,
+          _UIntType __f>
+    constexpr size_t
+    mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d,
+                           __s, __b, __t, __c, __l, __f>::tempering_l;
+
+  template<typename _UIntType,
+          size_t __w, size_t __n, size_t __m, size_t __r,
+          _UIntType __a, size_t __u, _UIntType __d, size_t __s,
+          _UIntType __b, size_t __t, _UIntType __c, size_t __l,
+          _UIntType __f>
+    constexpr _UIntType
+    mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d,
+                           __s, __b, __t, __c, __l, __f>::
+                                              initialization_multiplier;
+
+  template<typename _UIntType,
+          size_t __w, size_t __n, size_t __m, size_t __r,
+          _UIntType __a, size_t __u, _UIntType __d, size_t __s,
+          _UIntType __b, size_t __t, _UIntType __c, size_t __l,
+          _UIntType __f>
+    constexpr _UIntType
+    mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d,
+                           __s, __b, __t, __c, __l, __f>::default_seed;
+
   template<typename _UIntType,
           size_t __w, size_t __n, size_t __m, size_t __r,
           _UIntType __a, size_t __u, _UIntType __d, size_t __s,
@@ -199,43 +362,44 @@ namespace std
           _UIntType __a, size_t __u, _UIntType __d, size_t __s,
           _UIntType __b, size_t __t, _UIntType __c, size_t __l,
           _UIntType __f>
-    void
-    mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d,
+    template<typename _Sseq>
+      typename std::enable_if<std::is_class<_Sseq>::value>::type
+      mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d,
                              __s, __b, __t, __c, __l, __f>::
-    seed(seed_seq& __q)
-    {
-      const _UIntType __upper_mask = (~_UIntType()) << __r;
-      const size_t __k = (__w + 31) / 32;
-      uint_least32_t __arr[__n * __k];
-      __q.generate(__arr + 0, __arr + __n * __k);
-
-      bool __zero = true;
-      for (size_t __i = 0; __i < state_size; ++__i)
-        {
-          _UIntType __factor = 1u;
-          _UIntType __sum = 0u;
-          for (size_t __j = 0; __j < __k; ++__j)
-            {
-             __sum += __arr[__k * __i + __j] * __factor;
-             __factor *= __detail::_Shift<_UIntType, 32>::__value;
-            }
-          _M_x[__i] = __detail::__mod<_UIntType,
-           __detail::_Shift<_UIntType, __w>::__value>(__sum);
+      seed(_Sseq& __q)
+      {
+       const _UIntType __upper_mask = (~_UIntType()) << __r;
+       const size_t __k = (__w + 31) / 32;
+       uint_least32_t __arr[__n * __k];
+       __q.generate(__arr + 0, __arr + __n * __k);
 
-          if (__zero)
-            {
-             if (__i == 0)
-               {
-                 if ((_M_x[0] & __upper_mask) != 0u)
-                   __zero = false;
-               }
-             else if (_M_x[__i] != 0u)
-               __zero = false;
-            }
-        }
+       bool __zero = true;
+       for (size_t __i = 0; __i < state_size; ++__i)
+         {
+           _UIntType __factor = 1u;
+           _UIntType __sum = 0u;
+           for (size_t __j = 0; __j < __k; ++__j)
+             {
+               __sum += __arr[__k * __i + __j] * __factor;
+               __factor *= __detail::_Shift<_UIntType, 32>::__value;
+             }
+           _M_x[__i] = __detail::__mod<_UIntType,
+             __detail::_Shift<_UIntType, __w>::__value>(__sum);
+
+           if (__zero)
+             {
+               if (__i == 0)
+                 {
+                   if ((_M_x[0] & __upper_mask) != 0u)
+                     __zero = false;
+                 }
+               else if (_M_x[__i] != 0u)
+                 __zero = false;
+             }
+         }
         if (__zero)
           _M_x[0] = __detail::_Shift<_UIntType, __w - 1>::__value;
-    }
+      }
 
   template<typename _UIntType, size_t __w,
           size_t __n, size_t __m, size_t __r,
@@ -340,6 +504,22 @@ namespace std
     }
 
 
+  template<typename _UIntType, size_t __w, size_t __s, size_t __r>
+    constexpr size_t
+    subtract_with_carry_engine<_UIntType, __w, __s, __r>::word_size;
+
+  template<typename _UIntType, size_t __w, size_t __s, size_t __r>
+    constexpr size_t
+    subtract_with_carry_engine<_UIntType, __w, __s, __r>::short_lag;
+
+  template<typename _UIntType, size_t __w, size_t __s, size_t __r>
+    constexpr size_t
+    subtract_with_carry_engine<_UIntType, __w, __s, __r>::long_lag;
+
+  template<typename _UIntType, size_t __w, size_t __s, size_t __r>
+    constexpr _UIntType
+    subtract_with_carry_engine<_UIntType, __w, __s, __r>::default_seed;
+
   template<typename _UIntType, size_t __w, size_t __s, size_t __r>
     void
     subtract_with_carry_engine<_UIntType, __w, __s, __r>::
@@ -369,29 +549,30 @@ namespace std
     }
 
   template<typename _UIntType, size_t __w, size_t __s, size_t __r>
-    void
-    subtract_with_carry_engine<_UIntType, __w, __s, __r>::
-    seed(seed_seq& __q)
-    {
-      const size_t __k = (__w + 31) / 32;
-      uint_least32_t __arr[__r * __k];
-      __q.generate(__arr + 0, __arr + __r * __k);
+    template<typename _Sseq>
+      typename std::enable_if<std::is_class<_Sseq>::value>::type
+      subtract_with_carry_engine<_UIntType, __w, __s, __r>::
+      seed(_Sseq& __q)
+      {
+       const size_t __k = (__w + 31) / 32;
+       uint_least32_t __arr[__r * __k];
+       __q.generate(__arr + 0, __arr + __r * __k);
 
-      for (size_t __i = 0; __i < long_lag; ++__i)
-        {
-          _UIntType __sum = 0u;
-          _UIntType __factor = 1u;
-          for (size_t __j = 0; __j < __k; ++__j)
-            {
-             __sum += __arr[__k * __i + __j] * __factor;
-             __factor *= __detail::_Shift<_UIntType, 32>::__value;
-            }
-          _M_x[__i] = __detail::__mod<_UIntType,
-           __detail::_Shift<_UIntType, __w>::__value>(__sum);
-        }
-      _M_carry = (_M_x[long_lag - 1] == 0) ? 1 : 0;
-      _M_p = 0;
-    }
+       for (size_t __i = 0; __i < long_lag; ++__i)
+         {
+           _UIntType __sum = 0u;
+           _UIntType __factor = 1u;
+           for (size_t __j = 0; __j < __k; ++__j)
+             {
+               __sum += __arr[__k * __i + __j] * __factor;
+               __factor *= __detail::_Shift<_UIntType, 32>::__value;
+             }
+           _M_x[__i] = __detail::__mod<_UIntType,
+             __detail::_Shift<_UIntType, __w>::__value>(__sum);
+         }
+       _M_carry = (_M_x[long_lag - 1] == 0) ? 1 : 0;
+       _M_p = 0;
+      }
 
   template<typename _UIntType, size_t __w, size_t __s, size_t __r>
     typename subtract_with_carry_engine<_UIntType, __w, __s, __r>::
@@ -474,6 +655,14 @@ namespace std
     }
 
 
+  template<typename _RandomNumberEngine, size_t __p, size_t __r>
+    constexpr size_t
+    discard_block_engine<_RandomNumberEngine, __p, __r>::block_size;
+
+  template<typename _RandomNumberEngine, size_t __p, size_t __r>
+    constexpr size_t
+    discard_block_engine<_RandomNumberEngine, __p, __r>::used_block;
+
   template<typename _RandomNumberEngine, size_t __p, size_t __r>
     typename discard_block_engine<_RandomNumberEngine,
                           __p, __r>::result_type
@@ -576,6 +765,10 @@ namespace std
     }
 
 
+  template<typename _RandomNumberEngine, size_t __k>
+    constexpr size_t
+    shuffle_order_engine<_RandomNumberEngine, __k>::table_size;
+
   template<typename _RandomNumberEngine, size_t __k>
     typename shuffle_order_engine<_RandomNumberEngine, __k>::result_type
     shuffle_order_engine<_RandomNumberEngine, __k>::
@@ -643,29 +836,61 @@ namespace std
       operator()(_UniformRandomNumberGenerator& __urng,
                 const param_type& __param)
       {
-       // XXX Must be fixed to work well for *arbitrary* __urng.max(),
-       // __urng.min(), __param.b(), __param.a().  Currently works fine only
-       // in the most common case __urng.max() - __urng.min() >=
-       // __param.b() - __param.a(), with __urng.max() > __urng.min() >= 0.
-       typedef typename __gnu_cxx::__add_unsigned<typename
-         _UniformRandomNumberGenerator::result_type>::__type __urntype;
-       typedef typename __gnu_cxx::__add_unsigned<result_type>::__type
-                                                             __utype;
-       typedef typename __gnu_cxx::__conditional_type<(sizeof(__urntype)
-                                                       > sizeof(__utype)),
-         __urntype, __utype>::__type                         __uctype;
+       typedef typename std::make_unsigned<typename
+         _UniformRandomNumberGenerator::result_type>::type __urngtype;
+       typedef typename std::make_unsigned<result_type>::type __utype;
+       typedef typename std::conditional<(sizeof(__urngtype)
+                                          > sizeof(__utype)),
+         __urngtype, __utype>::type __uctype;
+
+       const __uctype __urngmin = __urng.min();
+       const __uctype __urngmax = __urng.max();
+       const __uctype __urngrange = __urngmax - __urngmin;
+       const __uctype __urange
+         = __uctype(__param.b()) - __uctype(__param.a());
+
+       __uctype __ret;
+
+       if (__urngrange > __urange)
+         {
+           // downscaling
+           const __uctype __uerange = __urange + 1; // __urange can be zero
+           const __uctype __scaling = __urngrange / __uerange;
+           const __uctype __past = __uerange * __scaling;
+           do
+             __ret = __uctype(__urng()) - __urngmin;
+           while (__ret >= __past);
+           __ret /= __scaling;
+         }
+       else if (__urngrange < __urange)
+         {
+           // upscaling
+           /*
+             Note that every value in [0, urange]
+             can be written uniquely as
 
-       result_type __ret;
+             (urngrange + 1) * high + low
 
-       const __urntype __urnmin = __urng.min();
-       const __urntype __urnmax = __urng.max();
-       const __urntype __urnrange = __urnmax - __urnmin;
-       const __uctype __urange = __param.b() - __param.a();
-       const __uctype __udenom = (__urnrange <= __urange
-                                  ? 1 : __urnrange / (__urange + 1));
-       do
-         __ret = (__urntype(__urng()) -  __urnmin) / __udenom;
-       while (__ret > __param.b() - __param.a());
+             where
+
+             high in [0, urange / (urngrange + 1)]
+
+             and
+       
+             low in [0, urngrange].
+           */
+           __uctype __tmp; // wraparound control
+           do
+             {
+               const __uctype __uerngrange = __urngrange + 1;
+               __tmp = (__uerngrange * operator()
+                        (__urng, param_type(0, __urange / __uerngrange)));
+               __ret = __tmp + (__uctype(__urng()) - __urngmin);
+             }
+           while (__ret > __urange || __ret < __tmp);
+         }
+       else
+         __ret = __uctype(__urng()) - __urngmin;
 
        return __ret + __param.a();
       }
@@ -726,7 +951,7 @@ namespace std
       const _CharT __space = __os.widen(' ');
       __os.flags(__ios_base::scientific | __ios_base::left);
       __os.fill(__space);
-      __os.precision(std::numeric_limits<_RealType>::digits10 + 1);
+      __os.precision(std::numeric_limits<_RealType>::max_digits10);
 
       __os << __x.a() << __space << __x.b();
 
@@ -770,7 +995,7 @@ namespace std
       const std::streamsize __precision = __os.precision();
       __os.flags(__ios_base::scientific | __ios_base::left);
       __os.fill(__os.widen(' '));
-      __os.precision(std::numeric_limits<double>::digits10 + 1);
+      __os.precision(std::numeric_limits<double>::max_digits10);
 
       __os << __x.p();
 
@@ -800,7 +1025,7 @@ namespace std
 
        double __cand;
        do
-         __cand = std::ceil(std::log(__aurng()) / __param._M_log_p);
+         __cand = std::floor(std::log(__aurng()) / __param._M_log_1_p);
        while (__cand >= __thr);
 
        return result_type(__cand + __naf);
@@ -820,7 +1045,7 @@ namespace std
       const std::streamsize __precision = __os.precision();
       __os.flags(__ios_base::scientific | __ios_base::left);
       __os.fill(__os.widen(' '));
-      __os.precision(std::numeric_limits<double>::digits10 + 1);
+      __os.precision(std::numeric_limits<double>::max_digits10);
 
       __os << __x.p();
 
@@ -850,7 +1075,7 @@ namespace std
       return __is;
     }
 
-
+  // This is Leger's algorithm, also in Devroye, Ch. X, Example 1.5.
   template<typename _IntType>
     template<typename _UniformRandomNumberGenerator>
       typename negative_binomial_distribution<_IntType>::result_type
@@ -875,7 +1100,7 @@ namespace std
          param_type;
        
        const double __y =
-         _M_gd(__urng, param_type(__p.k(), __p.p() / (1.0 - __p.p())));
+         _M_gd(__urng, param_type(__p.k(), (1.0 - __p.p()) / __p.p()));
 
        std::poisson_distribution<result_type> __poisson(__y);
        return __poisson(__urng);
@@ -895,7 +1120,7 @@ namespace std
       const _CharT __space = __os.widen(' ');
       __os.flags(__ios_base::scientific | __ios_base::left);
       __os.fill(__os.widen(' '));
-      __os.precision(std::numeric_limits<double>::digits10 + 1);
+      __os.precision(std::numeric_limits<double>::max_digits10);
 
       __os << __x.k() << __space << __x.p()
           << __space << __x._M_gd;
@@ -965,7 +1190,7 @@ namespace std
    * is defined.
    *
    * Reference:
-   * Devroye, L. "Non-Uniform Random Variates Generation." Springer-Verlag,
+   * Devroye, L. Non-Uniform Random Variates Generation. Springer-Verlag,
    * New York, 1986, Ch. X, Sects. 3.3 & 3.4 (+ Errata!).
    */
   template<typename _IntType>
@@ -1085,7 +1310,7 @@ namespace std
       const _CharT __space = __os.widen(' ');
       __os.flags(__ios_base::scientific | __ios_base::left);
       __os.fill(__space);
-      __os.precision(std::numeric_limits<double>::digits10 + 1);
+      __os.precision(std::numeric_limits<double>::max_digits10);
 
       __os << __x.mean() << __space << __x._M_nd;
 
@@ -1197,7 +1422,7 @@ namespace std
    * is defined.
    *
    * Reference:
-   * Devroye, L. "Non-Uniform Random Variates Generation." Springer-Verlag,
+   * Devroye, L. Non-Uniform Random Variates Generation. Springer-Verlag,
    * New York, 1986, Ch. X, Sect. 4 (+ Errata!).
    */
   template<typename _IntType>
@@ -1209,7 +1434,7 @@ namespace std
       {
        result_type __ret;
        const _IntType __t = __param.t();
-       const _IntType __p = __param.p();
+       const double __p = __param.p();
        const double __p12 = __p <= 0.5 ? __p : 1.0 - __p;
        __detail::_Adaptor<_UniformRandomNumberGenerator, double>
          __aurng(__urng);
@@ -1333,7 +1558,7 @@ namespace std
       const _CharT __space = __os.widen(' ');
       __os.flags(__ios_base::scientific | __ios_base::left);
       __os.fill(__space);
-      __os.precision(std::numeric_limits<double>::digits10 + 1);
+      __os.precision(std::numeric_limits<double>::max_digits10);
 
       __os << __x.t() << __space << __x.p()
           << __space << __x._M_nd;
@@ -1380,7 +1605,7 @@ namespace std
       const std::streamsize __precision = __os.precision();
       __os.flags(__ios_base::scientific | __ios_base::left);
       __os.fill(__os.widen(' '));
-      __os.precision(std::numeric_limits<_RealType>::digits10 + 1);
+      __os.precision(std::numeric_limits<_RealType>::max_digits10);
 
       __os << __x.lambda();
 
@@ -1414,7 +1639,7 @@ namespace std
   /**
    * Polar method due to Marsaglia.
    *
-   * Devroye, L. "Non-Uniform Random Variates Generation." Springer-Verlag,
+   * Devroye, L. Non-Uniform Random Variates Generation. Springer-Verlag,
    * New York, 1986, Ch. V, Sect. 4.4.
    */
   template<typename _RealType>
@@ -1454,6 +1679,26 @@ namespace std
        return __ret;
       }
 
+  template<typename _RealType>
+    bool
+    operator==(const std::normal_distribution<_RealType>& __d1,
+              const std::normal_distribution<_RealType>& __d2)
+    {
+      if (__d1._M_param == __d2._M_param
+         && __d1._M_saved_available == __d2._M_saved_available)
+       {
+         if (__d1._M_saved_available
+             && __d1._M_saved == __d2._M_saved)
+           return true;
+         else if(!__d1._M_saved_available)
+           return true;
+         else
+           return false;
+       }
+      else
+       return false;
+    }
+
   template<typename _RealType, typename _CharT, typename _Traits>
     std::basic_ostream<_CharT, _Traits>&
     operator<<(std::basic_ostream<_CharT, _Traits>& __os,
@@ -1468,7 +1713,7 @@ namespace std
       const _CharT __space = __os.widen(' ');
       __os.flags(__ios_base::scientific | __ios_base::left);
       __os.fill(__space);
-      __os.precision(std::numeric_limits<_RealType>::digits10 + 1);
+      __os.precision(std::numeric_limits<_RealType>::max_digits10);
 
       __os << __x.mean() << __space << __x.stddev()
           << __space << __x._M_saved_available;
@@ -1519,7 +1764,7 @@ namespace std
       const _CharT __space = __os.widen(' ');
       __os.flags(__ios_base::scientific | __ios_base::left);
       __os.fill(__space);
-      __os.precision(std::numeric_limits<_RealType>::digits10 + 1);
+      __os.precision(std::numeric_limits<_RealType>::max_digits10);
 
       __os << __x.m() << __space << __x.s()
           << __space << __x._M_nd;
@@ -1565,7 +1810,7 @@ namespace std
       const _CharT __space = __os.widen(' ');
       __os.flags(__ios_base::scientific | __ios_base::left);
       __os.fill(__space);
-      __os.precision(std::numeric_limits<_RealType>::digits10 + 1);
+      __os.precision(std::numeric_limits<_RealType>::max_digits10);
 
       __os << __x.n() << __space << __x._M_gd;
 
@@ -1628,7 +1873,7 @@ namespace std
       const _CharT __space = __os.widen(' ');
       __os.flags(__ios_base::scientific | __ios_base::left);
       __os.fill(__space);
-      __os.precision(std::numeric_limits<_RealType>::digits10 + 1);
+      __os.precision(std::numeric_limits<_RealType>::max_digits10);
 
       __os << __x.a() << __space << __x.b();
 
@@ -1673,7 +1918,7 @@ namespace std
       const _CharT __space = __os.widen(' ');
       __os.flags(__ios_base::scientific | __ios_base::left);
       __os.fill(__space);
-      __os.precision(std::numeric_limits<_RealType>::digits10 + 1);
+      __os.precision(std::numeric_limits<_RealType>::max_digits10);
 
       __os << __x.m() << __space << __x.n()
           << __space << __x._M_gd_x << __space << __x._M_gd_y;
@@ -1719,7 +1964,7 @@ namespace std
       const _CharT __space = __os.widen(' ');
       __os.flags(__ios_base::scientific | __ios_base::left);
       __os.fill(__space);
-      __os.precision(std::numeric_limits<_RealType>::digits10 + 1);
+      __os.precision(std::numeric_limits<_RealType>::max_digits10);
 
       __os << __x.n() << __space << __x._M_nd << __space << __x._M_gd;
 
@@ -1822,7 +2067,7 @@ namespace std
       const _CharT __space = __os.widen(' ');
       __os.flags(__ios_base::scientific | __ios_base::left);
       __os.fill(__space);
-      __os.precision(std::numeric_limits<_RealType>::digits10 + 1);
+      __os.precision(std::numeric_limits<_RealType>::max_digits10);
 
       __os << __x.alpha() << __space << __x.beta()
           << __space << __x._M_nd;
@@ -1881,7 +2126,7 @@ namespace std
       const _CharT __space = __os.widen(' ');
       __os.flags(__ios_base::scientific | __ios_base::left);
       __os.fill(__space);
-      __os.precision(std::numeric_limits<_RealType>::digits10 + 1);
+      __os.precision(std::numeric_limits<_RealType>::max_digits10);
 
       __os << __x.a() << __space << __x.b();
 
@@ -1938,7 +2183,7 @@ namespace std
       const _CharT __space = __os.widen(' ');
       __os.flags(__ios_base::scientific | __ios_base::left);
       __os.fill(__space);
-      __os.precision(std::numeric_limits<_RealType>::digits10 + 1);
+      __os.precision(std::numeric_limits<_RealType>::max_digits10);
 
       __os << __x.a() << __space << __x.b();
 
@@ -1977,15 +2222,14 @@ namespace std
       if (_M_prob.size() < 2)
        {
          _M_prob.clear();
-         _M_prob.push_back(1.0);
          return;
        }
 
       const double __sum = std::accumulate(_M_prob.begin(),
                                           _M_prob.end(), 0.0);
       // Now normalize the probabilites.
-      std::transform(_M_prob.begin(), _M_prob.end(), _M_prob.begin(),
-                    std::bind2nd(std::divides<double>(), __sum));
+      __detail::__transform(_M_prob.begin(), _M_prob.end(), _M_prob.begin(),
+                         std::bind2nd(std::divides<double>(), __sum));
       // Accumulate partial sums.
       _M_cp.reserve(_M_prob.size());
       std::partial_sum(_M_prob.begin(), _M_prob.end(),
@@ -2017,6 +2261,9 @@ namespace std
       operator()(_UniformRandomNumberGenerator& __urng,
                 const param_type& __param)
       {
+       if (__param._M_cp.empty())
+         return result_type(0);
+
        __detail::_Adaptor<_UniformRandomNumberGenerator, double>
          __aurng(__urng);
 
@@ -2041,7 +2288,7 @@ namespace std
       const _CharT __space = __os.widen(' ');
       __os.flags(__ios_base::scientific | __ios_base::left);
       __os.fill(__space);
-      __os.precision(std::numeric_limits<double>::digits10 + 1);
+      __os.precision(std::numeric_limits<double>::max_digits10);
 
       std::vector<double> __prob = __x.probabilities();
       __os << __prob.size();
@@ -2090,24 +2337,21 @@ namespace std
     piecewise_constant_distribution<_RealType>::param_type::
     _M_initialize()
     {
-      if (_M_int.size() < 2)
+      if (_M_int.size() < 2
+         || (_M_int.size() == 2
+             && _M_int[0] == _RealType(0)
+             && _M_int[1] == _RealType(1)))
        {
          _M_int.clear();
-         _M_int.reserve(2);
-         _M_int.push_back(_RealType(0));
-         _M_int.push_back(_RealType(1));
-
          _M_den.clear();
-         _M_den.push_back(1.0);
-
          return;
        }
 
       const double __sum = std::accumulate(_M_den.begin(),
                                           _M_den.end(), 0.0);
 
-      std::transform(_M_den.begin(), _M_den.end(), _M_den.begin(),
-                    std::bind2nd(std::divides<double>(), __sum));
+      __detail::__transform(_M_den.begin(), _M_den.end(), _M_den.begin(),
+                           std::bind2nd(std::divides<double>(), __sum));
 
       _M_cp.reserve(_M_den.size());
       std::partial_sum(_M_den.begin(), _M_den.end(),
@@ -2193,6 +2437,9 @@ namespace std
          __aurng(__urng);
 
        const double __p = __aurng();
+       if (__param._M_cp.empty())
+         return __p;
+
        auto __pos = std::lower_bound(__param._M_cp.begin(),
                                      __param._M_cp.end(), __p);
        const size_t __i = __pos - __param._M_cp.begin();
@@ -2216,7 +2463,7 @@ namespace std
       const _CharT __space = __os.widen(' ');
       __os.flags(__ios_base::scientific | __ios_base::left);
       __os.fill(__space);
-      __os.precision(std::numeric_limits<_RealType>::digits10 + 1);
+      __os.precision(std::numeric_limits<_RealType>::max_digits10);
 
       std::vector<_RealType> __int = __x.intervals();
       __os << __int.size() - 1;
@@ -2279,18 +2526,14 @@ namespace std
     piecewise_linear_distribution<_RealType>::param_type::
     _M_initialize()
     {
-      if (_M_int.size() < 2)
+      if (_M_int.size() < 2
+         || (_M_int.size() == 2
+             && _M_int[0] == _RealType(0)
+             && _M_int[1] == _RealType(1)
+             && _M_den[0] == _M_den[1]))
        {
          _M_int.clear();
-         _M_int.reserve(2);
-         _M_int.push_back(_RealType(0));
-         _M_int.push_back(_RealType(1));
-
          _M_den.clear();
-         _M_den.reserve(2);
-         _M_den.push_back(1.0);
-         _M_den.push_back(1.0);
-
          return;
        }
 
@@ -2306,14 +2549,14 @@ namespace std
        }
 
       //  Now normalize the densities...
-      std::transform(_M_den.begin(), _M_den.end(), _M_den.begin(),
-                    std::bind2nd(std::divides<double>(), __sum));
+      __detail::__transform(_M_den.begin(), _M_den.end(), _M_den.begin(),
+                         std::bind2nd(std::divides<double>(), __sum));
       //  ... and partial sums... 
-      std::transform(_M_cp.begin(), _M_cp.end(), _M_cp.begin(),
-                    std::bind2nd(std::divides<double>(), __sum));
+      __detail::__transform(_M_cp.begin(), _M_cp.end(), _M_cp.begin(),
+                           std::bind2nd(std::divides<double>(), __sum));
       //  ... and slopes.
-      std::transform(_M_m.begin(), _M_m.end(), _M_m.begin(),
-                    std::bind2nd(std::divides<double>(), __sum));
+      __detail::__transform(_M_m.begin(), _M_m.end(), _M_m.begin(),
+                           std::bind2nd(std::divides<double>(), __sum));
       //  Make sure the last cumulative probablility is one.
       _M_cp[_M_cp.size() - 1] = 1.0;
      }
@@ -2383,6 +2626,9 @@ namespace std
          __aurng(__urng);
 
        const double __p = __aurng();
+       if (__param._M_cp.empty())
+         return __p;
+
        auto __pos = std::lower_bound(__param._M_cp.begin(),
                                      __param._M_cp.end(), __p);
        const size_t __i = __pos - __param._M_cp.begin();
@@ -2419,7 +2665,7 @@ namespace std
       const _CharT __space = __os.widen(' ');
       __os.flags(__ios_base::scientific | __ios_base::left);
       __os.fill(__space);
-      __os.precision(std::numeric_limits<_RealType>::digits10 + 1);
+      __os.precision(std::numeric_limits<_RealType>::max_digits10);
 
       std::vector<_RealType> __int = __x.intervals();
       __os << __int.size() - 1;
@@ -2577,4 +2823,8 @@ namespace std
        }
       return __sum / __tmp;
     }
-}
+
+_GLIBCXX_END_NAMESPACE_VERSION
+} // namespace
+
+#endif