random.h (negative_binomial_distribution<>:: negative_binomial_distribution(_IntType...
[gcc.git] / libstdc++-v3 / include / bits / random.tcc
index eb2ce4a4500a95017ab650b00d5f3066133d8652..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.
     //
@@ -45,7 +49,7 @@ namespace std
     //
     // Preconditions:  a > 0, m > 0.
     //
-    template<typename _Tp, _Tp __a, _Tp __c, _Tp __m, bool>
+    template<typename _Tp, _Tp __m, _Tp __a, _Tp __c, bool>
       struct _Mod
       {
        static _Tp
@@ -80,15 +84,46 @@ namespace std
 
     // Special case for m == 0 -- use unsigned integer overflow as modulo
     // operator.
-    template<typename _Tp, _Tp __a, _Tp __c, _Tp __m>
-      struct _Mod<_Tp, __a, __c, __m, true>
+    template<typename _Tp, _Tp __m, _Tp __a, _Tp __c>
+      struct _Mod<_Tp, __m, __a, __c, true>
       {
        static _Tp
        __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.
@@ -98,35 +133,36 @@ namespace std
     linear_congruential_engine<_UIntType, __a, __c, __m>::
     seed(result_type __s)
     {
-      if ((__detail::__mod<_UIntType, 1U, 0U, __m>(__c) == 0U)
-         && (__detail::__mod<_UIntType, 1U, 0U, __m>(__s) == 0U))
-       _M_x = __detail::__mod<_UIntType, 1U, 0U, __m>(1U);
+      if ((__detail::__mod<_UIntType, __m>(__c) == 0)
+         && (__detail::__mod<_UIntType, __m>(__s) == 0))
+       _M_x = 1;
       else
-       _M_x = __detail::__mod<_UIntType, 1U, 0U, __m>(__s);
+       _M_x = __detail::__mod<_UIntType, __m>(__s);
     }
 
   /**
    * 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,
@@ -179,7 +342,7 @@ namespace std
                            __s, __b, __t, __c, __l, __f>::
     seed(result_type __sd)
     {
-      _M_x[0] = __detail::__mod<_UIntType, 1, 0,
+      _M_x[0] = __detail::__mod<_UIntType,
        __detail::_Shift<_UIntType, __w>::__value>(__sd);
 
       for (size_t __i = 1; __i < state_size; ++__i)
@@ -187,8 +350,8 @@ namespace std
          _UIntType __x = _M_x[__i - 1];
          __x ^= __x >> (__w - 2);
          __x *= __f;
-         __x += __i;
-         _M_x[__i] = __detail::__mod<_UIntType, 1, 0,
+         __x += __detail::__mod<_UIntType, __n>(__i);
+         _M_x[__i] = __detail::__mod<_UIntType,
            __detail::_Shift<_UIntType, __w>::__value>(__x);
        }
       _M_p = state_size;
@@ -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;
-      _UIntType __arr[__k * __n];
-      __q.generate(__arr + 0, __arr + __k * __n);
-
-      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[__i * __k + __j] * __factor;
-             __factor *= __detail::_Shift<_UIntType, 32>::__value;
-            }
-          _M_x[__i] = __detail::__mod<_UIntType, 1U, 0U,
-                     __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;
-            }
-        }
+      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);
+
+       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 - 1U>::__value;
-    }
+          _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,31 +504,44 @@ 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>::
     seed(result_type __value)
     {
-      if (__value == 0)
-       __value = default_seed;
-
-      std::linear_congruential_engine<result_type, 40014U, 0U, 2147483563U>
-       __lcg(__value);
+      std::linear_congruential_engine<result_type, 40014u, 0u, 2147483563u>
+       __lcg(__value == 0u ? default_seed : __value);
 
-      //  I hope this is right.  The "10000" tests work for the ranluxen.
-      const size_t __n = (word_size + 31) / 32;
+      const size_t __n = (__w + 31) / 32;
 
       for (size_t __i = 0; __i < long_lag; ++__i)
        {
-         _UIntType __sum = 0U;
-         _UIntType __factor = 1U;
+         _UIntType __sum = 0u;
+         _UIntType __factor = 1u;
          for (size_t __j = 0; __j < __n; ++__j)
            {
-             __sum += __detail::__mod<__detail::_UInt32Type, 1, 0, 0>
+             __sum += __detail::__mod<uint_least32_t,
+                      __detail::_Shift<uint_least32_t, 32>::__value>
                         (__lcg()) * __factor;
              __factor *= __detail::_Shift<_UIntType, 32>::__value;
            }
-         _M_x[__i] = __detail::__mod<_UIntType, 1, 0,
+         _M_x[__i] = __detail::__mod<_UIntType,
            __detail::_Shift<_UIntType, __w>::__value>(__sum);
        }
       _M_carry = (_M_x[long_lag - 1] == 0) ? 1 : 0;
@@ -372,30 +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 __n = (word_size + 31) / 32;
-      _UIntType __arr[long_lag + __n];
-      __q.generate(__arr + 0, __arr + long_lag + __n);
+    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 < __n; ++__j)
-            {
-             __sum += __detail::__mod<__detail::_UInt32Type, 1, 0, 0>
-                        (__arr[__i * __n + __j]) * __factor;
-             __factor *= __detail::_Shift<_UIntType, 32>::__value;
-            }
-          _M_x[__i] = __detail::__mod<_UIntType, 1, 0,
-           __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>::
@@ -478,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
@@ -580,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>::
@@ -644,34 +833,66 @@ namespace std
     template<typename _UniformRandomNumberGenerator>
       typename uniform_int_distribution<_IntType>::result_type
       uniform_int_distribution<_IntType>::
-      _M_call(_UniformRandomNumberGenerator& __urng,
-             result_type __min, result_type __max, true_type)
+      operator()(_UniformRandomNumberGenerator& __urng,
+                const param_type& __param)
       {
-       // XXX Must be fixed to work well for *arbitrary* __urng.max(),
-       // __urng.min(), __max, __min.  Currently works fine only in the
-       // most common case __urng.max() - __urng.min() >= __max - __min,
-       // 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 = __max - __min;
-       const __uctype __udenom = (__urnrange <= __urange
-                                  ? 1 : __urnrange / (__urange + 1));
-       do
-         __ret = (__urntype(__urng()) -  __urnmin) / __udenom;
-       while (__ret > __max - __min);
+             where
 
-       return __ret + __min;
+             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();
       }
 
   template<typename _IntType, typename _CharT, typename _Traits>
@@ -730,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();
 
@@ -774,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();
 
@@ -799,12 +1020,12 @@ namespace std
        // The largest _RealType convertible to _IntType.
        const double __thr =
          std::numeric_limits<_IntType>::max() + __naf;
-       __detail::_Adaptor<_UniformRandomNumberGenerator, result_type>
+       __detail::_Adaptor<_UniformRandomNumberGenerator, double>
          __aurng(__urng);
 
        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);
@@ -824,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();
 
@@ -854,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
@@ -879,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);
@@ -899,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;
@@ -969,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>
@@ -1089,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;
 
@@ -1201,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>
@@ -1213,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);
@@ -1337,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;
@@ -1384,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();
 
@@ -1418,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>
@@ -1458,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,
@@ -1472,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;
@@ -1523,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;
@@ -1569,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;
 
@@ -1632,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();
 
@@ -1677,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;
@@ -1723,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;
 
@@ -1826,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;
@@ -1885,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();
 
@@ -1942,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();
 
@@ -1981,34 +2222,34 @@ namespace std
       if (_M_prob.size() < 2)
        {
          _M_prob.clear();
-         _M_prob.push_back(1.0);
          return;
        }
 
-      double __sum = std::accumulate(_M_prob.begin(), _M_prob.end(), 0.0);
-      //  Now normalize the densities.
-      std::transform(_M_prob.begin(), _M_prob.end(), _M_prob.begin(),
-                    std::bind2nd(std::divides<double>(), __sum));
-      //  Accumulate partial sums.
+      const double __sum = std::accumulate(_M_prob.begin(),
+                                          _M_prob.end(), 0.0);
+      // Now normalize the probabilites.
+      __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(),
                       std::back_inserter(_M_cp));
-      //  Make sure the last cumulative probablility is one.
+      // Make sure the last cumulative probability is one.
       _M_cp[_M_cp.size() - 1] = 1.0;
     }
 
   template<typename _IntType>
     template<typename _Func>
       discrete_distribution<_IntType>::param_type::
-      param_type(size_t __nw, double __xmin, double __xmax,
-                _Func __fw)
+      param_type(size_t __nw, double __xmin, double __xmax, _Func __fw)
       : _M_prob(), _M_cp()
       {
-       for (size_t __i = 0; __i < __nw; ++__i)
-         {
-           const double __x = ((__nw - __i - 0.5) * __xmin
-                                    + (__i + 0.5) * __xmax) / __nw;
-           _M_prob.push_back(__fw(__x));
-         }
+       const size_t __n = __nw == 0 ? 1 : __nw;
+       const double __delta = (__xmax - __xmin) / __n;
+
+       _M_prob.reserve(__n);
+       for (size_t __k = 0; __k < __nw; ++__k)
+         _M_prob.push_back(__fw(__xmin + __k * __delta + 0.5 * __delta));
 
        _M_initialize();
       }
@@ -2020,17 +2261,17 @@ namespace std
       operator()(_UniformRandomNumberGenerator& __urng,
                 const param_type& __param)
       {
-       __detail::_Adaptor<_UniformRandomNumberGenerator, result_type>
+       if (__param._M_cp.empty())
+         return result_type(0);
+
+       __detail::_Adaptor<_UniformRandomNumberGenerator, double>
          __aurng(__urng);
 
        const double __p = __aurng();
        auto __pos = std::lower_bound(__param._M_cp.begin(),
                                      __param._M_cp.end(), __p);
-       if (__pos == __param._M_cp.end())
-         return 0;
-       const size_t __i = __pos - __param._M_cp.begin();
 
-       return __i;
+       return __pos - __param._M_cp.begin();
       }
 
   template<typename _IntType, typename _CharT, typename _Traits>
@@ -2047,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();
@@ -2075,6 +2316,7 @@ namespace std
       __is >> __n;
 
       std::vector<double> __prob_vec;
+      __prob_vec.reserve(__n);
       for (; __n != 0; --__n)
        {
          double __prob;
@@ -2095,33 +2337,31 @@ 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.push_back(_RealType(0));
-         _M_int.push_back(_RealType(1));
-
          _M_den.clear();
-         _M_den.push_back(1.0);
-
          return;
        }
 
-      double __sum = 0.0;
-      for (size_t __i = 0; __i < _M_den.size(); ++__i)
-       {
-         __sum += _M_den[__i] * (_M_int[__i + 1] - _M_int[__i]);
-         _M_cp.push_back(__sum);
-       }
+      const double __sum = std::accumulate(_M_den.begin(),
+                                          _M_den.end(), 0.0);
 
-      //  Now normalize the densities...
-      std::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));
-      //  Make sure the last cumulative probablility is one.
+      __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(),
+                      std::back_inserter(_M_cp));
+
+      // Make sure the last cumulative probability is one.
       _M_cp[_M_cp.size() - 1] = 1.0;
+
+      for (size_t __k = 0; __k < _M_den.size(); ++__k)
+       _M_den[__k] /= _M_int[__k + 1] - _M_int[__k];
     }
 
   template<typename _RealType>
@@ -2132,17 +2372,19 @@ namespace std
                 _InputIteratorW __wbegin)
       : _M_int(), _M_den(), _M_cp()
       {
-       do
+       if (__bbegin != __bend)
          {
-           _M_int.push_back(*__bbegin);
-           ++__bbegin;
-           if (__bbegin != __bend)
+           for (;;)
              {
+               _M_int.push_back(*__bbegin);
+               ++__bbegin;
+               if (__bbegin == __bend)
+                 break;
+
                _M_den.push_back(*__wbegin);
                ++__wbegin;
              }
          }
-       while (__bbegin != __bend);
 
        _M_initialize();
       }
@@ -2150,17 +2392,16 @@ namespace std
   template<typename _RealType>
     template<typename _Func>
       piecewise_constant_distribution<_RealType>::param_type::
-      param_type(initializer_list<_RealType> __bil, _Func __fw)
+      param_type(initializer_list<_RealType> __bl, _Func __fw)
       : _M_int(), _M_den(), _M_cp()
       {
-       for (auto __biter = __bil.begin(); __biter != __bil.end(); ++__biter)
+       _M_int.reserve(__bl.size());
+       for (auto __biter = __bl.begin(); __biter != __bl.end(); ++__biter)
          _M_int.push_back(*__biter);
 
-       for (size_t __i = 0; __i < _M_int.size() - 1; ++__i)
-         {
-           _RealType __x = 0.5 * (_M_int[__i] + _M_int[__i + 1]);
-           _M_den.push_back(__fw(__x));
-         }
+       _M_den.reserve(_M_int.size() - 1);
+       for (size_t __k = 0; __k < _M_int.size() - 1; ++__k)
+         _M_den.push_back(__fw(0.5 * (_M_int[__k + 1] + _M_int[__k])));
 
        _M_initialize();
       }
@@ -2171,18 +2412,16 @@ namespace std
       param_type(size_t __nw, _RealType __xmin, _RealType __xmax, _Func __fw)
       : _M_int(), _M_den(), _M_cp()
       {
-       for (size_t __i = 0; __i <= __nw; ++__i)
-         {
-           const _RealType __x = ((__nw - __i) * __xmin
-                                         + __i * __xmax) / __nw;
-           _M_int.push_back(__x);
-         }
-       for (size_t __i = 0; __i < __nw; ++__i)
-         {
-           const _RealType __x = ((__nw - __i - 0.5) * __xmin
-                                       + (__i + 0.5) * __xmax) / __nw;
-           _M_den.push_back(__fw(__x));
-         }
+       const size_t __n = __nw == 0 ? 1 : __nw;
+       const _RealType __delta = (__xmax - __xmin) / __n;
+
+       _M_int.reserve(__n + 1);
+       for (size_t __k = 0; __k <= __nw; ++__k)
+         _M_int.push_back(__xmin + __k * __delta);
+
+       _M_den.reserve(__n);
+       for (size_t __k = 0; __k < __nw; ++__k)
+         _M_den.push_back(__fw(_M_int[__k] + 0.5 * __delta));
 
        _M_initialize();
       }
@@ -2194,16 +2433,20 @@ namespace std
       operator()(_UniformRandomNumberGenerator& __urng,
                 const param_type& __param)
       {
-       __detail::_Adaptor<_UniformRandomNumberGenerator, result_type>
+       __detail::_Adaptor<_UniformRandomNumberGenerator, double>
          __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();
 
-       return __param._M_int[__i]
-            + (__p - __param._M_cp[__i]) / __param._M_den[__i];
+       const double __pref = __i > 0 ? __param._M_cp[__i - 1] : 0.0;
+
+       return __param._M_int[__i] + (__p - __pref) / __param._M_den[__i];
       }
 
   template<typename _RealType, typename _CharT, typename _Traits>
@@ -2220,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;
@@ -2253,6 +2496,7 @@ namespace std
       __is >> __n;
 
       std::vector<_RealType> __int_vec;
+      __int_vec.reserve(__n + 1);
       for (size_t __i = 0; __i <= __n; ++__i)
        {
          _RealType __int;
@@ -2261,6 +2505,7 @@ namespace std
        }
 
       std::vector<double> __den_vec;
+      __den_vec.reserve(__n);
       for (size_t __i = 0; __i < __n; ++__i)
        {
          double __den;
@@ -2281,46 +2526,40 @@ 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.push_back(_RealType(0));
-         _M_int.push_back(_RealType(1));
-
          _M_den.clear();
-         _M_den.push_back(1.0);
-         _M_den.push_back(1.0);
-
          return;
        }
 
       double __sum = 0.0;
-      for (size_t __i = 0; __i < _M_int.size() - 1; ++__i)
+      _M_cp.reserve(_M_int.size() - 1);
+      _M_m.reserve(_M_int.size() - 1);
+      for (size_t __k = 0; __k < _M_int.size() - 1; ++__k)
        {
-         const _RealType __delta = _M_int[__i + 1] - _M_int[__i];
-         __sum += 0.5 * (_M_den[__i + 1] + _M_den[__i]) * __delta;
+         const _RealType __delta = _M_int[__k + 1] - _M_int[__k];
+         __sum += 0.5 * (_M_den[__k + 1] + _M_den[__k]) * __delta;
          _M_cp.push_back(__sum);
-         _M_m.push_back((_M_den[__i + 1] - _M_den[__i]) / __delta);
+         _M_m.push_back((_M_den[__k + 1] - _M_den[__k]) / __delta);
        }
 
       //  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;
-    }
-
-  template<typename _RealType>
-    piecewise_linear_distribution<_RealType>::param_type::
-    param_type()
-    : _M_int(), _M_den(), _M_cp(), _M_m()
-    { _M_initialize(); }
+     }
 
   template<typename _RealType>
     template<typename _InputIteratorB, typename _InputIteratorW>
@@ -2342,10 +2581,12 @@ namespace std
   template<typename _RealType>
     template<typename _Func>
       piecewise_linear_distribution<_RealType>::param_type::
-      param_type(initializer_list<_RealType> __bil, _Func __fw)
+      param_type(initializer_list<_RealType> __bl, _Func __fw)
       : _M_int(), _M_den(), _M_cp(), _M_m()
       {
-       for (auto __biter = __bil.begin(); __biter != __bil.end(); ++__biter)
+       _M_int.reserve(__bl.size());
+       _M_den.reserve(__bl.size());
+       for (auto __biter = __bl.begin(); __biter != __bl.end(); ++__biter)
          {
            _M_int.push_back(*__biter);
            _M_den.push_back(__fw(*__biter));
@@ -2357,16 +2598,18 @@ namespace std
   template<typename _RealType>
     template<typename _Func>
       piecewise_linear_distribution<_RealType>::param_type::
-      param_type(size_t __nw, _RealType __xmin, _RealType __xmax,
-                _Func __fw)
+      param_type(size_t __nw, _RealType __xmin, _RealType __xmax, _Func __fw)
       : _M_int(), _M_den(), _M_cp(), _M_m()
       {
-       for (size_t __i = 0; __i <= __nw; ++__i)
+       const size_t __n = __nw == 0 ? 1 : __nw;
+       const _RealType __delta = (__xmax - __xmin) / __n;
+
+       _M_int.reserve(__n + 1);
+       _M_den.reserve(__n + 1);
+       for (size_t __k = 0; __k <= __nw; ++__k)
          {
-           const _RealType __x = ((__nw - __i) * __xmin
-                                         + __i * __xmax) / __nw;
-           _M_int.push_back(__x);
-           _M_den.push_back(__fw(__x));
+           _M_int.push_back(__xmin + __k * __delta);
+           _M_den.push_back(__fw(_M_int[__k] + __delta));
          }
 
        _M_initialize();
@@ -2379,31 +2622,33 @@ namespace std
       operator()(_UniformRandomNumberGenerator& __urng,
                 const param_type& __param)
       {
-       result_type __x;
-       __detail::_Adaptor<_UniformRandomNumberGenerator, result_type>
+       __detail::_Adaptor<_UniformRandomNumberGenerator, double>
          __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();
+
+       const double __pref = __i > 0 ? __param._M_cp[__i - 1] : 0.0;
+
        const double __a = 0.5 * __param._M_m[__i];
        const double __b = __param._M_den[__i];
-       const double __c = __param._M_cp[__i];
-       const double __q = -0.5 * (__b
-#if _GLIBCXX_USE_C99_MATH_TR1
-                        + std::copysign(std::sqrt(__b * __b
-                                                - 4.0 * __a * __c), __b));
-#else
-                        + (__b < 0.0 ? -1.0 : 1.0)
-                        * std::sqrt(__b * __b - 4.0 * __a * __c));
-#endif
-       const double __x0 = __param._M_int[__i];
-       const double __x1 = __q / __a;
-       const double __x2 = __c / __q;
-       __x = std::max(__x0 + __x1, __x0 + __x2);
+       const double __cm = __p - __pref;
+
+       _RealType __x = __param._M_int[__i];
+       if (__a == 0)
+         __x += __cm / __b;
+       else
+         {
+           const double __d = __b * __b + 4.0 * __a * __cm;
+           __x += 0.5 * (std::sqrt(__d) - __b) / __a;
+          }
 
-       return __x;
+        return __x;
       }
 
   template<typename _RealType, typename _CharT, typename _Traits>
@@ -2420,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;
@@ -2453,6 +2698,7 @@ namespace std
       __is >> __n;
 
       std::vector<_RealType> __int_vec;
+      __int_vec.reserve(__n + 1);
       for (size_t __i = 0; __i <= __n; ++__i)
        {
          _RealType __int;
@@ -2461,6 +2707,7 @@ namespace std
        }
 
       std::vector<double> __den_vec;
+      __den_vec.reserve(__n + 1);
       for (size_t __i = 0; __i <= __n; ++__i)
        {
          double __den;
@@ -2480,7 +2727,7 @@ namespace std
     seed_seq::seed_seq(std::initializer_list<_IntType> __il)
     {
       for (auto __iter = __il.begin(); __iter != __il.end(); ++__iter)
-       _M_v.push_back(__detail::__mod<result_type, 1, 0,
+       _M_v.push_back(__detail::__mod<result_type,
                       __detail::_Shift<result_type, 32>::__value>(*__iter));
     }
 
@@ -2488,7 +2735,7 @@ namespace std
     seed_seq::seed_seq(_InputIterator __begin, _InputIterator __end)
     {
       for (_InputIterator __iter = __begin; __iter != __end; ++__iter)
-       _M_v.push_back(__detail::__mod<result_type, 1, 0,
+       _M_v.push_back(__detail::__mod<result_type,
                       __detail::_Shift<result_type, 32>::__value>(*__iter));
     }
 
@@ -2503,7 +2750,7 @@ namespace std
       if (__begin == __end)
        return;
 
-      std::fill(__begin, __end, _Type(0x8b8b8b8bU));
+      std::fill(__begin, __end, _Type(0x8b8b8b8bu));
 
       const size_t __n = __end - __begin;
       const size_t __s = _M_v.size();
@@ -2522,8 +2769,8 @@ namespace std
                         ^ __begin[(__k + __p) % __n]
                         ^ __begin[(__k - 1) % __n]);
          _Type __r1 = __arg ^ (__arg << 27);
-         __r1 = __detail::__mod<_Type, 1664525U, 0U,
-                  __detail::_Shift<_Type, 32>::__value>(__r1);
+         __r1 = __detail::__mod<_Type, __detail::_Shift<_Type, 32>::__value,
+                                1664525u, 0u>(__r1);
          _Type __r2 = __r1;
          if (__k == 0)
            __r2 += __s;
@@ -2531,8 +2778,8 @@ namespace std
            __r2 += __k % __n + _M_v[__k - 1];
          else
            __r2 += __k % __n;
-         __r2 = __detail::__mod<_Type, 1U, 0U,
-                  __detail::_Shift<_Type, 32>::__value>(__r2);
+         __r2 = __detail::__mod<_Type,
+                  __detail::_Shift<_Type, 32>::__value>(__r2);
          __begin[(__k + __p) % __n] += __r1;
          __begin[(__k + __q) % __n] += __r2;
          __begin[__k % __n] = __r2;
@@ -2544,11 +2791,11 @@ namespace std
                         + __begin[(__k + __p) % __n]
                         + __begin[(__k - 1) % __n]);
          _Type __r3 = __arg ^ (__arg << 27);
-         __r3 = __detail::__mod<_Type, 1566083941U, 0U,
-                  __detail::_Shift<_Type, 32>::__value>(__r3);
+         __r3 = __detail::__mod<_Type, __detail::_Shift<_Type, 32>::__value,
+                                1566083941u, 0u>(__r3);
          _Type __r4 = __r3 - __k % __n;
-         __r4 = __detail::__mod<_Type, 1U, 0U,
-                  __detail::_Shift<_Type, 32>::__value>(__r4);
+         __r4 = __detail::__mod<_Type,
+                  __detail::_Shift<_Type, 32>::__value>(__r4);
          __begin[(__k + __p) % __n] ^= __r4;
          __begin[(__k + __q) % __n] ^= __r3;
          __begin[__k % __n] = __r4;
@@ -2576,4 +2823,8 @@ namespace std
        }
       return __sum / __tmp;
     }
-}
+
+_GLIBCXX_END_NAMESPACE_VERSION
+} // namespace
+
+#endif