diff --git a/include/random b/include/random
index dac6dd1..6d17c2c 100644
--- a/include/random
+++ b/include/random
@@ -3434,9 +3434,7 @@
     void param(const param_type& __p) {__p_ = __p;}
 
     result_type min() const {return 0;}
-    result_type max() const
-        {return -std::log(1-std::nextafter(result_type(1), result_type(-1))) /
-                 __p_.lambda();}
+    result_type max() const {return numeric_limits<result_type>::infinity();}
 
     friend bool operator==(const exponential_distribution& __x,
                            const exponential_distribution& __y)
diff --git a/test/numerics/rand/rand.dis/rand.dist.bern/rand.dist.bern.bernoulli/eval.pass.cpp b/test/numerics/rand/rand.dis/rand.dist.bern/rand.dist.bern.bernoulli/eval.pass.cpp
index c8e9db1..6764efa 100644
--- a/test/numerics/rand/rand.dis/rand.dist.bern/rand.dist.bern.bernoulli/eval.pass.cpp
+++ b/test/numerics/rand/rand.dis/rand.dist.bern/rand.dist.bern.bernoulli/eval.pass.cpp
@@ -40,13 +40,29 @@
         double mean = std::accumulate(u.begin(), u.end(),
                                               double(0)) / u.size();
         double var = 0;
+        double skew = 0;
+        double kurtosis = 0;
         for (int i = 0; i < u.size(); ++i)
-            var += sqr(u[i] - mean);
+        {
+            double d = (u[i] - mean);
+            double d2 = sqr(d);
+            var += d2;
+            skew += d * d2;
+            kurtosis += d2 * d2;
+        }
         var /= u.size();
+        double dev = std::sqrt(var);
+        skew /= u.size() * dev * var;
+        kurtosis /= u.size() * var * var;
+        kurtosis -= 3;
         double x_mean = d.p();
         double x_var = d.p()*(1-d.p());
+        double x_skew = (1 - 2 * d.p())/std::sqrt(x_var);
+        double x_kurtosis = (6 * sqr(d.p()) - 6 * d.p() + 1)/x_var;
         assert(std::abs(mean - x_mean) / x_mean < 0.01);
         assert(std::abs(var - x_var) / x_var < 0.01);
+        assert(std::abs(skew - x_skew) / x_skew < 0.01);
+        assert(std::abs(kurtosis - x_kurtosis) / x_kurtosis < 0.01);
     }
     {
         typedef std::bernoulli_distribution D;
@@ -60,12 +76,28 @@
         double mean = std::accumulate(u.begin(), u.end(),
                                               double(0)) / u.size();
         double var = 0;
+        double skew = 0;
+        double kurtosis = 0;
         for (int i = 0; i < u.size(); ++i)
-            var += sqr(u[i] - mean);
+        {
+            double d = (u[i] - mean);
+            double d2 = sqr(d);
+            var += d2;
+            skew += d * d2;
+            kurtosis += d2 * d2;
+        }
         var /= u.size();
+        double dev = std::sqrt(var);
+        skew /= u.size() * dev * var;
+        kurtosis /= u.size() * var * var;
+        kurtosis -= 3;
         double x_mean = d.p();
         double x_var = d.p()*(1-d.p());
+        double x_skew = (1 - 2 * d.p())/std::sqrt(x_var);
+        double x_kurtosis = (6 * sqr(d.p()) - 6 * d.p() + 1)/x_var;
         assert(std::abs(mean - x_mean) / x_mean < 0.01);
         assert(std::abs(var - x_var) / x_var < 0.01);
+        assert(std::abs(skew - x_skew) / x_skew < 0.01);
+        assert(std::abs(kurtosis - x_kurtosis) / x_kurtosis < 0.01);
     }
 }
diff --git a/test/numerics/rand/rand.dis/rand.dist.bern/rand.dist.bern.bernoulli/eval_param.pass.cpp b/test/numerics/rand/rand.dis/rand.dist.bern/rand.dist.bern.bernoulli/eval_param.pass.cpp
index 03eacb6..bb4be85 100644
--- a/test/numerics/rand/rand.dis/rand.dist.bern/rand.dist.bern.bernoulli/eval_param.pass.cpp
+++ b/test/numerics/rand/rand.dis/rand.dist.bern/rand.dist.bern.bernoulli/eval_param.pass.cpp
@@ -42,13 +42,29 @@
         double mean = std::accumulate(u.begin(), u.end(),
                                               double(0)) / u.size();
         double var = 0;
+        double skew = 0;
+        double kurtosis = 0;
         for (int i = 0; i < u.size(); ++i)
-            var += sqr(u[i] - mean);
+        {
+            double d = (u[i] - mean);
+            double d2 = sqr(d);
+            var += d2;
+            skew += d * d2;
+            kurtosis += d2 * d2;
+        }
         var /= u.size();
+        double dev = std::sqrt(var);
+        skew /= u.size() * dev * var;
+        kurtosis /= u.size() * var * var;
+        kurtosis -= 3;
         double x_mean = p.p();
         double x_var = p.p()*(1-p.p());
+        double x_skew = (1 - 2 * p.p())/std::sqrt(x_var);
+        double x_kurtosis = (6 * sqr(p.p()) - 6 * p.p() + 1)/x_var;
         assert(std::abs(mean - x_mean) / x_mean < 0.01);
         assert(std::abs(var - x_var) / x_var < 0.01);
+        assert(std::abs(skew - x_skew) / x_skew < 0.01);
+        assert(std::abs(kurtosis - x_kurtosis) / x_kurtosis < 0.01);
     }
     {
         typedef std::bernoulli_distribution D;
@@ -64,12 +80,28 @@
         double mean = std::accumulate(u.begin(), u.end(),
                                               double(0)) / u.size();
         double var = 0;
+        double skew = 0;
+        double kurtosis = 0;
         for (int i = 0; i < u.size(); ++i)
-            var += sqr(u[i] - mean);
+        {
+            double d = (u[i] - mean);
+            double d2 = sqr(d);
+            var += d2;
+            skew += d * d2;
+            kurtosis += d2 * d2;
+        }
         var /= u.size();
+        double dev = std::sqrt(var);
+        skew /= u.size() * dev * var;
+        kurtosis /= u.size() * var * var;
+        kurtosis -= 3;
         double x_mean = p.p();
         double x_var = p.p()*(1-p.p());
+        double x_skew = (1 - 2 * p.p())/std::sqrt(x_var);
+        double x_kurtosis = (6 * sqr(p.p()) - 6 * p.p() + 1)/x_var;
         assert(std::abs(mean - x_mean) / x_mean < 0.01);
         assert(std::abs(var - x_var) / x_var < 0.01);
+        assert(std::abs(skew - x_skew) / x_skew < 0.01);
+        assert(std::abs(kurtosis - x_kurtosis) / x_kurtosis < 0.01);
     }
 }
diff --git a/test/numerics/rand/rand.dis/rand.dist.bern/rand.dist.bern.bin/eval.pass.cpp b/test/numerics/rand/rand.dis/rand.dist.bern/rand.dist.bern.bin/eval.pass.cpp
index 95e267c..7ff25d5 100644
--- a/test/numerics/rand/rand.dis/rand.dist.bern/rand.dist.bern.bin/eval.pass.cpp
+++ b/test/numerics/rand/rand.dis/rand.dist.bern/rand.dist.bern.bin/eval.pass.cpp
@@ -14,6 +14,8 @@
 
 // template<class _URNG> result_type operator()(_URNG& g);
 
+#include <iostream>
+
 #include <random>
 #include <numeric>
 #include <vector>
@@ -37,195 +39,385 @@
         const int N = 100000;
         std::vector<D::result_type> u;
         for (int i = 0; i < N; ++i)
-            u.push_back(d(g));
+        {
+            D::result_type v = d(g);
+            assert(d.min() <= v && v <= d.max());
+            u.push_back(v);
+        }
         double mean = std::accumulate(u.begin(), u.end(),
                                               double(0)) / u.size();
         double var = 0;
+        double skew = 0;
+        double kurtosis = 0;
         for (int i = 0; i < u.size(); ++i)
-            var += sqr(u[i] - mean);
+        {
+            double d = (u[i] - mean);
+            double d2 = sqr(d);
+            var += d2;
+            skew += d * d2;
+            kurtosis += d2 * d2;
+        }
         var /= u.size();
+        double dev = std::sqrt(var);
+        skew /= u.size() * dev * var;
+        kurtosis /= u.size() * var * var;
+        kurtosis -= 3;
         double x_mean = d.t() * d.p();
         double x_var = x_mean*(1-d.p());
+        double x_skew = (1-2*d.p()) / std::sqrt(x_var);
+        double x_kurtosis = (1-6*d.p()*(1-d.p())) / x_var;
         assert(std::abs(mean - x_mean) / x_mean < 0.01);
         assert(std::abs(var - x_var) / x_var < 0.01);
+        assert(std::abs(skew - x_skew) / x_skew < 0.01);
+        assert(std::abs(kurtosis - x_kurtosis) / x_kurtosis < 0.01);
     }
     {
         typedef std::binomial_distribution<> D;
-        typedef std::minstd_rand G;
+        typedef std::mt19937 G;
         G g;
         D d(30, .03125);
         const int N = 100000;
         std::vector<D::result_type> u;
         for (int i = 0; i < N; ++i)
-            u.push_back(d(g));
+        {
+            D::result_type v = d(g);
+            assert(d.min() <= v && v <= d.max());
+            u.push_back(v);
+        }
         double mean = std::accumulate(u.begin(), u.end(),
                                               double(0)) / u.size();
         double var = 0;
+        double skew = 0;
+        double kurtosis = 0;
         for (int i = 0; i < u.size(); ++i)
-            var += sqr(u[i] - mean);
+        {
+            double d = (u[i] - mean);
+            double d2 = sqr(d);
+            var += d2;
+            skew += d * d2;
+            kurtosis += d2 * d2;
+        }
         var /= u.size();
+        double dev = std::sqrt(var);
+        skew /= u.size() * dev * var;
+        kurtosis /= u.size() * var * var;
+        kurtosis -= 3;
         double x_mean = d.t() * d.p();
         double x_var = x_mean*(1-d.p());
+        double x_skew = (1-2*d.p()) / std::sqrt(x_var);
+        double x_kurtosis = (1-6*d.p()*(1-d.p())) / x_var;
         assert(std::abs(mean - x_mean) / x_mean < 0.01);
         assert(std::abs(var - x_var) / x_var < 0.01);
+        assert(std::abs(skew - x_skew) / x_skew < 0.01);
+        assert(std::abs(kurtosis - x_kurtosis) / x_kurtosis < 0.01);
     }
     {
         typedef std::binomial_distribution<> D;
-        typedef std::minstd_rand G;
+        typedef std::mt19937 G;
         G g;
         D d(40, .25);
         const int N = 100000;
         std::vector<D::result_type> u;
         for (int i = 0; i < N; ++i)
-            u.push_back(d(g));
+        {
+            D::result_type v = d(g);
+            assert(d.min() <= v && v <= d.max());
+            u.push_back(v);
+        }
         double mean = std::accumulate(u.begin(), u.end(),
                                               double(0)) / u.size();
         double var = 0;
+        double skew = 0;
+        double kurtosis = 0;
         for (int i = 0; i < u.size(); ++i)
-            var += sqr(u[i] - mean);
+        {
+            double d = (u[i] - mean);
+            double d2 = sqr(d);
+            var += d2;
+            skew += d * d2;
+            kurtosis += d2 * d2;
+        }
         var /= u.size();
+        double dev = std::sqrt(var);
+        skew /= u.size() * dev * var;
+        kurtosis /= u.size() * var * var;
+        kurtosis -= 3;
         double x_mean = d.t() * d.p();
         double x_var = x_mean*(1-d.p());
+        double x_skew = (1-2*d.p()) / std::sqrt(x_var);
+        double x_kurtosis = (1-6*d.p()*(1-d.p())) / x_var;
         assert(std::abs(mean - x_mean) / x_mean < 0.01);
         assert(std::abs(var - x_var) / x_var < 0.01);
+        assert(std::abs(skew - x_skew) / x_skew < 0.03);
+        assert(std::abs(kurtosis - x_kurtosis) / x_kurtosis < 0.01);
     }
     {
         typedef std::binomial_distribution<> D;
-        typedef std::minstd_rand G;
+        typedef std::mt19937 G;
         G g;
         D d(40, 0);
         const int N = 100000;
         std::vector<D::result_type> u;
         for (int i = 0; i < N; ++i)
-            u.push_back(d(g));
+        {
+            D::result_type v = d(g);
+            assert(d.min() <= v && v <= d.max());
+            u.push_back(v);
+        }
         double mean = std::accumulate(u.begin(), u.end(),
                                               double(0)) / u.size();
         double var = 0;
+        double skew = 0;
+        double kurtosis = 0;
         for (int i = 0; i < u.size(); ++i)
-            var += sqr(u[i] - mean);
+        {
+            double d = (u[i] - mean);
+            double d2 = sqr(d);
+            var += d2;
+            skew += d * d2;
+            kurtosis += d2 * d2;
+        }
         var /= u.size();
+        double dev = std::sqrt(var);
+        skew /= u.size() * dev * var;
+        kurtosis /= u.size() * var * var;
+        kurtosis -= 3;
         double x_mean = d.t() * d.p();
         double x_var = x_mean*(1-d.p());
+        double x_skew = (1-2*d.p()) / std::sqrt(x_var);
+        double x_kurtosis = (1-6*d.p()*(1-d.p())) / x_var;
         assert(mean == x_mean);
         assert(var == x_var);
     }
     {
         typedef std::binomial_distribution<> D;
-        typedef std::minstd_rand G;
+        typedef std::mt19937 G;
         G g;
         D d(40, 1);
         const int N = 100000;
         std::vector<D::result_type> u;
         for (int i = 0; i < N; ++i)
-            u.push_back(d(g));
+        {
+            D::result_type v = d(g);
+            assert(d.min() <= v && v <= d.max());
+            u.push_back(v);
+        }
         double mean = std::accumulate(u.begin(), u.end(),
                                               double(0)) / u.size();
         double var = 0;
+        double skew = 0;
+        double kurtosis = 0;
         for (int i = 0; i < u.size(); ++i)
-            var += sqr(u[i] - mean);
+        {
+            double d = (u[i] - mean);
+            double d2 = sqr(d);
+            var += d2;
+            skew += d * d2;
+            kurtosis += d2 * d2;
+        }
         var /= u.size();
+        double dev = std::sqrt(var);
+        skew /= u.size() * dev * var;
+        kurtosis /= u.size() * var * var;
+        kurtosis -= 3;
         double x_mean = d.t() * d.p();
         double x_var = x_mean*(1-d.p());
+        double x_skew = (1-2*d.p()) / std::sqrt(x_var);
+        double x_kurtosis = (1-6*d.p()*(1-d.p())) / x_var;
         assert(mean == x_mean);
         assert(var == x_var);
     }
     {
         typedef std::binomial_distribution<> D;
-        typedef std::minstd_rand G;
+        typedef std::mt19937 G;
         G g;
         D d(400, 0.5);
         const int N = 100000;
         std::vector<D::result_type> u;
         for (int i = 0; i < N; ++i)
-            u.push_back(d(g));
+        {
+            D::result_type v = d(g);
+            assert(d.min() <= v && v <= d.max());
+            u.push_back(v);
+        }
         double mean = std::accumulate(u.begin(), u.end(),
                                               double(0)) / u.size();
         double var = 0;
+        double skew = 0;
+        double kurtosis = 0;
         for (int i = 0; i < u.size(); ++i)
-            var += sqr(u[i] - mean);
+        {
+            double d = (u[i] - mean);
+            double d2 = sqr(d);
+            var += d2;
+            skew += d * d2;
+            kurtosis += d2 * d2;
+        }
         var /= u.size();
+        double dev = std::sqrt(var);
+        skew /= u.size() * dev * var;
+        kurtosis /= u.size() * var * var;
+        kurtosis -= 3;
         double x_mean = d.t() * d.p();
         double x_var = x_mean*(1-d.p());
+        double x_skew = (1-2*d.p()) / std::sqrt(x_var);
+        double x_kurtosis = (1-6*d.p()*(1-d.p())) / x_var;
         assert(std::abs(mean - x_mean) / x_mean < 0.01);
         assert(std::abs(var - x_var) / x_var < 0.01);
+        assert(std::abs(skew - x_skew) < 0.01);
+        assert(std::abs(kurtosis - x_kurtosis) / x_kurtosis < 0.01);
     }
     {
         typedef std::binomial_distribution<> D;
-        typedef std::minstd_rand G;
+        typedef std::mt19937 G;
         G g;
         D d(1, 0.5);
         const int N = 100000;
         std::vector<D::result_type> u;
         for (int i = 0; i < N; ++i)
-            u.push_back(d(g));
+        {
+            D::result_type v = d(g);
+            assert(d.min() <= v && v <= d.max());
+            u.push_back(v);
+        }
         double mean = std::accumulate(u.begin(), u.end(),
                                               double(0)) / u.size();
         double var = 0;
+        double skew = 0;
+        double kurtosis = 0;
         for (int i = 0; i < u.size(); ++i)
-            var += sqr(u[i] - mean);
+        {
+            double d = (u[i] - mean);
+            double d2 = sqr(d);
+            var += d2;
+            skew += d * d2;
+            kurtosis += d2 * d2;
+        }
         var /= u.size();
+        double dev = std::sqrt(var);
+        skew /= u.size() * dev * var;
+        kurtosis /= u.size() * var * var;
+        kurtosis -= 3;
         double x_mean = d.t() * d.p();
         double x_var = x_mean*(1-d.p());
+        double x_skew = (1-2*d.p()) / std::sqrt(x_var);
+        double x_kurtosis = (1-6*d.p()*(1-d.p())) / x_var;
         assert(std::abs(mean - x_mean) / x_mean < 0.01);
         assert(std::abs(var - x_var) / x_var < 0.01);
+        assert(std::abs(skew - x_skew) < 0.01);
+        assert(std::abs(kurtosis - x_kurtosis) / x_kurtosis < 0.01);
     }
     {
         typedef std::binomial_distribution<> D;
-        typedef std::minstd_rand G;
+        typedef std::mt19937 G;
         G g;
         D d(0, 0.005);
         const int N = 100000;
         std::vector<D::result_type> u;
         for (int i = 0; i < N; ++i)
-            u.push_back(d(g));
+        {
+            D::result_type v = d(g);
+            assert(d.min() <= v && v <= d.max());
+            u.push_back(v);
+        }
         double mean = std::accumulate(u.begin(), u.end(),
                                               double(0)) / u.size();
         double var = 0;
+        double skew = 0;
+        double kurtosis = 0;
         for (int i = 0; i < u.size(); ++i)
-            var += sqr(u[i] - mean);
+        {
+            double d = (u[i] - mean);
+            double d2 = sqr(d);
+            var += d2;
+            skew += d * d2;
+            kurtosis += d2 * d2;
+        }
         var /= u.size();
+        double dev = std::sqrt(var);
+        skew /= u.size() * dev * var;
+        kurtosis /= u.size() * var * var;
+        kurtosis -= 3;
         double x_mean = d.t() * d.p();
         double x_var = x_mean*(1-d.p());
+        double x_skew = (1-2*d.p()) / std::sqrt(x_var);
+        double x_kurtosis = (1-6*d.p()*(1-d.p())) / x_var;
         assert(mean == x_mean);
         assert(var == x_var);
     }
     {
         typedef std::binomial_distribution<> D;
-        typedef std::minstd_rand G;
+        typedef std::mt19937 G;
         G g;
         D d(0, 0);
         const int N = 100000;
         std::vector<D::result_type> u;
         for (int i = 0; i < N; ++i)
-            u.push_back(d(g));
+        {
+            D::result_type v = d(g);
+            assert(d.min() <= v && v <= d.max());
+            u.push_back(v);
+        }
         double mean = std::accumulate(u.begin(), u.end(),
                                               double(0)) / u.size();
         double var = 0;
+        double skew = 0;
+        double kurtosis = 0;
         for (int i = 0; i < u.size(); ++i)
-            var += sqr(u[i] - mean);
+        {
+            double d = (u[i] - mean);
+            double d2 = sqr(d);
+            var += d2;
+            skew += d * d2;
+            kurtosis += d2 * d2;
+        }
         var /= u.size();
+        double dev = std::sqrt(var);
+        skew /= u.size() * dev * var;
+        kurtosis /= u.size() * var * var;
+        kurtosis -= 3;
         double x_mean = d.t() * d.p();
         double x_var = x_mean*(1-d.p());
+        double x_skew = (1-2*d.p()) / std::sqrt(x_var);
+        double x_kurtosis = (1-6*d.p()*(1-d.p())) / x_var;
         assert(mean == x_mean);
         assert(var == x_var);
     }
     {
         typedef std::binomial_distribution<> D;
-        typedef std::minstd_rand G;
+        typedef std::mt19937 G;
         G g;
         D d(0, 1);
         const int N = 100000;
         std::vector<D::result_type> u;
         for (int i = 0; i < N; ++i)
-            u.push_back(d(g));
+        {
+            D::result_type v = d(g);
+            assert(d.min() <= v && v <= d.max());
+            u.push_back(v);
+        }
         double mean = std::accumulate(u.begin(), u.end(),
                                               double(0)) / u.size();
         double var = 0;
+        double skew = 0;
+        double kurtosis = 0;
         for (int i = 0; i < u.size(); ++i)
-            var += sqr(u[i] - mean);
+        {
+            double d = (u[i] - mean);
+            double d2 = sqr(d);
+            var += d2;
+            skew += d * d2;
+            kurtosis += d2 * d2;
+        }
         var /= u.size();
+        double dev = std::sqrt(var);
+        skew /= u.size() * dev * var;
+        kurtosis /= u.size() * var * var;
+        kurtosis -= 3;
         double x_mean = d.t() * d.p();
         double x_var = x_mean*(1-d.p());
+        double x_skew = (1-2*d.p()) / std::sqrt(x_var);
+        double x_kurtosis = (1-6*d.p()*(1-d.p())) / x_var;
         assert(mean == x_mean);
         assert(var == x_var);
     }
diff --git a/test/numerics/rand/rand.dis/rand.dist.bern/rand.dist.bern.bin/eval_param.pass.cpp b/test/numerics/rand/rand.dis/rand.dist.bern/rand.dist.bern.bin/eval_param.pass.cpp
index 4bb9ea1..2632ab5 100644
--- a/test/numerics/rand/rand.dis/rand.dist.bern/rand.dist.bern.bin/eval_param.pass.cpp
+++ b/test/numerics/rand/rand.dis/rand.dist.bern/rand.dist.bern.bin/eval_param.pass.cpp
@@ -39,60 +39,120 @@
         const int N = 100000;
         std::vector<D::result_type> u;
         for (int i = 0; i < N; ++i)
-            u.push_back(d(g, p));
+        {
+            D::result_type v = d(g, p);
+            assert(0 <= v && v <= p.t());
+            u.push_back(v);
+        }
         double mean = std::accumulate(u.begin(), u.end(),
                                               double(0)) / u.size();
         double var = 0;
+        double skew = 0;
+        double kurtosis = 0;
         for (int i = 0; i < u.size(); ++i)
-            var += sqr(u[i] - mean);
+        {
+            double d = (u[i] - mean);
+            double d2 = sqr(d);
+            var += d2;
+            skew += d * d2;
+            kurtosis += d2 * d2;
+        }
         var /= u.size();
+        double dev = std::sqrt(var);
+        skew /= u.size() * dev * var;
+        kurtosis /= u.size() * var * var;
+        kurtosis -= 3;
         double x_mean = p.t() * p.p();
         double x_var = x_mean*(1-p.p());
+        double x_skew = (1-2*p.p()) / std::sqrt(x_var);
+        double x_kurtosis = (1-6*p.p()*(1-p.p())) / x_var;
         assert(std::abs(mean - x_mean) / x_mean < 0.01);
         assert(std::abs(var - x_var) / x_var < 0.01);
+        assert(std::abs(skew - x_skew) / x_skew < 0.01);
+        assert(std::abs(kurtosis - x_kurtosis) / x_kurtosis < 0.01);
     }
     {
         typedef std::binomial_distribution<> D;
         typedef D::param_type P;
-        typedef std::minstd_rand G;
+        typedef std::mt19937 G;
         G g;
         D d(16, .75);
         P p(30, .03125);
         const int N = 100000;
         std::vector<D::result_type> u;
         for (int i = 0; i < N; ++i)
-            u.push_back(d(g, p));
+        {
+            D::result_type v = d(g, p);
+            assert(0 <= v && v <= p.t());
+            u.push_back(v);
+        }
         double mean = std::accumulate(u.begin(), u.end(),
                                               double(0)) / u.size();
         double var = 0;
+        double skew = 0;
+        double kurtosis = 0;
         for (int i = 0; i < u.size(); ++i)
-            var += sqr(u[i] - mean);
+        {
+            double d = (u[i] - mean);
+            double d2 = sqr(d);
+            var += d2;
+            skew += d * d2;
+            kurtosis += d2 * d2;
+        }
         var /= u.size();
+        double dev = std::sqrt(var);
+        skew /= u.size() * dev * var;
+        kurtosis /= u.size() * var * var;
+        kurtosis -= 3;
         double x_mean = p.t() * p.p();
         double x_var = x_mean*(1-p.p());
+        double x_skew = (1-2*p.p()) / std::sqrt(x_var);
+        double x_kurtosis = (1-6*p.p()*(1-p.p())) / x_var;
         assert(std::abs(mean - x_mean) / x_mean < 0.01);
         assert(std::abs(var - x_var) / x_var < 0.01);
+        assert(std::abs(skew - x_skew) / x_skew < 0.01);
+        assert(std::abs(kurtosis - x_kurtosis) / x_kurtosis < 0.01);
     }
     {
         typedef std::binomial_distribution<> D;
         typedef D::param_type P;
-        typedef std::minstd_rand G;
+        typedef std::mt19937 G;
         G g;
         D d(16, .75);
         P p(40, .25);
         const int N = 100000;
         std::vector<D::result_type> u;
         for (int i = 0; i < N; ++i)
-            u.push_back(d(g, p));
+        {
+            D::result_type v = d(g, p);
+            assert(0 <= v && v <= p.t());
+            u.push_back(v);
+        }
         double mean = std::accumulate(u.begin(), u.end(),
                                               double(0)) / u.size();
         double var = 0;
+        double skew = 0;
+        double kurtosis = 0;
         for (int i = 0; i < u.size(); ++i)
-            var += sqr(u[i] - mean);
+        {
+            double d = (u[i] - mean);
+            double d2 = sqr(d);
+            var += d2;
+            skew += d * d2;
+            kurtosis += d2 * d2;
+        }
         var /= u.size();
+        double dev = std::sqrt(var);
+        skew /= u.size() * dev * var;
+        kurtosis /= u.size() * var * var;
+        kurtosis -= 3;
         double x_mean = p.t() * p.p();
         double x_var = x_mean*(1-p.p());
+        double x_skew = (1-2*p.p()) / std::sqrt(x_var);
+        double x_kurtosis = (1-6*p.p()*(1-p.p())) / x_var;
         assert(std::abs(mean - x_mean) / x_mean < 0.01);
         assert(std::abs(var - x_var) / x_var < 0.01);
+        assert(std::abs(skew - x_skew) / x_skew < 0.03);
+        assert(std::abs(kurtosis - x_kurtosis) / x_kurtosis < 0.01);
     }
 }
diff --git a/test/numerics/rand/rand.dis/rand.dist.norm/rand.dist.norm.chisq/eval.pass.cpp b/test/numerics/rand/rand.dis/rand.dist.norm/rand.dist.norm.chisq/eval.pass.cpp
index 5c24f91..5713643 100644
--- a/test/numerics/rand/rand.dis/rand.dist.norm/rand.dist.norm.chisq/eval.pass.cpp
+++ b/test/numerics/rand/rand.dis/rand.dist.norm/rand.dist.norm.chisq/eval.pass.cpp
@@ -35,20 +35,39 @@
         typedef std::minstd_rand G;
         G g;
         D d(0.5);
-        const int N = 100000;
+        const int N = 1000000;
         std::vector<D::result_type> u;
         for (int i = 0; i < N; ++i)
-            u.push_back(d(g));
-        D::result_type mean = std::accumulate(u.begin(), u.end(),
-                                              D::result_type(0)) / u.size();
-        D::result_type var = 0;
+        {
+            D::result_type v = d(g);
+            assert(d.min() < v);
+            u.push_back(v);
+        }
+        double mean = std::accumulate(u.begin(), u.end(), 0.0) / u.size();
+        double var = 0;
+        double skew = 0;
+        double kurtosis = 0;
         for (int i = 0; i < u.size(); ++i)
-            var += sqr(u[i] - mean);
+        {
+            double d = (u[i] - mean);
+            double d2 = sqr(d);
+            var += d2;
+            skew += d * d2;
+            kurtosis += d2 * d2;
+        }
         var /= u.size();
-        D::result_type x_mean = d.n();
-        D::result_type x_var = 2*d.n();
-        assert(std::abs(mean - x_mean) / x_mean < 0.02);
-        assert(std::abs(var - x_var) / x_var < 0.02);
+        double dev = std::sqrt(var);
+        skew /= u.size() * dev * var;
+        kurtosis /= u.size() * var * var;
+        kurtosis -= 3;
+        double x_mean = d.n();
+        double x_var = 2 * d.n();
+        double x_skew = std::sqrt(8 / d.n());
+        double x_kurtosis = 12 / d.n();
+        assert(std::abs(mean - x_mean) / x_mean < 0.01);
+        assert(std::abs(var - x_var) / x_var < 0.01);
+        assert(std::abs(skew - x_skew) / x_skew < 0.01);
+        assert(std::abs(kurtosis - x_kurtosis) / x_kurtosis < 0.01);
     }
     {
         typedef std::chi_squared_distribution<> D;
@@ -56,40 +75,78 @@
         typedef std::minstd_rand G;
         G g;
         D d(1);
-        const int N = 100000;
+        const int N = 1000000;
         std::vector<D::result_type> u;
         for (int i = 0; i < N; ++i)
-            u.push_back(d(g));
-        D::result_type mean = std::accumulate(u.begin(), u.end(),
-                                              D::result_type(0)) / u.size();
-        D::result_type var = 0;
+        {
+            D::result_type v = d(g);
+            assert(d.min() < v);
+            u.push_back(v);
+        }
+        double mean = std::accumulate(u.begin(), u.end(), 0.0) / u.size();
+        double var = 0;
+        double skew = 0;
+        double kurtosis = 0;
         for (int i = 0; i < u.size(); ++i)
-            var += sqr(u[i] - mean);
+        {
+            double d = (u[i] - mean);
+            double d2 = sqr(d);
+            var += d2;
+            skew += d * d2;
+            kurtosis += d2 * d2;
+        }
         var /= u.size();
-        D::result_type x_mean = d.n();
-        D::result_type x_var = 2*d.n();
-        assert(std::abs(mean - x_mean) / x_mean < 0.02);
-        assert(std::abs(var - x_var) / x_var < 0.02);
+        double dev = std::sqrt(var);
+        skew /= u.size() * dev * var;
+        kurtosis /= u.size() * var * var;
+        kurtosis -= 3;
+        double x_mean = d.n();
+        double x_var = 2 * d.n();
+        double x_skew = std::sqrt(8 / d.n());
+        double x_kurtosis = 12 / d.n();
+        assert(std::abs(mean - x_mean) / x_mean < 0.01);
+        assert(std::abs(var - x_var) / x_var < 0.01);
+        assert(std::abs(skew - x_skew) / x_skew < 0.01);
+        assert(std::abs(kurtosis - x_kurtosis) / x_kurtosis < 0.01);
     }
     {
         typedef std::chi_squared_distribution<> D;
         typedef D::param_type P;
-        typedef std::minstd_rand G;
+        typedef std::mt19937 G;
         G g;
         D d(2);
-        const int N = 100000;
+        const int N = 1000000;
         std::vector<D::result_type> u;
         for (int i = 0; i < N; ++i)
-            u.push_back(d(g));
-        D::result_type mean = std::accumulate(u.begin(), u.end(),
-                                              D::result_type(0)) / u.size();
-        D::result_type var = 0;
+        {
+            D::result_type v = d(g);
+            assert(d.min() < v);
+            u.push_back(v);
+        }
+        double mean = std::accumulate(u.begin(), u.end(), 0.0) / u.size();
+        double var = 0;
+        double skew = 0;
+        double kurtosis = 0;
         for (int i = 0; i < u.size(); ++i)
-            var += sqr(u[i] - mean);
+        {
+            double d = (u[i] - mean);
+            double d2 = sqr(d);
+            var += d2;
+            skew += d * d2;
+            kurtosis += d2 * d2;
+        }
         var /= u.size();
-        D::result_type x_mean = d.n();
-        D::result_type x_var = 2*d.n();
-        assert(std::abs(mean - x_mean) / x_mean < 0.02);
-        assert(std::abs(var - x_var) / x_var < 0.02);
+        double dev = std::sqrt(var);
+        skew /= u.size() * dev * var;
+        kurtosis /= u.size() * var * var;
+        kurtosis -= 3;
+        double x_mean = d.n();
+        double x_var = 2 * d.n();
+        double x_skew = std::sqrt(8 / d.n());
+        double x_kurtosis = 12 / d.n();
+        assert(std::abs(mean - x_mean) / x_mean < 0.01);
+        assert(std::abs(var - x_var) / x_var < 0.01);
+        assert(std::abs(skew - x_skew) / x_skew < 0.01);
+        assert(std::abs(kurtosis - x_kurtosis) / x_kurtosis < 0.01);
     }
 }
diff --git a/test/numerics/rand/rand.dis/rand.dist.norm/rand.dist.norm.chisq/eval_param.pass.cpp b/test/numerics/rand/rand.dis/rand.dist.norm/rand.dist.norm.chisq/eval_param.pass.cpp
index b771225..561f94d 100644
--- a/test/numerics/rand/rand.dis/rand.dist.norm/rand.dist.norm.chisq/eval_param.pass.cpp
+++ b/test/numerics/rand/rand.dis/rand.dist.norm/rand.dist.norm.chisq/eval_param.pass.cpp
@@ -36,42 +36,80 @@
         G g;
         D d(0.5);
         P p(1);
-        const int N = 100000;
+        const int N = 1000000;
         std::vector<D::result_type> u;
         for (int i = 0; i < N; ++i)
-            u.push_back(d(g, p));
-        D::result_type mean = std::accumulate(u.begin(), u.end(),
-                                              D::result_type(0)) / u.size();
-        D::result_type var = 0;
+        {
+            D::result_type v = d(g, p);
+            assert(d.min() < v);
+            u.push_back(v);
+        }
+        double mean = std::accumulate(u.begin(), u.end(), 0.0) / u.size();
+        double var = 0;
+        double skew = 0;
+        double kurtosis = 0;
         for (int i = 0; i < u.size(); ++i)
-            var += sqr(u[i] - mean);
+        {
+            double d = (u[i] - mean);
+            double d2 = sqr(d);
+            var += d2;
+            skew += d * d2;
+            kurtosis += d2 * d2;
+        }
         var /= u.size();
-        D::result_type x_mean = p.n();
-        D::result_type x_var = 2*p.n();
-        assert(std::abs(mean - x_mean) / x_mean < 0.02);
-        assert(std::abs(var - x_var) / x_var < 0.02);
+        double dev = std::sqrt(var);
+        skew /= u.size() * dev * var;
+        kurtosis /= u.size() * var * var;
+        kurtosis -= 3;
+        double x_mean = p.n();
+        double x_var = 2 * p.n();
+        double x_skew = std::sqrt(8 / p.n());
+        double x_kurtosis = 12 / p.n();
+        assert(std::abs(mean - x_mean) / x_mean < 0.01);
+        assert(std::abs(var - x_var) / x_var < 0.01);
+        assert(std::abs(skew - x_skew) / x_skew < 0.01);
+        assert(std::abs(kurtosis - x_kurtosis) / x_kurtosis < 0.01);
     }
     {
         typedef std::chi_squared_distribution<> D;
         typedef D::param_type P;
-        typedef std::minstd_rand G;
+        typedef std::mt19937 G;
         G g;
         D d(1);
         P p(2);
-        const int N = 100000;
+        const int N = 1000000;
         std::vector<D::result_type> u;
         for (int i = 0; i < N; ++i)
-            u.push_back(d(g, p));
-        D::result_type mean = std::accumulate(u.begin(), u.end(),
-                                              D::result_type(0)) / u.size();
-        D::result_type var = 0;
+        {
+            D::result_type v = d(g, p);
+            assert(d.min() < v);
+            u.push_back(v);
+        }
+        double mean = std::accumulate(u.begin(), u.end(), 0.0) / u.size();
+        double var = 0;
+        double skew = 0;
+        double kurtosis = 0;
         for (int i = 0; i < u.size(); ++i)
-            var += sqr(u[i] - mean);
+        {
+            double d = (u[i] - mean);
+            double d2 = sqr(d);
+            var += d2;
+            skew += d * d2;
+            kurtosis += d2 * d2;
+        }
         var /= u.size();
-        D::result_type x_mean = p.n();
-        D::result_type x_var = 2*p.n();
-        assert(std::abs(mean - x_mean) / x_mean < 0.02);
-        assert(std::abs(var - x_var) / x_var < 0.02);
+        double dev = std::sqrt(var);
+        skew /= u.size() * dev * var;
+        kurtosis /= u.size() * var * var;
+        kurtosis -= 3;
+        double x_mean = p.n();
+        double x_var = 2 * p.n();
+        double x_skew = std::sqrt(8 / p.n());
+        double x_kurtosis = 12 / p.n();
+        assert(std::abs(mean - x_mean) / x_mean < 0.01);
+        assert(std::abs(var - x_var) / x_var < 0.01);
+        assert(std::abs(skew - x_skew) / x_skew < 0.01);
+        assert(std::abs(kurtosis - x_kurtosis) / x_kurtosis < 0.01);
     }
     {
         typedef std::chi_squared_distribution<> D;
@@ -80,19 +118,38 @@
         G g;
         D d(2);
         P p(.5);
-        const int N = 100000;
+        const int N = 1000000;
         std::vector<D::result_type> u;
         for (int i = 0; i < N; ++i)
-            u.push_back(d(g, p));
-        D::result_type mean = std::accumulate(u.begin(), u.end(),
-                                              D::result_type(0)) / u.size();
-        D::result_type var = 0;
+        {
+            D::result_type v = d(g, p);
+            assert(d.min() < v);
+            u.push_back(v);
+        }
+        double mean = std::accumulate(u.begin(), u.end(), 0.0) / u.size();
+        double var = 0;
+        double skew = 0;
+        double kurtosis = 0;
         for (int i = 0; i < u.size(); ++i)
-            var += sqr(u[i] - mean);
+        {
+            double d = (u[i] - mean);
+            double d2 = sqr(d);
+            var += d2;
+            skew += d * d2;
+            kurtosis += d2 * d2;
+        }
         var /= u.size();
-        D::result_type x_mean = p.n();
-        D::result_type x_var = 2*p.n();
-        assert(std::abs(mean - x_mean) / x_mean < 0.02);
-        assert(std::abs(var - x_var) / x_var < 0.02);
+        double dev = std::sqrt(var);
+        skew /= u.size() * dev * var;
+        kurtosis /= u.size() * var * var;
+        kurtosis -= 3;
+        double x_mean = p.n();
+        double x_var = 2 * p.n();
+        double x_skew = std::sqrt(8 / p.n());
+        double x_kurtosis = 12 / p.n();
+        assert(std::abs(mean - x_mean) / x_mean < 0.01);
+        assert(std::abs(var - x_var) / x_var < 0.01);
+        assert(std::abs(skew - x_skew) / x_skew < 0.01);
+        assert(std::abs(kurtosis - x_kurtosis) / x_kurtosis < 0.01);
     }
 }
diff --git a/test/numerics/rand/rand.dis/rand.dist.norm/rand.dist.norm.normal/eval.pass.cpp b/test/numerics/rand/rand.dis/rand.dist.norm/rand.dist.norm.normal/eval.pass.cpp
index 71f9014..abfb8aa 100644
--- a/test/numerics/rand/rand.dis/rand.dist.norm/rand.dist.norm.normal/eval.pass.cpp
+++ b/test/numerics/rand/rand.dis/rand.dist.norm/rand.dist.norm.normal/eval.pass.cpp
@@ -35,19 +35,34 @@
         typedef std::minstd_rand G;
         G g;
         D d(5, 4);
-        const int N = 10000;
+        const int N = 1000000;
         std::vector<D::result_type> u;
         for (int i = 0; i < N; ++i)
             u.push_back(d(g));
-        D::result_type mean = std::accumulate(u.begin(), u.end(),
-                                              D::result_type(0)) / u.size();
-        D::result_type var = 0;
+        double mean = std::accumulate(u.begin(), u.end(), 0.0) / u.size();
+        double var = 0;
+        double skew = 0;
+        double kurtosis = 0;
         for (int i = 0; i < u.size(); ++i)
-            var += sqr(u[i] - mean);
+        {
+            double d = (u[i] - mean);
+            double d2 = sqr(d);
+            var += d2;
+            skew += d * d2;
+            kurtosis += d2 * d2;
+        }
         var /= u.size();
-        D::result_type x_mean = d.mean();
-        D::result_type x_var = sqr(d.stddev());
+        double dev = std::sqrt(var);
+        skew /= u.size() * dev * var;
+        kurtosis /= u.size() * var * var;
+        kurtosis -= 3;
+        double x_mean = d.mean();
+        double x_var = sqr(d.stddev());
+        double x_skew = 0;
+        double x_kurtosis = 0;
         assert(std::abs(mean - x_mean) / x_mean < 0.01);
-        assert(std::abs(var - x_var) / x_var < 0.02);
+        assert(std::abs(var - x_var) / x_var < 0.01);
+        assert(std::abs(skew - x_skew) < 0.01);
+        assert(std::abs(kurtosis - x_kurtosis) < 0.01);
     }
 }
diff --git a/test/numerics/rand/rand.dis/rand.dist.norm/rand.dist.norm.normal/eval_param.pass.cpp b/test/numerics/rand/rand.dis/rand.dist.norm/rand.dist.norm.normal/eval_param.pass.cpp
index 652c4c1..e55c75f 100644
--- a/test/numerics/rand/rand.dis/rand.dist.norm/rand.dist.norm.normal/eval_param.pass.cpp
+++ b/test/numerics/rand/rand.dis/rand.dist.norm/rand.dist.norm.normal/eval_param.pass.cpp
@@ -36,19 +36,34 @@
         G g;
         D d(5, 4);
         P p(50, .5);
-        const int N = 1000;
+        const int N = 1000000;
         std::vector<D::result_type> u;
         for (int i = 0; i < N; ++i)
             u.push_back(d(g, p));
-        D::result_type mean = std::accumulate(u.begin(), u.end(),
-                                              D::result_type(0)) / u.size();
-        D::result_type var = 0;
+        double mean = std::accumulate(u.begin(), u.end(), 0.0) / u.size();
+        double var = 0;
+        double skew = 0;
+        double kurtosis = 0;
         for (int i = 0; i < u.size(); ++i)
-            var += sqr(u[i] - mean);
+        {
+            double d = (u[i] - mean);
+            double d2 = sqr(d);
+            var += d2;
+            skew += d * d2;
+            kurtosis += d2 * d2;
+        }
         var /= u.size();
-        D::result_type x_mean = p.mean();
-        D::result_type x_var = sqr(p.stddev());
+        double dev = std::sqrt(var);
+        skew /= u.size() * dev * var;
+        kurtosis /= u.size() * var * var;
+        kurtosis -= 3;
+        double x_mean = p.mean();
+        double x_var = sqr(p.stddev());
+        double x_skew = 0;
+        double x_kurtosis = 0;
         assert(std::abs(mean - x_mean) / x_mean < 0.01);
         assert(std::abs(var - x_var) / x_var < 0.01);
+        assert(std::abs(skew - x_skew) < 0.01);
+        assert(std::abs(kurtosis - x_kurtosis) < 0.01);
     }
 }
diff --git a/test/numerics/rand/rand.dis/rand.dist.pois/rand.dist.pois.exp/eval.pass.cpp b/test/numerics/rand/rand.dis/rand.dist.pois/rand.dist.pois.exp/eval.pass.cpp
index 48cfc7e..b4e8928 100644
--- a/test/numerics/rand/rand.dis/rand.dist.pois/rand.dist.pois.exp/eval.pass.cpp
+++ b/test/numerics/rand/rand.dis/rand.dist.pois/rand.dist.pois.exp/eval.pass.cpp
@@ -32,22 +32,121 @@
     {
         typedef std::exponential_distribution<> D;
         typedef D::param_type P;
-        typedef std::knuth_b G;
+        typedef std::mt19937 G;
         G g;
         D d(.75);
-        const int N = 1000;
+        const int N = 1000000;
         std::vector<D::result_type> u;
         for (int i = 0; i < N; ++i)
-            u.push_back(d(g));
-        D::result_type mean = std::accumulate(u.begin(), u.end(),
-                                              D::result_type(0)) / u.size();
-        D::result_type var = 0;
+        {
+            D::result_type v = d(g);
+            assert(d.min() < v);
+            u.push_back(v);
+        }
+        double mean = std::accumulate(u.begin(), u.end(), 0.0) / u.size();
+        double var = 0;
+        double skew = 0;
+        double kurtosis = 0;
         for (int i = 0; i < u.size(); ++i)
-            var += sqr(u[i] - mean);
+        {
+            double d = (u[i] - mean);
+            double d2 = sqr(d);
+            var += d2;
+            skew += d * d2;
+            kurtosis += d2 * d2;
+        }
         var /= u.size();
-        D::result_type x_mean = 1/d.lambda();
-        D::result_type x_var = 1/sqr(d.lambda());
+        double dev = std::sqrt(var);
+        skew /= u.size() * dev * var;
+        kurtosis /= u.size() * var * var;
+        kurtosis -= 3;
+        double x_mean = 1/d.lambda();
+        double x_var = 1/sqr(d.lambda());
+        double x_skew = 2;
+        double x_kurtosis = 6;
         assert(std::abs(mean - x_mean) / x_mean < 0.01);
         assert(std::abs(var - x_var) / x_var < 0.01);
+        assert(std::abs(skew - x_skew) / x_skew < 0.01);
+        assert(std::abs(kurtosis - x_kurtosis) / x_kurtosis < 0.01);
+    }
+    {
+        typedef std::exponential_distribution<> D;
+        typedef D::param_type P;
+        typedef std::mt19937 G;
+        G g;
+        D d(1);
+        const int N = 1000000;
+        std::vector<D::result_type> u;
+        for (int i = 0; i < N; ++i)
+        {
+            D::result_type v = d(g);
+            assert(d.min() < v);
+            u.push_back(v);
+        }
+        double mean = std::accumulate(u.begin(), u.end(), 0.0) / u.size();
+        double var = 0;
+        double skew = 0;
+        double kurtosis = 0;
+        for (int i = 0; i < u.size(); ++i)
+        {
+            double d = (u[i] - mean);
+            double d2 = sqr(d);
+            var += d2;
+            skew += d * d2;
+            kurtosis += d2 * d2;
+        }
+        var /= u.size();
+        double dev = std::sqrt(var);
+        skew /= u.size() * dev * var;
+        kurtosis /= u.size() * var * var;
+        kurtosis -= 3;
+        double x_mean = 1/d.lambda();
+        double x_var = 1/sqr(d.lambda());
+        double x_skew = 2;
+        double x_kurtosis = 6;
+        assert(std::abs(mean - x_mean) / x_mean < 0.01);
+        assert(std::abs(var - x_var) / x_var < 0.01);
+        assert(std::abs(skew - x_skew) / x_skew < 0.01);
+        assert(std::abs(kurtosis - x_kurtosis) / x_kurtosis < 0.01);
+    }
+    {
+        typedef std::exponential_distribution<> D;
+        typedef D::param_type P;
+        typedef std::mt19937 G;
+        G g;
+        D d(10);
+        const int N = 1000000;
+        std::vector<D::result_type> u;
+        for (int i = 0; i < N; ++i)
+        {
+            D::result_type v = d(g);
+            assert(d.min() < v);
+            u.push_back(v);
+        }
+        double mean = std::accumulate(u.begin(), u.end(), 0.0) / u.size();
+        double var = 0;
+        double skew = 0;
+        double kurtosis = 0;
+        for (int i = 0; i < u.size(); ++i)
+        {
+            double d = (u[i] - mean);
+            double d2 = sqr(d);
+            var += d2;
+            skew += d * d2;
+            kurtosis += d2 * d2;
+        }
+        var /= u.size();
+        double dev = std::sqrt(var);
+        skew /= u.size() * dev * var;
+        kurtosis /= u.size() * var * var;
+        kurtosis -= 3;
+        double x_mean = 1/d.lambda();
+        double x_var = 1/sqr(d.lambda());
+        double x_skew = 2;
+        double x_kurtosis = 6;
+        assert(std::abs(mean - x_mean) / x_mean < 0.01);
+        assert(std::abs(var - x_var) / x_var < 0.01);
+        assert(std::abs(skew - x_skew) / x_skew < 0.01);
+        assert(std::abs(kurtosis - x_kurtosis) / x_kurtosis < 0.01);
     }
 }
diff --git a/test/numerics/rand/rand.dis/rand.dist.pois/rand.dist.pois.exp/eval_param.pass.cpp b/test/numerics/rand/rand.dis/rand.dist.pois/rand.dist.pois.exp/eval_param.pass.cpp
index e51b45e..e5805f1 100644
--- a/test/numerics/rand/rand.dis/rand.dist.pois/rand.dist.pois.exp/eval_param.pass.cpp
+++ b/test/numerics/rand/rand.dis/rand.dist.pois/rand.dist.pois.exp/eval_param.pass.cpp
@@ -32,23 +32,42 @@
     {
         typedef std::exponential_distribution<> D;
         typedef D::param_type P;
-        typedef std::knuth_b G;
+        typedef std::mt19937 G;
         G g;
         D d(.75);
         P p(2);
-        const int N = 1000;
+        const int N = 1000000;
         std::vector<D::result_type> u;
         for (int i = 0; i < N; ++i)
-            u.push_back(d(g, p));
-        D::result_type mean = std::accumulate(u.begin(), u.end(),
-                                              D::result_type(0)) / u.size();
-        D::result_type var = 0;
+        {
+            D::result_type v = d(g, p);
+            assert(d.min() < v);
+            u.push_back(v);
+        }
+        double mean = std::accumulate(u.begin(), u.end(), 0.0) / u.size();
+        double var = 0;
+        double skew = 0;
+        double kurtosis = 0;
         for (int i = 0; i < u.size(); ++i)
-            var += sqr(u[i] - mean);
+        {
+            double d = (u[i] - mean);
+            double d2 = sqr(d);
+            var += d2;
+            skew += d * d2;
+            kurtosis += d2 * d2;
+        }
         var /= u.size();
-        D::result_type x_mean = 1/p.lambda();
-        D::result_type x_var = 1/sqr(p.lambda());
+        double dev = std::sqrt(var);
+        skew /= u.size() * dev * var;
+        kurtosis /= u.size() * var * var;
+        kurtosis -= 3;
+        double x_mean = 1/p.lambda();
+        double x_var = 1/sqr(p.lambda());
+        double x_skew = 2;
+        double x_kurtosis = 6;
         assert(std::abs(mean - x_mean) / x_mean < 0.01);
         assert(std::abs(var - x_var) / x_var < 0.01);
+        assert(std::abs(skew - x_skew) / x_skew < 0.01);
+        assert(std::abs(kurtosis - x_kurtosis) / x_kurtosis < 0.01);
     }
 }
diff --git a/test/numerics/rand/rand.dis/rand.dist.pois/rand.dist.pois.exp/max.pass.cpp b/test/numerics/rand/rand.dis/rand.dist.pois/rand.dist.pois.exp/max.pass.cpp
index 004b878..c80c74f 100644
--- a/test/numerics/rand/rand.dis/rand.dist.pois/rand.dist.pois.exp/max.pass.cpp
+++ b/test/numerics/rand/rand.dis/rand.dist.pois/rand.dist.pois.exp/max.pass.cpp
@@ -23,6 +23,6 @@
         typedef std::exponential_distribution<> D;
         D d(.25);
         D::result_type m = d.max();
-        assert(146 < m && m < 147);
+        assert(m == std::numeric_limits<D::result_type>::infinity());
     }
 }
diff --git a/test/numerics/rand/rand.dis/rand.dist.pois/rand.dist.pois.gamma/eval.pass.cpp b/test/numerics/rand/rand.dis/rand.dist.pois/rand.dist.pois.gamma/eval.pass.cpp
index 13f7f36..0d803f3 100644
--- a/test/numerics/rand/rand.dis/rand.dist.pois/rand.dist.pois.gamma/eval.pass.cpp
+++ b/test/numerics/rand/rand.dis/rand.dist.pois/rand.dist.pois.gamma/eval.pass.cpp
@@ -32,64 +32,121 @@
     {
         typedef std::gamma_distribution<> D;
         typedef D::param_type P;
-        typedef std::minstd_rand G;
+        typedef std::mt19937 G;
         G g;
         D d(0.5, 2);
-        const int N = 100000;
+        const int N = 1000000;
         std::vector<D::result_type> u;
         for (int i = 0; i < N; ++i)
-            u.push_back(d(g));
-        D::result_type mean = std::accumulate(u.begin(), u.end(),
-                                              D::result_type(0)) / u.size();
-        D::result_type var = 0;
+        {
+            D::result_type v = d(g);
+            assert(d.min() < v);
+            u.push_back(v);
+        }
+        double mean = std::accumulate(u.begin(), u.end(), 0.0) / u.size();
+        double var = 0;
+        double skew = 0;
+        double kurtosis = 0;
         for (int i = 0; i < u.size(); ++i)
-            var += sqr(u[i] - mean);
+        {
+            double d = (u[i] - mean);
+            double d2 = sqr(d);
+            var += d2;
+            skew += d * d2;
+            kurtosis += d2 * d2;
+        }
         var /= u.size();
-        D::result_type x_mean = d.alpha() * d.beta();
-        D::result_type x_var = d.alpha() * sqr(d.beta());
-        assert(std::abs(mean - x_mean) / x_mean < 0.02);
-        assert(std::abs(var - x_var) / x_var < 0.02);
+        double dev = std::sqrt(var);
+        skew /= u.size() * dev * var;
+        kurtosis /= u.size() * var * var;
+        kurtosis -= 3;
+        double x_mean = d.alpha() * d.beta();
+        double x_var = d.alpha() * sqr(d.beta());
+        double x_skew = 2 / std::sqrt(d.alpha());
+        double x_kurtosis = 6 / d.alpha();
+        assert(std::abs(mean - x_mean) / x_mean < 0.01);
+        assert(std::abs(var - x_var) / x_var < 0.01);
+        assert(std::abs(skew - x_skew) / x_skew < 0.01);
+        assert(std::abs(kurtosis - x_kurtosis) / x_kurtosis < 0.01);
     }
     {
         typedef std::gamma_distribution<> D;
         typedef D::param_type P;
-        typedef std::minstd_rand G;
+        typedef std::mt19937 G;
         G g;
         D d(1, .5);
-        const int N = 100000;
+        const int N = 1000000;
         std::vector<D::result_type> u;
         for (int i = 0; i < N; ++i)
-            u.push_back(d(g));
-        D::result_type mean = std::accumulate(u.begin(), u.end(),
-                                              D::result_type(0)) / u.size();
-        D::result_type var = 0;
+        {
+            D::result_type v = d(g);
+            assert(d.min() < v);
+            u.push_back(v);
+        }
+        double mean = std::accumulate(u.begin(), u.end(), 0.0) / u.size();
+        double var = 0;
+        double skew = 0;
+        double kurtosis = 0;
         for (int i = 0; i < u.size(); ++i)
-            var += sqr(u[i] - mean);
+        {
+            double d = (u[i] - mean);
+            double d2 = sqr(d);
+            var += d2;
+            skew += d * d2;
+            kurtosis += d2 * d2;
+        }
         var /= u.size();
-        D::result_type x_mean = d.alpha() * d.beta();
-        D::result_type x_var = d.alpha() * sqr(d.beta());
-        assert(std::abs(mean - x_mean) / x_mean < 0.02);
-        assert(std::abs(var - x_var) / x_var < 0.02);
+        double dev = std::sqrt(var);
+        skew /= u.size() * dev * var;
+        kurtosis /= u.size() * var * var;
+        kurtosis -= 3;
+        double x_mean = d.alpha() * d.beta();
+        double x_var = d.alpha() * sqr(d.beta());
+        double x_skew = 2 / std::sqrt(d.alpha());
+        double x_kurtosis = 6 / d.alpha();
+        assert(std::abs(mean - x_mean) / x_mean < 0.01);
+        assert(std::abs(var - x_var) / x_var < 0.01);
+        assert(std::abs(skew - x_skew) / x_skew < 0.01);
+        assert(std::abs(kurtosis - x_kurtosis) / x_kurtosis < 0.01);
     }
     {
         typedef std::gamma_distribution<> D;
         typedef D::param_type P;
-        typedef std::minstd_rand G;
+        typedef std::mt19937 G;
         G g;
         D d(2, 3);
-        const int N = 100000;
+        const int N = 1000000;
         std::vector<D::result_type> u;
         for (int i = 0; i < N; ++i)
-            u.push_back(d(g));
-        D::result_type mean = std::accumulate(u.begin(), u.end(),
-                                              D::result_type(0)) / u.size();
-        D::result_type var = 0;
+        {
+            D::result_type v = d(g);
+            assert(d.min() < v);
+            u.push_back(v);
+        }
+        double mean = std::accumulate(u.begin(), u.end(), 0.0) / u.size();
+        double var = 0;
+        double skew = 0;
+        double kurtosis = 0;
         for (int i = 0; i < u.size(); ++i)
-            var += sqr(u[i] - mean);
+        {
+            double d = (u[i] - mean);
+            double d2 = sqr(d);
+            var += d2;
+            skew += d * d2;
+            kurtosis += d2 * d2;
+        }
         var /= u.size();
-        D::result_type x_mean = d.alpha() * d.beta();
-        D::result_type x_var = d.alpha() * sqr(d.beta());
-        assert(std::abs(mean - x_mean) / x_mean < 0.02);
-        assert(std::abs(var - x_var) / x_var < 0.02);
+        double dev = std::sqrt(var);
+        skew /= u.size() * dev * var;
+        kurtosis /= u.size() * var * var;
+        kurtosis -= 3;
+        double x_mean = d.alpha() * d.beta();
+        double x_var = d.alpha() * sqr(d.beta());
+        double x_skew = 2 / std::sqrt(d.alpha());
+        double x_kurtosis = 6 / d.alpha();
+        assert(std::abs(mean - x_mean) / x_mean < 0.01);
+        assert(std::abs(var - x_var) / x_var < 0.01);
+        assert(std::abs(skew - x_skew) / x_skew < 0.01);
+        assert(std::abs(kurtosis - x_kurtosis) / x_kurtosis < 0.01);
     }
 }
diff --git a/test/numerics/rand/rand.dis/rand.dist.pois/rand.dist.pois.gamma/eval_param.pass.cpp b/test/numerics/rand/rand.dis/rand.dist.pois/rand.dist.pois.gamma/eval_param.pass.cpp
index 1ab0531..cdb3fc5 100644
--- a/test/numerics/rand/rand.dis/rand.dist.pois/rand.dist.pois.gamma/eval_param.pass.cpp
+++ b/test/numerics/rand/rand.dis/rand.dist.pois/rand.dist.pois.gamma/eval_param.pass.cpp
@@ -32,67 +32,124 @@
     {
         typedef std::gamma_distribution<> D;
         typedef D::param_type P;
-        typedef std::minstd_rand G;
+        typedef std::mt19937 G;
         G g;
         D d(0.5, 2);
         P p(1, .5);
-        const int N = 100000;
+        const int N = 1000000;
         std::vector<D::result_type> u;
         for (int i = 0; i < N; ++i)
-            u.push_back(d(g, p));
-        D::result_type mean = std::accumulate(u.begin(), u.end(),
-                                              D::result_type(0)) / u.size();
-        D::result_type var = 0;
+        {
+            D::result_type v = d(g, p);
+            assert(d.min() < v);
+            u.push_back(v);
+        }
+        double mean = std::accumulate(u.begin(), u.end(), 0.0) / u.size();
+        double var = 0;
+        double skew = 0;
+        double kurtosis = 0;
         for (int i = 0; i < u.size(); ++i)
-            var += sqr(u[i] - mean);
+        {
+            double d = (u[i] - mean);
+            double d2 = sqr(d);
+            var += d2;
+            skew += d * d2;
+            kurtosis += d2 * d2;
+        }
         var /= u.size();
-        D::result_type x_mean = p.alpha() * p.beta();
-        D::result_type x_var = p.alpha() * sqr(p.beta());
-        assert(std::abs(mean - x_mean) / x_mean < 0.02);
-        assert(std::abs(var - x_var) / x_var < 0.02);
+        double dev = std::sqrt(var);
+        skew /= u.size() * dev * var;
+        kurtosis /= u.size() * var * var;
+        kurtosis -= 3;
+        double x_mean = p.alpha() * p.beta();
+        double x_var = p.alpha() * sqr(p.beta());
+        double x_skew = 2 / std::sqrt(p.alpha());
+        double x_kurtosis = 6 / p.alpha();
+        assert(std::abs(mean - x_mean) / x_mean < 0.01);
+        assert(std::abs(var - x_var) / x_var < 0.01);
+        assert(std::abs(skew - x_skew) / x_skew < 0.01);
+        assert(std::abs(kurtosis - x_kurtosis) / x_kurtosis < 0.01);
     }
     {
         typedef std::gamma_distribution<> D;
         typedef D::param_type P;
-        typedef std::minstd_rand G;
+        typedef std::mt19937 G;
         G g;
         D d(1, .5);
         P p(2, 3);
-        const int N = 100000;
+        const int N = 1000000;
         std::vector<D::result_type> u;
         for (int i = 0; i < N; ++i)
-            u.push_back(d(g, p));
-        D::result_type mean = std::accumulate(u.begin(), u.end(),
-                                              D::result_type(0)) / u.size();
-        D::result_type var = 0;
+        {
+            D::result_type v = d(g, p);
+            assert(d.min() < v);
+            u.push_back(v);
+        }
+        double mean = std::accumulate(u.begin(), u.end(), 0.0) / u.size();
+        double var = 0;
+        double skew = 0;
+        double kurtosis = 0;
         for (int i = 0; i < u.size(); ++i)
-            var += sqr(u[i] - mean);
+        {
+            double d = (u[i] - mean);
+            double d2 = sqr(d);
+            var += d2;
+            skew += d * d2;
+            kurtosis += d2 * d2;
+        }
         var /= u.size();
-        D::result_type x_mean = p.alpha() * p.beta();
-        D::result_type x_var = p.alpha() * sqr(p.beta());
-        assert(std::abs(mean - x_mean) / x_mean < 0.02);
-        assert(std::abs(var - x_var) / x_var < 0.02);
+        double dev = std::sqrt(var);
+        skew /= u.size() * dev * var;
+        kurtosis /= u.size() * var * var;
+        kurtosis -= 3;
+        double x_mean = p.alpha() * p.beta();
+        double x_var = p.alpha() * sqr(p.beta());
+        double x_skew = 2 / std::sqrt(p.alpha());
+        double x_kurtosis = 6 / p.alpha();
+        assert(std::abs(mean - x_mean) / x_mean < 0.01);
+        assert(std::abs(var - x_var) / x_var < 0.01);
+        assert(std::abs(skew - x_skew) / x_skew < 0.01);
+        assert(std::abs(kurtosis - x_kurtosis) / x_kurtosis < 0.01);
     }
     {
         typedef std::gamma_distribution<> D;
         typedef D::param_type P;
-        typedef std::minstd_rand G;
+        typedef std::mt19937 G;
         G g;
         D d(2, 3);
         P p(.5, 2);
-        const int N = 100000;
+        const int N = 1000000;
         std::vector<D::result_type> u;
         for (int i = 0; i < N; ++i)
-            u.push_back(d(g, p));
-        D::result_type mean = std::accumulate(u.begin(), u.end(),
-                                              D::result_type(0)) / u.size();
-        D::result_type var = 0;
+        {
+            D::result_type v = d(g, p);
+            assert(d.min() < v);
+            u.push_back(v);
+        }
+        double mean = std::accumulate(u.begin(), u.end(), 0.0) / u.size();
+        double var = 0;
+        double skew = 0;
+        double kurtosis = 0;
         for (int i = 0; i < u.size(); ++i)
-            var += sqr(u[i] - mean);
+        {
+            double d = (u[i] - mean);
+            double d2 = sqr(d);
+            var += d2;
+            skew += d * d2;
+            kurtosis += d2 * d2;
+        }
         var /= u.size();
-        D::result_type x_mean = p.alpha() * p.beta();
-        D::result_type x_var = p.alpha() * sqr(p.beta());
-        assert(std::abs(mean - x_mean) / x_mean < 0.02);
-        assert(std::abs(var - x_var) / x_var < 0.02);
+        double dev = std::sqrt(var);
+        skew /= u.size() * dev * var;
+        kurtosis /= u.size() * var * var;
+        kurtosis -= 3;
+        double x_mean = p.alpha() * p.beta();
+        double x_var = p.alpha() * sqr(p.beta());
+        double x_skew = 2 / std::sqrt(p.alpha());
+        double x_kurtosis = 6 / p.alpha();
+        assert(std::abs(mean - x_mean) / x_mean < 0.01);
+        assert(std::abs(var - x_var) / x_var < 0.01);
+        assert(std::abs(skew - x_skew) / x_skew < 0.01);
+        assert(std::abs(kurtosis - x_kurtosis) / x_kurtosis < 0.01);
     }
 }
diff --git a/test/numerics/rand/rand.dis/rand.dist.pois/rand.dist.pois.poisson/eval.pass.cpp b/test/numerics/rand/rand.dis/rand.dist.pois/rand.dist.pois.poisson/eval.pass.cpp
index 4c8a1fd..2bc5c23 100644
--- a/test/numerics/rand/rand.dis/rand.dist.pois/rand.dist.pois.poisson/eval.pass.cpp
+++ b/test/numerics/rand/rand.dis/rand.dist.pois/rand.dist.pois.poisson/eval.pass.cpp
@@ -37,16 +37,36 @@
         const int N = 100000;
         std::vector<double> u;
         for (int i = 0; i < N; ++i)
-            u.push_back(d(g));
+        {
+            D::result_type v = d(g);
+            assert(d.min() <= v && v <= d.max());
+            u.push_back(v);
+        }
         double mean = std::accumulate(u.begin(), u.end(), 0.0) / u.size();
         double var = 0;
+        double skew = 0;
+        double kurtosis = 0;
         for (int i = 0; i < u.size(); ++i)
-            var += sqr(u[i] - mean);
+        {
+            double d = (u[i] - mean);
+            double d2 = sqr(d);
+            var += d2;
+            skew += d * d2;
+            kurtosis += d2 * d2;
+        }
         var /= u.size();
+        double dev = std::sqrt(var);
+        skew /= u.size() * dev * var;
+        kurtosis /= u.size() * var * var;
+        kurtosis -= 3;
         double x_mean = d.mean();
         double x_var = d.mean();
+        double x_skew = 1 / std::sqrt(x_var);
+        double x_kurtosis = 1 / x_var;
         assert(std::abs(mean - x_mean) / x_mean < 0.01);
         assert(std::abs(var - x_var) / x_var < 0.01);
+        assert(std::abs(skew - x_skew) / x_skew < 0.01);
+        assert(std::abs(kurtosis - x_kurtosis) / x_kurtosis < 0.03);
     }
     {
         typedef std::poisson_distribution<> D;
@@ -56,34 +76,74 @@
         const int N = 100000;
         std::vector<double> u;
         for (int i = 0; i < N; ++i)
-            u.push_back(d(g));
+        {
+            D::result_type v = d(g);
+            assert(d.min() <= v && v <= d.max());
+            u.push_back(v);
+        }
         double mean = std::accumulate(u.begin(), u.end(), 0.0) / u.size();
         double var = 0;
+        double skew = 0;
+        double kurtosis = 0;
         for (int i = 0; i < u.size(); ++i)
-            var += sqr(u[i] - mean);
+        {
+            double d = (u[i] - mean);
+            double d2 = sqr(d);
+            var += d2;
+            skew += d * d2;
+            kurtosis += d2 * d2;
+        }
         var /= u.size();
+        double dev = std::sqrt(var);
+        skew /= u.size() * dev * var;
+        kurtosis /= u.size() * var * var;
+        kurtosis -= 3;
         double x_mean = d.mean();
         double x_var = d.mean();
+        double x_skew = 1 / std::sqrt(x_var);
+        double x_kurtosis = 1 / x_var;
         assert(std::abs(mean - x_mean) / x_mean < 0.01);
         assert(std::abs(var - x_var) / x_var < 0.01);
+        assert(std::abs(skew - x_skew) / x_skew < 0.01);
+        assert(std::abs(kurtosis - x_kurtosis) / x_kurtosis < 0.04);
     }
     {
         typedef std::poisson_distribution<> D;
-        typedef std::minstd_rand G;
+        typedef std::mt19937 G;
         G g;
         D d(20);
-        const int N = 100000;
+        const int N = 1000000;
         std::vector<double> u;
         for (int i = 0; i < N; ++i)
-            u.push_back(d(g));
+        {
+            D::result_type v = d(g);
+            assert(d.min() <= v && v <= d.max());
+            u.push_back(v);
+        }
         double mean = std::accumulate(u.begin(), u.end(), 0.0) / u.size();
         double var = 0;
+        double skew = 0;
+        double kurtosis = 0;
         for (int i = 0; i < u.size(); ++i)
-            var += sqr(u[i] - mean);
+        {
+            double d = (u[i] - mean);
+            double d2 = sqr(d);
+            var += d2;
+            skew += d * d2;
+            kurtosis += d2 * d2;
+        }
         var /= u.size();
+        double dev = std::sqrt(var);
+        skew /= u.size() * dev * var;
+        kurtosis /= u.size() * var * var;
+        kurtosis -= 3;
         double x_mean = d.mean();
         double x_var = d.mean();
+        double x_skew = 1 / std::sqrt(x_var);
+        double x_kurtosis = 1 / x_var;
         assert(std::abs(mean - x_mean) / x_mean < 0.01);
         assert(std::abs(var - x_var) / x_var < 0.01);
+        assert(std::abs(skew - x_skew) / x_skew < 0.01);
+        assert(std::abs(kurtosis - x_kurtosis) / x_kurtosis < 0.01);
     }
 }
diff --git a/test/numerics/rand/rand.dis/rand.dist.pois/rand.dist.pois.poisson/eval_param.pass.cpp b/test/numerics/rand/rand.dis/rand.dist.pois/rand.dist.pois.poisson/eval_param.pass.cpp
index 0c094bd..ba22729 100644
--- a/test/numerics/rand/rand.dis/rand.dist.pois/rand.dist.pois.poisson/eval_param.pass.cpp
+++ b/test/numerics/rand/rand.dis/rand.dist.pois/rand.dist.pois.poisson/eval_param.pass.cpp
@@ -39,16 +39,36 @@
         const int N = 100000;
         std::vector<double> u;
         for (int i = 0; i < N; ++i)
-            u.push_back(d(g, p));
+        {
+            D::result_type v = d(g, p);
+            assert(d.min() <= v && v <= d.max());
+            u.push_back(v);
+        }
         double mean = std::accumulate(u.begin(), u.end(), 0.0) / u.size();
         double var = 0;
+        double skew = 0;
+        double kurtosis = 0;
         for (int i = 0; i < u.size(); ++i)
-            var += sqr(u[i] - mean);
+        {
+            double d = (u[i] - mean);
+            double d2 = sqr(d);
+            var += d2;
+            skew += d * d2;
+            kurtosis += d2 * d2;
+        }
         var /= u.size();
+        double dev = std::sqrt(var);
+        skew /= u.size() * dev * var;
+        kurtosis /= u.size() * var * var;
+        kurtosis -= 3;
         double x_mean = p.mean();
         double x_var = p.mean();
+        double x_skew = 1 / std::sqrt(x_var);
+        double x_kurtosis = 1 / x_var;
         assert(std::abs(mean - x_mean) / x_mean < 0.01);
         assert(std::abs(var - x_var) / x_var < 0.01);
+        assert(std::abs(skew - x_skew) / x_skew < 0.01);
+        assert(std::abs(kurtosis - x_kurtosis) / x_kurtosis < 0.03);
     }
     {
         typedef std::poisson_distribution<> D;
@@ -60,36 +80,76 @@
         const int N = 100000;
         std::vector<double> u;
         for (int i = 0; i < N; ++i)
-            u.push_back(d(g, p));
+        {
+            D::result_type v = d(g, p);
+            assert(d.min() <= v && v <= d.max());
+            u.push_back(v);
+        }
         double mean = std::accumulate(u.begin(), u.end(), 0.0) / u.size();
         double var = 0;
+        double skew = 0;
+        double kurtosis = 0;
         for (int i = 0; i < u.size(); ++i)
-            var += sqr(u[i] - mean);
+        {
+            double d = (u[i] - mean);
+            double d2 = sqr(d);
+            var += d2;
+            skew += d * d2;
+            kurtosis += d2 * d2;
+        }
         var /= u.size();
+        double dev = std::sqrt(var);
+        skew /= u.size() * dev * var;
+        kurtosis /= u.size() * var * var;
+        kurtosis -= 3;
         double x_mean = p.mean();
         double x_var = p.mean();
+        double x_skew = 1 / std::sqrt(x_var);
+        double x_kurtosis = 1 / x_var;
         assert(std::abs(mean - x_mean) / x_mean < 0.01);
         assert(std::abs(var - x_var) / x_var < 0.01);
+        assert(std::abs(skew - x_skew) / x_skew < 0.01);
+        assert(std::abs(kurtosis - x_kurtosis) / x_kurtosis < 0.04);
     }
     {
         typedef std::poisson_distribution<> D;
         typedef D::param_type P;
-        typedef std::minstd_rand G;
+        typedef std::mt19937 G;
         G g;
         D d(2);
         P p(20);
-        const int N = 100000;
+        const int N = 1000000;
         std::vector<double> u;
         for (int i = 0; i < N; ++i)
-            u.push_back(d(g, p));
+        {
+            D::result_type v = d(g, p);
+            assert(d.min() <= v && v <= d.max());
+            u.push_back(v);
+        }
         double mean = std::accumulate(u.begin(), u.end(), 0.0) / u.size();
         double var = 0;
+        double skew = 0;
+        double kurtosis = 0;
         for (int i = 0; i < u.size(); ++i)
-            var += sqr(u[i] - mean);
+        {
+            double d = (u[i] - mean);
+            double d2 = sqr(d);
+            var += d2;
+            skew += d * d2;
+            kurtosis += d2 * d2;
+        }
         var /= u.size();
+        double dev = std::sqrt(var);
+        skew /= u.size() * dev * var;
+        kurtosis /= u.size() * var * var;
+        kurtosis -= 3;
         double x_mean = p.mean();
         double x_var = p.mean();
+        double x_skew = 1 / std::sqrt(x_var);
+        double x_kurtosis = 1 / x_var;
         assert(std::abs(mean - x_mean) / x_mean < 0.01);
         assert(std::abs(var - x_var) / x_var < 0.01);
+        assert(std::abs(skew - x_skew) / x_skew < 0.01);
+        assert(std::abs(kurtosis - x_kurtosis) / x_kurtosis < 0.01);
     }
 }
diff --git a/test/numerics/rand/rand.dis/rand.dist.pois/rand.dist.pois.weibull/eval.pass.cpp b/test/numerics/rand/rand.dis/rand.dist.pois/rand.dist.pois.weibull/eval.pass.cpp
index 27e2b52..1381a43 100644
--- a/test/numerics/rand/rand.dis/rand.dist.pois/rand.dist.pois.weibull/eval.pass.cpp
+++ b/test/numerics/rand/rand.dis/rand.dist.pois/rand.dist.pois.weibull/eval.pass.cpp
@@ -35,20 +35,43 @@
         typedef std::mt19937 G;
         G g;
         D d(0.5, 2);
-        const int N = 100000;
+        const int N = 1000000;
         std::vector<D::result_type> u;
         for (int i = 0; i < N; ++i)
-            u.push_back(d(g));
-        D::result_type mean = std::accumulate(u.begin(), u.end(),
-                                              D::result_type(0)) / u.size();
-        D::result_type var = 0;
+        {
+            D::result_type v = d(g);
+            assert(d.min() <= v);
+            u.push_back(v);
+        }
+        double mean = std::accumulate(u.begin(), u.end(), 0.0) / u.size();
+        double var = 0;
+        double skew = 0;
+        double kurtosis = 0;
         for (int i = 0; i < u.size(); ++i)
-            var += sqr(u[i] - mean);
+        {
+            double d = (u[i] - mean);
+            double d2 = sqr(d);
+            var += d2;
+            skew += d * d2;
+            kurtosis += d2 * d2;
+        }
         var /= u.size();
-        D::result_type x_mean = d.b() * std::tgamma(1 + 1/d.a());
-        D::result_type x_var = sqr(d.b()) * std::tgamma(1 + 2/d.a()) - sqr(x_mean);
-        assert(std::abs(mean - x_mean) / x_mean < 0.02);
-        assert(std::abs(var - x_var) / x_var < 0.02);
+        double dev = std::sqrt(var);
+        skew /= u.size() * dev * var;
+        kurtosis /= u.size() * var * var;
+        kurtosis -= 3;
+        double x_mean = d.b() * std::tgamma(1 + 1/d.a());
+        double x_var = sqr(d.b()) * std::tgamma(1 + 2/d.a()) - sqr(x_mean);
+        double x_skew = (sqr(d.b())*d.b() * std::tgamma(1 + 3/d.a()) -
+                        3*x_mean*x_var - sqr(x_mean)*x_mean) /
+                        (std::sqrt(x_var)*x_var);
+        double x_kurtosis = (sqr(sqr(d.b())) * std::tgamma(1 + 4/d.a()) -
+                       4*x_skew*x_var*sqrt(x_var)*x_mean -
+                       6*sqr(x_mean)*x_var - sqr(sqr(x_mean))) / sqr(x_var) - 3;
+        assert(std::abs(mean - x_mean) / x_mean < 0.01);
+        assert(std::abs(var - x_var) / x_var < 0.01);
+        assert(std::abs(skew - x_skew) / x_skew < 0.01);
+        assert(std::abs(kurtosis - x_kurtosis) / x_kurtosis < 0.03);
     }
     {
         typedef std::weibull_distribution<> D;
@@ -56,20 +79,43 @@
         typedef std::mt19937 G;
         G g;
         D d(1, .5);
-        const int N = 100000;
+        const int N = 1000000;
         std::vector<D::result_type> u;
         for (int i = 0; i < N; ++i)
-            u.push_back(d(g));
-        D::result_type mean = std::accumulate(u.begin(), u.end(),
-                                              D::result_type(0)) / u.size();
-        D::result_type var = 0;
+        {
+            D::result_type v = d(g);
+            assert(d.min() <= v);
+            u.push_back(v);
+        }
+        double mean = std::accumulate(u.begin(), u.end(), 0.0) / u.size();
+        double var = 0;
+        double skew = 0;
+        double kurtosis = 0;
         for (int i = 0; i < u.size(); ++i)
-            var += sqr(u[i] - mean);
+        {
+            double d = (u[i] - mean);
+            double d2 = sqr(d);
+            var += d2;
+            skew += d * d2;
+            kurtosis += d2 * d2;
+        }
         var /= u.size();
-        D::result_type x_mean = d.b() * std::tgamma(1 + 1/d.a());
-        D::result_type x_var = sqr(d.b()) * std::tgamma(1 + 2/d.a()) - sqr(x_mean);
-        assert(std::abs(mean - x_mean) / x_mean < 0.02);
-        assert(std::abs(var - x_var) / x_var < 0.02);
+        double dev = std::sqrt(var);
+        skew /= u.size() * dev * var;
+        kurtosis /= u.size() * var * var;
+        kurtosis -= 3;
+        double x_mean = d.b() * std::tgamma(1 + 1/d.a());
+        double x_var = sqr(d.b()) * std::tgamma(1 + 2/d.a()) - sqr(x_mean);
+        double x_skew = (sqr(d.b())*d.b() * std::tgamma(1 + 3/d.a()) -
+                        3*x_mean*x_var - sqr(x_mean)*x_mean) /
+                        (std::sqrt(x_var)*x_var);
+        double x_kurtosis = (sqr(sqr(d.b())) * std::tgamma(1 + 4/d.a()) -
+                       4*x_skew*x_var*sqrt(x_var)*x_mean -
+                       6*sqr(x_mean)*x_var - sqr(sqr(x_mean))) / sqr(x_var) - 3;
+        assert(std::abs(mean - x_mean) / x_mean < 0.01);
+        assert(std::abs(var - x_var) / x_var < 0.01);
+        assert(std::abs(skew - x_skew) / x_skew < 0.01);
+        assert(std::abs(kurtosis - x_kurtosis) / x_kurtosis < 0.01);
     }
     {
         typedef std::weibull_distribution<> D;
@@ -77,19 +123,42 @@
         typedef std::mt19937 G;
         G g;
         D d(2, 3);
-        const int N = 100000;
+        const int N = 1000000;
         std::vector<D::result_type> u;
         for (int i = 0; i < N; ++i)
-            u.push_back(d(g));
-        D::result_type mean = std::accumulate(u.begin(), u.end(),
-                                              D::result_type(0)) / u.size();
-        D::result_type var = 0;
+        {
+            D::result_type v = d(g);
+            assert(d.min() <= v);
+            u.push_back(v);
+        }
+        double mean = std::accumulate(u.begin(), u.end(), 0.0) / u.size();
+        double var = 0;
+        double skew = 0;
+        double kurtosis = 0;
         for (int i = 0; i < u.size(); ++i)
-            var += sqr(u[i] - mean);
+        {
+            double d = (u[i] - mean);
+            double d2 = sqr(d);
+            var += d2;
+            skew += d * d2;
+            kurtosis += d2 * d2;
+        }
         var /= u.size();
-        D::result_type x_mean = d.b() * std::tgamma(1 + 1/d.a());
-        D::result_type x_var = sqr(d.b()) * std::tgamma(1 + 2/d.a()) - sqr(x_mean);
-        assert(std::abs(mean - x_mean) / x_mean < 0.02);
-        assert(std::abs(var - x_var) / x_var < 0.02);
+        double dev = std::sqrt(var);
+        skew /= u.size() * dev * var;
+        kurtosis /= u.size() * var * var;
+        kurtosis -= 3;
+        double x_mean = d.b() * std::tgamma(1 + 1/d.a());
+        double x_var = sqr(d.b()) * std::tgamma(1 + 2/d.a()) - sqr(x_mean);
+        double x_skew = (sqr(d.b())*d.b() * std::tgamma(1 + 3/d.a()) -
+                        3*x_mean*x_var - sqr(x_mean)*x_mean) /
+                        (std::sqrt(x_var)*x_var);
+        double x_kurtosis = (sqr(sqr(d.b())) * std::tgamma(1 + 4/d.a()) -
+                       4*x_skew*x_var*sqrt(x_var)*x_mean -
+                       6*sqr(x_mean)*x_var - sqr(sqr(x_mean))) / sqr(x_var) - 3;
+        assert(std::abs(mean - x_mean) / x_mean < 0.01);
+        assert(std::abs(var - x_var) / x_var < 0.01);
+        assert(std::abs(skew - x_skew) / x_skew < 0.01);
+        assert(std::abs(kurtosis - x_kurtosis) / x_kurtosis < 0.03);
     }
 }
diff --git a/test/numerics/rand/rand.dis/rand.dist.pois/rand.dist.pois.weibull/eval_param.pass.cpp b/test/numerics/rand/rand.dis/rand.dist.pois/rand.dist.pois.weibull/eval_param.pass.cpp
index 48610f8..3c1076c 100644
--- a/test/numerics/rand/rand.dis/rand.dist.pois/rand.dist.pois.weibull/eval_param.pass.cpp
+++ b/test/numerics/rand/rand.dis/rand.dist.pois/rand.dist.pois.weibull/eval_param.pass.cpp
@@ -36,20 +36,43 @@
         G g;
         D d(0.5, 2);
         P p(1, .5);
-        const int N = 100000;
+        const int N = 1000000;
         std::vector<D::result_type> u;
         for (int i = 0; i < N; ++i)
-            u.push_back(d(g, p));
-        D::result_type mean = std::accumulate(u.begin(), u.end(),
-                                              D::result_type(0)) / u.size();
-        D::result_type var = 0;
+        {
+            D::result_type v = d(g, p);
+            assert(d.min() <= v);
+            u.push_back(v);
+        }
+        double mean = std::accumulate(u.begin(), u.end(), 0.0) / u.size();
+        double var = 0;
+        double skew = 0;
+        double kurtosis = 0;
         for (int i = 0; i < u.size(); ++i)
-            var += sqr(u[i] - mean);
+        {
+            double d = (u[i] - mean);
+            double d2 = sqr(d);
+            var += d2;
+            skew += d * d2;
+            kurtosis += d2 * d2;
+        }
         var /= u.size();
-        D::result_type x_mean = p.b() * std::tgamma(1 + 1/p.a());
-        D::result_type x_var = sqr(p.b()) * std::tgamma(1 + 2/p.a()) - sqr(x_mean);
-        assert(std::abs(mean - x_mean) / x_mean < 0.02);
-        assert(std::abs(var - x_var) / x_var < 0.02);
+        double dev = std::sqrt(var);
+        skew /= u.size() * dev * var;
+        kurtosis /= u.size() * var * var;
+        kurtosis -= 3;
+        double x_mean = p.b() * std::tgamma(1 + 1/p.a());
+        double x_var = sqr(p.b()) * std::tgamma(1 + 2/p.a()) - sqr(x_mean);
+        double x_skew = (sqr(p.b())*p.b() * std::tgamma(1 + 3/p.a()) -
+                        3*x_mean*x_var - sqr(x_mean)*x_mean) /
+                        (std::sqrt(x_var)*x_var);
+        double x_kurtosis = (sqr(sqr(p.b())) * std::tgamma(1 + 4/p.a()) -
+                       4*x_skew*x_var*sqrt(x_var)*x_mean -
+                       6*sqr(x_mean)*x_var - sqr(sqr(x_mean))) / sqr(x_var) - 3;
+        assert(std::abs(mean - x_mean) / x_mean < 0.01);
+        assert(std::abs(var - x_var) / x_var < 0.01);
+        assert(std::abs(skew - x_skew) / x_skew < 0.01);
+        assert(std::abs(kurtosis - x_kurtosis) / x_kurtosis < 0.01);
     }
     {
         typedef std::weibull_distribution<> D;
@@ -58,20 +81,43 @@
         G g;
         D d(1, .5);
         P p(2, 3);
-        const int N = 100000;
+        const int N = 1000000;
         std::vector<D::result_type> u;
         for (int i = 0; i < N; ++i)
-            u.push_back(d(g, p));
-        D::result_type mean = std::accumulate(u.begin(), u.end(),
-                                              D::result_type(0)) / u.size();
-        D::result_type var = 0;
+        {
+            D::result_type v = d(g, p);
+            assert(d.min() <= v);
+            u.push_back(v);
+        }
+        double mean = std::accumulate(u.begin(), u.end(), 0.0) / u.size();
+        double var = 0;
+        double skew = 0;
+        double kurtosis = 0;
         for (int i = 0; i < u.size(); ++i)
-            var += sqr(u[i] - mean);
+        {
+            double d = (u[i] - mean);
+            double d2 = sqr(d);
+            var += d2;
+            skew += d * d2;
+            kurtosis += d2 * d2;
+        }
         var /= u.size();
-        D::result_type x_mean = p.b() * std::tgamma(1 + 1/p.a());
-        D::result_type x_var = sqr(p.b()) * std::tgamma(1 + 2/p.a()) - sqr(x_mean);
-        assert(std::abs(mean - x_mean) / x_mean < 0.02);
-        assert(std::abs(var - x_var) / x_var < 0.02);
+        double dev = std::sqrt(var);
+        skew /= u.size() * dev * var;
+        kurtosis /= u.size() * var * var;
+        kurtosis -= 3;
+        double x_mean = p.b() * std::tgamma(1 + 1/p.a());
+        double x_var = sqr(p.b()) * std::tgamma(1 + 2/p.a()) - sqr(x_mean);
+        double x_skew = (sqr(p.b())*p.b() * std::tgamma(1 + 3/p.a()) -
+                        3*x_mean*x_var - sqr(x_mean)*x_mean) /
+                        (std::sqrt(x_var)*x_var);
+        double x_kurtosis = (sqr(sqr(p.b())) * std::tgamma(1 + 4/p.a()) -
+                       4*x_skew*x_var*sqrt(x_var)*x_mean -
+                       6*sqr(x_mean)*x_var - sqr(sqr(x_mean))) / sqr(x_var) - 3;
+        assert(std::abs(mean - x_mean) / x_mean < 0.01);
+        assert(std::abs(var - x_var) / x_var < 0.01);
+        assert(std::abs(skew - x_skew) / x_skew < 0.01);
+        assert(std::abs(kurtosis - x_kurtosis) / x_kurtosis < 0.03);
     }
     {
         typedef std::weibull_distribution<> D;
@@ -80,19 +126,42 @@
         G g;
         D d(2, 3);
         P p(.5, 2);
-        const int N = 100000;
+        const int N = 1000000;
         std::vector<D::result_type> u;
         for (int i = 0; i < N; ++i)
-            u.push_back(d(g, p));
-        D::result_type mean = std::accumulate(u.begin(), u.end(),
-                                              D::result_type(0)) / u.size();
-        D::result_type var = 0;
+        {
+            D::result_type v = d(g, p);
+            assert(d.min() <= v);
+            u.push_back(v);
+        }
+        double mean = std::accumulate(u.begin(), u.end(), 0.0) / u.size();
+        double var = 0;
+        double skew = 0;
+        double kurtosis = 0;
         for (int i = 0; i < u.size(); ++i)
-            var += sqr(u[i] - mean);
+        {
+            double d = (u[i] - mean);
+            double d2 = sqr(d);
+            var += d2;
+            skew += d * d2;
+            kurtosis += d2 * d2;
+        }
         var /= u.size();
-        D::result_type x_mean = p.b() * std::tgamma(1 + 1/p.a());
-        D::result_type x_var = sqr(p.b()) * std::tgamma(1 + 2/p.a()) - sqr(x_mean);
-        assert(std::abs(mean - x_mean) / x_mean < 0.02);
-        assert(std::abs(var - x_var) / x_var < 0.02);
+        double dev = std::sqrt(var);
+        skew /= u.size() * dev * var;
+        kurtosis /= u.size() * var * var;
+        kurtosis -= 3;
+        double x_mean = p.b() * std::tgamma(1 + 1/p.a());
+        double x_var = sqr(p.b()) * std::tgamma(1 + 2/p.a()) - sqr(x_mean);
+        double x_skew = (sqr(p.b())*p.b() * std::tgamma(1 + 3/p.a()) -
+                        3*x_mean*x_var - sqr(x_mean)*x_mean) /
+                        (std::sqrt(x_var)*x_var);
+        double x_kurtosis = (sqr(sqr(p.b())) * std::tgamma(1 + 4/p.a()) -
+                       4*x_skew*x_var*sqrt(x_var)*x_mean -
+                       6*sqr(x_mean)*x_var - sqr(sqr(x_mean))) / sqr(x_var) - 3;
+        assert(std::abs(mean - x_mean) / x_mean < 0.01);
+        assert(std::abs(var - x_var) / x_var < 0.01);
+        assert(std::abs(skew - x_skew) / x_skew < 0.01);
+        assert(std::abs(kurtosis - x_kurtosis) / x_kurtosis < 0.03);
     }
 }
diff --git a/test/numerics/rand/rand.dis/rand.dist.uni/rand.dist.uni.int/eval.pass.cpp b/test/numerics/rand/rand.dis/rand.dist.uni/rand.dist.uni.int/eval.pass.cpp
index 192242a..3443e16 100644
--- a/test/numerics/rand/rand.dis/rand.dist.uni/rand.dist.uni.int/eval.pass.cpp
+++ b/test/numerics/rand/rand.dis/rand.dist.uni/rand.dist.uni.int/eval.pass.cpp
@@ -16,6 +16,16 @@
 
 #include <random>
 #include <cassert>
+#include <vector>
+#include <numeric>
+
+template <class T>
+inline
+T
+sqr(T x)
+{
+    return x * x;
+}
 
 int main()
 {
@@ -23,6 +33,375 @@
         typedef std::uniform_int_distribution<> D;
         typedef std::minstd_rand0 G;
         G g;
+        D d;
+        const int N = 100000;
+        std::vector<D::result_type> u;
+        for (int i = 0; i < N; ++i)
+        {
+            D::result_type v = d(g);
+            assert(d.a() <= v && v <= d.b());
+            u.push_back(v);
+        }
+        double mean = std::accumulate(u.begin(), u.end(),
+                                              double(0)) / u.size();
+        double var = 0;
+        double skew = 0;
+        double kurtosis = 0;
+        for (int i = 0; i < u.size(); ++i)
+        {
+            double d = (u[i] - mean);
+            double d2 = sqr(d);
+            var += d2;
+            skew += d * d2;
+            kurtosis += d2 * d2;
+        }
+        var /= u.size();
+        double dev = std::sqrt(var);
+        skew /= u.size() * dev * var;
+        kurtosis /= u.size() * var * var;
+        kurtosis -= 3;
+        double x_mean = ((double)d.a() + d.b()) / 2;
+        double x_var = (sqr((double)d.b() - d.a() + 1) - 1) / 12;
+        double x_skew = 0;
+        double x_kurtosis = -6. * (sqr((double)d.b() - d.a() + 1) + 1) /
+                            (5. * (sqr((double)d.b() - d.a() + 1) - 1));
+        assert(std::abs(mean - x_mean) / x_mean < 0.01);
+        assert(std::abs(var - x_var) / x_var < 0.01);
+        assert(std::abs(skew - x_skew) < 0.01);
+        assert(std::abs(kurtosis - x_kurtosis) / x_kurtosis < 0.01);
+    }
+    {
+        typedef std::uniform_int_distribution<> D;
+        typedef std::minstd_rand G;
+        G g;
+        D d;
+        const int N = 100000;
+        std::vector<D::result_type> u;
+        for (int i = 0; i < N; ++i)
+        {
+            D::result_type v = d(g);
+            assert(d.a() <= v && v <= d.b());
+            u.push_back(v);
+        }
+        double mean = std::accumulate(u.begin(), u.end(),
+                                              double(0)) / u.size();
+        double var = 0;
+        double skew = 0;
+        double kurtosis = 0;
+        for (int i = 0; i < u.size(); ++i)
+        {
+            double d = (u[i] - mean);
+            double d2 = sqr(d);
+            var += d2;
+            skew += d * d2;
+            kurtosis += d2 * d2;
+        }
+        var /= u.size();
+        double dev = std::sqrt(var);
+        skew /= u.size() * dev * var;
+        kurtosis /= u.size() * var * var;
+        kurtosis -= 3;
+        double x_mean = ((double)d.a() + d.b()) / 2;
+        double x_var = (sqr((double)d.b() - d.a() + 1) - 1) / 12;
+        double x_skew = 0;
+        double x_kurtosis = -6. * (sqr((double)d.b() - d.a() + 1) + 1) /
+                            (5. * (sqr((double)d.b() - d.a() + 1) - 1));
+        assert(std::abs(mean - x_mean) / x_mean < 0.01);
+        assert(std::abs(var - x_var) / x_var < 0.01);
+        assert(std::abs(skew - x_skew) < 0.01);
+        assert(std::abs(kurtosis - x_kurtosis) / x_kurtosis < 0.01);
+    }
+    {
+        typedef std::uniform_int_distribution<> D;
+        typedef std::mt19937 G;
+        G g;
+        D d;
+        const int N = 100000;
+        std::vector<D::result_type> u;
+        for (int i = 0; i < N; ++i)
+        {
+            D::result_type v = d(g);
+            assert(d.a() <= v && v <= d.b());
+            u.push_back(v);
+        }
+        double mean = std::accumulate(u.begin(), u.end(),
+                                              double(0)) / u.size();
+        double var = 0;
+        double skew = 0;
+        double kurtosis = 0;
+        for (int i = 0; i < u.size(); ++i)
+        {
+            double d = (u[i] - mean);
+            double d2 = sqr(d);
+            var += d2;
+            skew += d * d2;
+            kurtosis += d2 * d2;
+        }
+        var /= u.size();
+        double dev = std::sqrt(var);
+        skew /= u.size() * dev * var;
+        kurtosis /= u.size() * var * var;
+        kurtosis -= 3;
+        double x_mean = ((double)d.a() + d.b()) / 2;
+        double x_var = (sqr((double)d.b() - d.a() + 1) - 1) / 12;
+        double x_skew = 0;
+        double x_kurtosis = -6. * (sqr((double)d.b() - d.a() + 1) + 1) /
+                            (5. * (sqr((double)d.b() - d.a() + 1) - 1));
+        assert(std::abs(mean - x_mean) / x_mean < 0.01);
+        assert(std::abs(var - x_var) / x_var < 0.01);
+        assert(std::abs(skew - x_skew) < 0.01);
+        assert(std::abs(kurtosis - x_kurtosis) / x_kurtosis < 0.01);
+    }
+    {
+        typedef std::uniform_int_distribution<> D;
+        typedef std::mt19937_64 G;
+        G g;
+        D d;
+        const int N = 100000;
+        std::vector<D::result_type> u;
+        for (int i = 0; i < N; ++i)
+        {
+            D::result_type v = d(g);
+            assert(d.a() <= v && v <= d.b());
+            u.push_back(v);
+        }
+        double mean = std::accumulate(u.begin(), u.end(),
+                                              double(0)) / u.size();
+        double var = 0;
+        double skew = 0;
+        double kurtosis = 0;
+        for (int i = 0; i < u.size(); ++i)
+        {
+            double d = (u[i] - mean);
+            double d2 = sqr(d);
+            var += d2;
+            skew += d * d2;
+            kurtosis += d2 * d2;
+        }
+        var /= u.size();
+        double dev = std::sqrt(var);
+        skew /= u.size() * dev * var;
+        kurtosis /= u.size() * var * var;
+        kurtosis -= 3;
+        double x_mean = ((double)d.a() + d.b()) / 2;
+        double x_var = (sqr((double)d.b() - d.a() + 1) - 1) / 12;
+        double x_skew = 0;
+        double x_kurtosis = -6. * (sqr((double)d.b() - d.a() + 1) + 1) /
+                            (5. * (sqr((double)d.b() - d.a() + 1) - 1));
+        assert(std::abs(mean - x_mean) / x_mean < 0.01);
+        assert(std::abs(var - x_var) / x_var < 0.01);
+        assert(std::abs(skew - x_skew) < 0.01);
+        assert(std::abs(kurtosis - x_kurtosis) / x_kurtosis < 0.01);
+    }
+    {
+        typedef std::uniform_int_distribution<> D;
+        typedef std::ranlux24_base G;
+        G g;
+        D d;
+        const int N = 100000;
+        std::vector<D::result_type> u;
+        for (int i = 0; i < N; ++i)
+        {
+            D::result_type v = d(g);
+            assert(d.a() <= v && v <= d.b());
+            u.push_back(v);
+        }
+        double mean = std::accumulate(u.begin(), u.end(),
+                                              double(0)) / u.size();
+        double var = 0;
+        double skew = 0;
+        double kurtosis = 0;
+        for (int i = 0; i < u.size(); ++i)
+        {
+            double d = (u[i] - mean);
+            double d2 = sqr(d);
+            var += d2;
+            skew += d * d2;
+            kurtosis += d2 * d2;
+        }
+        var /= u.size();
+        double dev = std::sqrt(var);
+        skew /= u.size() * dev * var;
+        kurtosis /= u.size() * var * var;
+        kurtosis -= 3;
+        double x_mean = ((double)d.a() + d.b()) / 2;
+        double x_var = (sqr((double)d.b() - d.a() + 1) - 1) / 12;
+        double x_skew = 0;
+        double x_kurtosis = -6. * (sqr((double)d.b() - d.a() + 1) + 1) /
+                            (5. * (sqr((double)d.b() - d.a() + 1) - 1));
+        assert(std::abs(mean - x_mean) / x_mean < 0.01);
+        assert(std::abs(var - x_var) / x_var < 0.01);
+        assert(std::abs(skew - x_skew) < 0.01);
+        assert(std::abs(kurtosis - x_kurtosis) / x_kurtosis < 0.01);
+    }
+    {
+        typedef std::uniform_int_distribution<> D;
+        typedef std::ranlux48_base G;
+        G g;
+        D d;
+        const int N = 100000;
+        std::vector<D::result_type> u;
+        for (int i = 0; i < N; ++i)
+        {
+            D::result_type v = d(g);
+            assert(d.a() <= v && v <= d.b());
+            u.push_back(v);
+        }
+        double mean = std::accumulate(u.begin(), u.end(),
+                                              double(0)) / u.size();
+        double var = 0;
+        double skew = 0;
+        double kurtosis = 0;
+        for (int i = 0; i < u.size(); ++i)
+        {
+            double d = (u[i] - mean);
+            double d2 = sqr(d);
+            var += d2;
+            skew += d * d2;
+            kurtosis += d2 * d2;
+        }
+        var /= u.size();
+        double dev = std::sqrt(var);
+        skew /= u.size() * dev * var;
+        kurtosis /= u.size() * var * var;
+        kurtosis -= 3;
+        double x_mean = ((double)d.a() + d.b()) / 2;
+        double x_var = (sqr((double)d.b() - d.a() + 1) - 1) / 12;
+        double x_skew = 0;
+        double x_kurtosis = -6. * (sqr((double)d.b() - d.a() + 1) + 1) /
+                            (5. * (sqr((double)d.b() - d.a() + 1) - 1));
+        assert(std::abs(mean - x_mean) / x_mean < 0.01);
+        assert(std::abs(var - x_var) / x_var < 0.01);
+        assert(std::abs(skew - x_skew) < 0.01);
+        assert(std::abs(kurtosis - x_kurtosis) / x_kurtosis < 0.01);
+    }
+    {
+        typedef std::uniform_int_distribution<> D;
+        typedef std::ranlux24 G;
+        G g;
+        D d;
+        const int N = 100000;
+        std::vector<D::result_type> u;
+        for (int i = 0; i < N; ++i)
+        {
+            D::result_type v = d(g);
+            assert(d.a() <= v && v <= d.b());
+            u.push_back(v);
+        }
+        double mean = std::accumulate(u.begin(), u.end(),
+                                              double(0)) / u.size();
+        double var = 0;
+        double skew = 0;
+        double kurtosis = 0;
+        for (int i = 0; i < u.size(); ++i)
+        {
+            double d = (u[i] - mean);
+            double d2 = sqr(d);
+            var += d2;
+            skew += d * d2;
+            kurtosis += d2 * d2;
+        }
+        var /= u.size();
+        double dev = std::sqrt(var);
+        skew /= u.size() * dev * var;
+        kurtosis /= u.size() * var * var;
+        kurtosis -= 3;
+        double x_mean = ((double)d.a() + d.b()) / 2;
+        double x_var = (sqr((double)d.b() - d.a() + 1) - 1) / 12;
+        double x_skew = 0;
+        double x_kurtosis = -6. * (sqr((double)d.b() - d.a() + 1) + 1) /
+                            (5. * (sqr((double)d.b() - d.a() + 1) - 1));
+        assert(std::abs(mean - x_mean) / x_mean < 0.01);
+        assert(std::abs(var - x_var) / x_var < 0.01);
+        assert(std::abs(skew - x_skew) < 0.01);
+        assert(std::abs(kurtosis - x_kurtosis) / x_kurtosis < 0.01);
+    }
+    {
+        typedef std::uniform_int_distribution<> D;
+        typedef std::ranlux48 G;
+        G g;
+        D d;
+        const int N = 100000;
+        std::vector<D::result_type> u;
+        for (int i = 0; i < N; ++i)
+        {
+            D::result_type v = d(g);
+            assert(d.a() <= v && v <= d.b());
+            u.push_back(v);
+        }
+        double mean = std::accumulate(u.begin(), u.end(),
+                                              double(0)) / u.size();
+        double var = 0;
+        double skew = 0;
+        double kurtosis = 0;
+        for (int i = 0; i < u.size(); ++i)
+        {
+            double d = (u[i] - mean);
+            double d2 = sqr(d);
+            var += d2;
+            skew += d * d2;
+            kurtosis += d2 * d2;
+        }
+        var /= u.size();
+        double dev = std::sqrt(var);
+        skew /= u.size() * dev * var;
+        kurtosis /= u.size() * var * var;
+        kurtosis -= 3;
+        double x_mean = ((double)d.a() + d.b()) / 2;
+        double x_var = (sqr((double)d.b() - d.a() + 1) - 1) / 12;
+        double x_skew = 0;
+        double x_kurtosis = -6. * (sqr((double)d.b() - d.a() + 1) + 1) /
+                            (5. * (sqr((double)d.b() - d.a() + 1) - 1));
+        assert(std::abs(mean - x_mean) / x_mean < 0.01);
+        assert(std::abs(var - x_var) / x_var < 0.01);
+        assert(std::abs(skew - x_skew) < 0.01);
+        assert(std::abs(kurtosis - x_kurtosis) / x_kurtosis < 0.01);
+    }
+    {
+        typedef std::uniform_int_distribution<> D;
+        typedef std::knuth_b G;
+        G g;
+        D d;
+        const int N = 100000;
+        std::vector<D::result_type> u;
+        for (int i = 0; i < N; ++i)
+        {
+            D::result_type v = d(g);
+            assert(d.a() <= v && v <= d.b());
+            u.push_back(v);
+        }
+        double mean = std::accumulate(u.begin(), u.end(),
+                                              double(0)) / u.size();
+        double var = 0;
+        double skew = 0;
+        double kurtosis = 0;
+        for (int i = 0; i < u.size(); ++i)
+        {
+            double d = (u[i] - mean);
+            double d2 = sqr(d);
+            var += d2;
+            skew += d * d2;
+            kurtosis += d2 * d2;
+        }
+        var /= u.size();
+        double dev = std::sqrt(var);
+        skew /= u.size() * dev * var;
+        kurtosis /= u.size() * var * var;
+        kurtosis -= 3;
+        double x_mean = ((double)d.a() + d.b()) / 2;
+        double x_var = (sqr((double)d.b() - d.a() + 1) - 1) / 12;
+        double x_skew = 0;
+        double x_kurtosis = -6. * (sqr((double)d.b() - d.a() + 1) + 1) /
+                            (5. * (sqr((double)d.b() - d.a() + 1) - 1));
+        assert(std::abs(mean - x_mean) / x_mean < 0.01);
+        assert(std::abs(var - x_var) / x_var < 0.01);
+        assert(std::abs(skew - x_skew) < 0.01);
+        assert(std::abs(kurtosis - x_kurtosis) / x_kurtosis < 0.01);
+    }
+    {
+        typedef std::uniform_int_distribution<> D;
+        typedef std::minstd_rand0 G;
+        G g;
         D d(-6, 106);
         for (int i = 0; i < 10000; ++i)
         {
@@ -30,4 +409,45 @@
             assert(-6 <= u && u <= 106);
         }
     }
+    {
+        typedef std::uniform_int_distribution<> D;
+        typedef std::minstd_rand G;
+        G g;
+        D d(5, 100);
+        const int N = 100000;
+        std::vector<D::result_type> u;
+        for (int i = 0; i < N; ++i)
+        {
+            D::result_type v = d(g);
+            assert(d.a() <= v && v <= d.b());
+            u.push_back(v);
+        }
+        double mean = std::accumulate(u.begin(), u.end(),
+                                              double(0)) / u.size();
+        double var = 0;
+        double skew = 0;
+        double kurtosis = 0;
+        for (int i = 0; i < u.size(); ++i)
+        {
+            double d = (u[i] - mean);
+            double d2 = sqr(d);
+            var += d2;
+            skew += d * d2;
+            kurtosis += d2 * d2;
+        }
+        var /= u.size();
+        double dev = std::sqrt(var);
+        skew /= u.size() * dev * var;
+        kurtosis /= u.size() * var * var;
+        kurtosis -= 3;
+        double x_mean = ((double)d.a() + d.b()) / 2;
+        double x_var = (sqr((double)d.b() - d.a() + 1) - 1) / 12;
+        double x_skew = 0;
+        double x_kurtosis = -6. * (sqr((double)d.b() - d.a() + 1) + 1) /
+                            (5. * (sqr((double)d.b() - d.a() + 1) - 1));
+        assert(std::abs(mean - x_mean) / x_mean < 0.01);
+        assert(std::abs(var - x_var) / x_var < 0.01);
+        assert(std::abs(skew - x_skew) < 0.01);
+        assert(std::abs(kurtosis - x_kurtosis) / x_kurtosis < 0.01);
+    }
 }
diff --git a/test/numerics/rand/rand.dis/rand.dist.uni/rand.dist.uni.int/eval_param.pass.cpp b/test/numerics/rand/rand.dis/rand.dist.uni/rand.dist.uni.int/eval_param.pass.cpp
index 62507d7..6bbd4fa 100644
--- a/test/numerics/rand/rand.dis/rand.dist.uni/rand.dist.uni.int/eval_param.pass.cpp
+++ b/test/numerics/rand/rand.dis/rand.dist.uni/rand.dist.uni.int/eval_param.pass.cpp
@@ -16,20 +16,60 @@
 
 #include <random>
 #include <cassert>
+#include <vector>
+#include <numeric>
+
+template <class T>
+inline
+T
+sqr(T x)
+{
+    return x * x;
+}
 
 int main()
 {
     {
         typedef std::uniform_int_distribution<> D;
+        typedef std::minstd_rand G;
         typedef D::param_type P;
-        typedef std::minstd_rand0 G;
         G g;
-        D d(-6, 106);
+        D d(5, 100);
         P p(-10, 20);
-        for (int i = 0; i < 10000; ++i)
+        const int N = 100000;
+        std::vector<D::result_type> u;
+        for (int i = 0; i < N; ++i)
         {
-            int u = d(g, p);
-            assert(-10 <= u && u <= 20);
+            D::result_type v = d(g, p);
+            assert(p.a() <= v && v <= p.b());
+            u.push_back(v);
         }
+        double mean = std::accumulate(u.begin(), u.end(),
+                                              double(0)) / u.size();
+        double var = 0;
+        double skew = 0;
+        double kurtosis = 0;
+        for (int i = 0; i < u.size(); ++i)
+        {
+            double d = (u[i] - mean);
+            double d2 = sqr(d);
+            var += d2;
+            skew += d * d2;
+            kurtosis += d2 * d2;
+        }
+        var /= u.size();
+        double dev = std::sqrt(var);
+        skew /= u.size() * dev * var;
+        kurtosis /= u.size() * var * var;
+        kurtosis -= 3;
+        double x_mean = ((double)p.a() + p.b()) / 2;
+        double x_var = (sqr((double)p.b() - p.a() + 1) - 1) / 12;
+        double x_skew = 0;
+        double x_kurtosis = -6. * (sqr((double)p.b() - p.a() + 1) + 1) /
+                            (5. * (sqr((double)p.b() - p.a() + 1) - 1));
+        assert(std::abs(mean - x_mean) / x_mean < 0.01);
+        assert(std::abs(var - x_var) / x_var < 0.01);
+        assert(std::abs(skew - x_skew) < 0.01);
+        assert(std::abs(kurtosis - x_kurtosis) / x_kurtosis < 0.01);
     }
 }
diff --git a/test/numerics/rand/rand.dis/rand.dist.uni/rand.dist.uni.real/eval.pass.cpp b/test/numerics/rand/rand.dis/rand.dist.uni/rand.dist.uni.real/eval.pass.cpp
index 09ac9d4..b6de395 100644
--- a/test/numerics/rand/rand.dis/rand.dist.uni/rand.dist.uni.real/eval.pass.cpp
+++ b/test/numerics/rand/rand.dis/rand.dist.uni/rand.dist.uni.real/eval.pass.cpp
@@ -16,7 +16,16 @@
 
 #include <random>
 #include <cassert>
-#include <iostream>
+#include <vector>
+#include <numeric>
+
+template <class T>
+inline
+T
+sqr(T x)
+{
+    return x * x;
+}
 
 int main()
 {
@@ -24,11 +33,440 @@
         typedef std::uniform_real_distribution<> D;
         typedef std::minstd_rand0 G;
         G g;
-        D d(-6.5, 106.75);
-        for (int i = 0; i < 10000; ++i)
+        D d;
+        const int N = 100000;
+        std::vector<D::result_type> u;
+        for (int i = 0; i < N; ++i)
         {
-            D::result_type u = d(g);
-            assert(d.min() <= u && u < d.max());
+            D::result_type v = d(g);
+            assert(d.a() <= v && v < d.b());
+            u.push_back(v);
         }
+        D::result_type mean = std::accumulate(u.begin(), u.end(),
+                                              D::result_type(0)) / u.size();
+        D::result_type var = 0;
+        D::result_type skew = 0;
+        D::result_type kurtosis = 0;
+        for (int i = 0; i < u.size(); ++i)
+        {
+            D::result_type d = (u[i] - mean);
+            D::result_type d2 = sqr(d);
+            var += d2;
+            skew += d * d2;
+            kurtosis += d2 * d2;
+        }
+        var /= u.size();
+        D::result_type dev = std::sqrt(var);
+        skew /= u.size() * dev * var;
+        kurtosis /= u.size() * var * var;
+        kurtosis -= 3;
+        D::result_type x_mean = (d.a() + d.b()) / 2;
+        D::result_type x_var = sqr(d.b() - d.a()) / 12;
+        D::result_type x_skew = 0;
+        D::result_type x_kurtosis = -6./5;
+        assert(std::abs(mean - x_mean) / x_mean < 0.01);
+        assert(std::abs(var - x_var) / x_var < 0.01);
+        assert(std::abs(skew - x_skew) < 0.01);
+        assert(std::abs(kurtosis - x_kurtosis) / x_kurtosis < 0.01);
+    }
+    {
+        typedef std::uniform_real_distribution<> D;
+        typedef std::minstd_rand G;
+        G g;
+        D d;
+        const int N = 100000;
+        std::vector<D::result_type> u;
+        for (int i = 0; i < N; ++i)
+        {
+            D::result_type v = d(g);
+            assert(d.a() <= v && v < d.b());
+            u.push_back(v);
+        }
+        D::result_type mean = std::accumulate(u.begin(), u.end(),
+                                              D::result_type(0)) / u.size();
+        D::result_type var = 0;
+        D::result_type skew = 0;
+        D::result_type kurtosis = 0;
+        for (int i = 0; i < u.size(); ++i)
+        {
+            D::result_type d = (u[i] - mean);
+            D::result_type d2 = sqr(d);
+            var += d2;
+            skew += d * d2;
+            kurtosis += d2 * d2;
+        }
+        var /= u.size();
+        D::result_type dev = std::sqrt(var);
+        skew /= u.size() * dev * var;
+        kurtosis /= u.size() * var * var;
+        kurtosis -= 3;
+        D::result_type x_mean = (d.a() + d.b()) / 2;
+        D::result_type x_var = sqr(d.b() - d.a()) / 12;
+        D::result_type x_skew = 0;
+        D::result_type x_kurtosis = -6./5;
+        assert(std::abs(mean - x_mean) / x_mean < 0.01);
+        assert(std::abs(var - x_var) / x_var < 0.01);
+        assert(std::abs(skew - x_skew) < 0.01);
+        assert(std::abs(kurtosis - x_kurtosis) / x_kurtosis < 0.01);
+    }
+    {
+        typedef std::uniform_real_distribution<> D;
+        typedef std::mt19937 G;
+        G g;
+        D d;
+        const int N = 100000;
+        std::vector<D::result_type> u;
+        for (int i = 0; i < N; ++i)
+        {
+            D::result_type v = d(g);
+            assert(d.a() <= v && v < d.b());
+            u.push_back(v);
+        }
+        D::result_type mean = std::accumulate(u.begin(), u.end(),
+                                              D::result_type(0)) / u.size();
+        D::result_type var = 0;
+        D::result_type skew = 0;
+        D::result_type kurtosis = 0;
+        for (int i = 0; i < u.size(); ++i)
+        {
+            D::result_type d = (u[i] - mean);
+            D::result_type d2 = sqr(d);
+            var += d2;
+            skew += d * d2;
+            kurtosis += d2 * d2;
+        }
+        var /= u.size();
+        D::result_type dev = std::sqrt(var);
+        skew /= u.size() * dev * var;
+        kurtosis /= u.size() * var * var;
+        kurtosis -= 3;
+        D::result_type x_mean = (d.a() + d.b()) / 2;
+        D::result_type x_var = sqr(d.b() - d.a()) / 12;
+        D::result_type x_skew = 0;
+        D::result_type x_kurtosis = -6./5;
+        assert(std::abs(mean - x_mean) / x_mean < 0.01);
+        assert(std::abs(var - x_var) / x_var < 0.01);
+        assert(std::abs(skew - x_skew) < 0.01);
+        assert(std::abs(kurtosis - x_kurtosis) / x_kurtosis < 0.01);
+    }
+    {
+        typedef std::uniform_real_distribution<> D;
+        typedef std::mt19937_64 G;
+        G g;
+        D d;
+        const int N = 100000;
+        std::vector<D::result_type> u;
+        for (int i = 0; i < N; ++i)
+        {
+            D::result_type v = d(g);
+            assert(d.a() <= v && v < d.b());
+            u.push_back(v);
+        }
+        D::result_type mean = std::accumulate(u.begin(), u.end(),
+                                              D::result_type(0)) / u.size();
+        D::result_type var = 0;
+        D::result_type skew = 0;
+        D::result_type kurtosis = 0;
+        for (int i = 0; i < u.size(); ++i)
+        {
+            D::result_type d = (u[i] - mean);
+            D::result_type d2 = sqr(d);
+            var += d2;
+            skew += d * d2;
+            kurtosis += d2 * d2;
+        }
+        var /= u.size();
+        D::result_type dev = std::sqrt(var);
+        skew /= u.size() * dev * var;
+        kurtosis /= u.size() * var * var;
+        kurtosis -= 3;
+        D::result_type x_mean = (d.a() + d.b()) / 2;
+        D::result_type x_var = sqr(d.b() - d.a()) / 12;
+        D::result_type x_skew = 0;
+        D::result_type x_kurtosis = -6./5;
+        assert(std::abs(mean - x_mean) / x_mean < 0.01);
+        assert(std::abs(var - x_var) / x_var < 0.01);
+        assert(std::abs(skew - x_skew) < 0.01);
+        assert(std::abs(kurtosis - x_kurtosis) / x_kurtosis < 0.01);
+    }
+    {
+        typedef std::uniform_real_distribution<> D;
+        typedef std::ranlux24_base G;
+        G g;
+        D d;
+        const int N = 100000;
+        std::vector<D::result_type> u;
+        for (int i = 0; i < N; ++i)
+        {
+            D::result_type v = d(g);
+            assert(d.a() <= v && v < d.b());
+            u.push_back(v);
+        }
+        D::result_type mean = std::accumulate(u.begin(), u.end(),
+                                              D::result_type(0)) / u.size();
+        D::result_type var = 0;
+        D::result_type skew = 0;
+        D::result_type kurtosis = 0;
+        for (int i = 0; i < u.size(); ++i)
+        {
+            D::result_type d = (u[i] - mean);
+            D::result_type d2 = sqr(d);
+            var += d2;
+            skew += d * d2;
+            kurtosis += d2 * d2;
+        }
+        var /= u.size();
+        D::result_type dev = std::sqrt(var);
+        skew /= u.size() * dev * var;
+        kurtosis /= u.size() * var * var;
+        kurtosis -= 3;
+        D::result_type x_mean = (d.a() + d.b()) / 2;
+        D::result_type x_var = sqr(d.b() - d.a()) / 12;
+        D::result_type x_skew = 0;
+        D::result_type x_kurtosis = -6./5;
+        assert(std::abs(mean - x_mean) / x_mean < 0.01);
+        assert(std::abs(var - x_var) / x_var < 0.01);
+        assert(std::abs(skew - x_skew) < 0.02);
+        assert(std::abs(kurtosis - x_kurtosis) / x_kurtosis < 0.01);
+    }
+    {
+        typedef std::uniform_real_distribution<> D;
+        typedef std::ranlux48_base G;
+        G g;
+        D d;
+        const int N = 100000;
+        std::vector<D::result_type> u;
+        for (int i = 0; i < N; ++i)
+        {
+            D::result_type v = d(g);
+            assert(d.a() <= v && v < d.b());
+            u.push_back(v);
+        }
+        D::result_type mean = std::accumulate(u.begin(), u.end(),
+                                              D::result_type(0)) / u.size();
+        D::result_type var = 0;
+        D::result_type skew = 0;
+        D::result_type kurtosis = 0;
+        for (int i = 0; i < u.size(); ++i)
+        {
+            D::result_type d = (u[i] - mean);
+            D::result_type d2 = sqr(d);
+            var += d2;
+            skew += d * d2;
+            kurtosis += d2 * d2;
+        }
+        var /= u.size();
+        D::result_type dev = std::sqrt(var);
+        skew /= u.size() * dev * var;
+        kurtosis /= u.size() * var * var;
+        kurtosis -= 3;
+        D::result_type x_mean = (d.a() + d.b()) / 2;
+        D::result_type x_var = sqr(d.b() - d.a()) / 12;
+        D::result_type x_skew = 0;
+        D::result_type x_kurtosis = -6./5;
+        assert(std::abs(mean - x_mean) / x_mean < 0.01);
+        assert(std::abs(var - x_var) / x_var < 0.01);
+        assert(std::abs(skew - x_skew) < 0.01);
+        assert(std::abs(kurtosis - x_kurtosis) / x_kurtosis < 0.01);
+    }
+    {
+        typedef std::uniform_real_distribution<> D;
+        typedef std::ranlux24 G;
+        G g;
+        D d;
+        const int N = 100000;
+        std::vector<D::result_type> u;
+        for (int i = 0; i < N; ++i)
+        {
+            D::result_type v = d(g);
+            assert(d.a() <= v && v < d.b());
+            u.push_back(v);
+        }
+        D::result_type mean = std::accumulate(u.begin(), u.end(),
+                                              D::result_type(0)) / u.size();
+        D::result_type var = 0;
+        D::result_type skew = 0;
+        D::result_type kurtosis = 0;
+        for (int i = 0; i < u.size(); ++i)
+        {
+            D::result_type d = (u[i] - mean);
+            D::result_type d2 = sqr(d);
+            var += d2;
+            skew += d * d2;
+            kurtosis += d2 * d2;
+        }
+        var /= u.size();
+        D::result_type dev = std::sqrt(var);
+        skew /= u.size() * dev * var;
+        kurtosis /= u.size() * var * var;
+        kurtosis -= 3;
+        D::result_type x_mean = (d.a() + d.b()) / 2;
+        D::result_type x_var = sqr(d.b() - d.a()) / 12;
+        D::result_type x_skew = 0;
+        D::result_type x_kurtosis = -6./5;
+        assert(std::abs(mean - x_mean) / x_mean < 0.01);
+        assert(std::abs(var - x_var) / x_var < 0.01);
+        assert(std::abs(skew - x_skew) < 0.01);
+        assert(std::abs(kurtosis - x_kurtosis) / x_kurtosis < 0.01);
+    }
+    {
+        typedef std::uniform_real_distribution<> D;
+        typedef std::ranlux48 G;
+        G g;
+        D d;
+        const int N = 100000;
+        std::vector<D::result_type> u;
+        for (int i = 0; i < N; ++i)
+        {
+            D::result_type v = d(g);
+            assert(d.a() <= v && v < d.b());
+            u.push_back(v);
+        }
+        D::result_type mean = std::accumulate(u.begin(), u.end(),
+                                              D::result_type(0)) / u.size();
+        D::result_type var = 0;
+        D::result_type skew = 0;
+        D::result_type kurtosis = 0;
+        for (int i = 0; i < u.size(); ++i)
+        {
+            D::result_type d = (u[i] - mean);
+            D::result_type d2 = sqr(d);
+            var += d2;
+            skew += d * d2;
+            kurtosis += d2 * d2;
+        }
+        var /= u.size();
+        D::result_type dev = std::sqrt(var);
+        skew /= u.size() * dev * var;
+        kurtosis /= u.size() * var * var;
+        kurtosis -= 3;
+        D::result_type x_mean = (d.a() + d.b()) / 2;
+        D::result_type x_var = sqr(d.b() - d.a()) / 12;
+        D::result_type x_skew = 0;
+        D::result_type x_kurtosis = -6./5;
+        assert(std::abs(mean - x_mean) / x_mean < 0.01);
+        assert(std::abs(var - x_var) / x_var < 0.01);
+        assert(std::abs(skew - x_skew) < 0.01);
+        assert(std::abs(kurtosis - x_kurtosis) / x_kurtosis < 0.01);
+    }
+    {
+        typedef std::uniform_real_distribution<> D;
+        typedef std::knuth_b G;
+        G g;
+        D d;
+        const int N = 100000;
+        std::vector<D::result_type> u;
+        for (int i = 0; i < N; ++i)
+        {
+            D::result_type v = d(g);
+            assert(d.a() <= v && v < d.b());
+            u.push_back(v);
+        }
+        D::result_type mean = std::accumulate(u.begin(), u.end(),
+                                              D::result_type(0)) / u.size();
+        D::result_type var = 0;
+        D::result_type skew = 0;
+        D::result_type kurtosis = 0;
+        for (int i = 0; i < u.size(); ++i)
+        {
+            D::result_type d = (u[i] - mean);
+            D::result_type d2 = sqr(d);
+            var += d2;
+            skew += d * d2;
+            kurtosis += d2 * d2;
+        }
+        var /= u.size();
+        D::result_type dev = std::sqrt(var);
+        skew /= u.size() * dev * var;
+        kurtosis /= u.size() * var * var;
+        kurtosis -= 3;
+        D::result_type x_mean = (d.a() + d.b()) / 2;
+        D::result_type x_var = sqr(d.b() - d.a()) / 12;
+        D::result_type x_skew = 0;
+        D::result_type x_kurtosis = -6./5;
+        assert(std::abs(mean - x_mean) / x_mean < 0.01);
+        assert(std::abs(var - x_var) / x_var < 0.01);
+        assert(std::abs(skew - x_skew) < 0.01);
+        assert(std::abs(kurtosis - x_kurtosis) / x_kurtosis < 0.01);
+    }
+    {
+        typedef std::uniform_real_distribution<> D;
+        typedef std::minstd_rand G;
+        G g;
+        D d(-1, 1);
+        const int N = 100000;
+        std::vector<D::result_type> u;
+        for (int i = 0; i < N; ++i)
+        {
+            D::result_type v = d(g);
+            assert(d.a() <= v && v < d.b());
+            u.push_back(v);
+        }
+        D::result_type mean = std::accumulate(u.begin(), u.end(),
+                                              D::result_type(0)) / u.size();
+        D::result_type var = 0;
+        D::result_type skew = 0;
+        D::result_type kurtosis = 0;
+        for (int i = 0; i < u.size(); ++i)
+        {
+            D::result_type d = (u[i] - mean);
+            D::result_type d2 = sqr(d);
+            var += d2;
+            skew += d * d2;
+            kurtosis += d2 * d2;
+        }
+        var /= u.size();
+        D::result_type dev = std::sqrt(var);
+        skew /= u.size() * dev * var;
+        kurtosis /= u.size() * var * var;
+        kurtosis -= 3;
+        D::result_type x_mean = (d.a() + d.b()) / 2;
+        D::result_type x_var = sqr(d.b() - d.a()) / 12;
+        D::result_type x_skew = 0;
+        D::result_type x_kurtosis = -6./5;
+        assert(std::abs(mean - x_mean) < 0.01);
+        assert(std::abs(var - x_var) / x_var < 0.01);
+        assert(std::abs(skew - x_skew) < 0.01);
+        assert(std::abs(kurtosis - x_kurtosis) / x_kurtosis < 0.01);
+    }
+    {
+        typedef std::uniform_real_distribution<> D;
+        typedef std::minstd_rand G;
+        G g;
+        D d(5.5, 25);
+        const int N = 100000;
+        std::vector<D::result_type> u;
+        for (int i = 0; i < N; ++i)
+        {
+            D::result_type v = d(g);
+            assert(d.a() <= v && v < d.b());
+            u.push_back(v);
+        }
+        D::result_type mean = std::accumulate(u.begin(), u.end(),
+                                              D::result_type(0)) / u.size();
+        D::result_type var = 0;
+        D::result_type skew = 0;
+        D::result_type kurtosis = 0;
+        for (int i = 0; i < u.size(); ++i)
+        {
+            D::result_type d = (u[i] - mean);
+            D::result_type d2 = sqr(d);
+            var += d2;
+            skew += d * d2;
+            kurtosis += d2 * d2;
+        }
+        var /= u.size();
+        D::result_type dev = std::sqrt(var);
+        skew /= u.size() * dev * var;
+        kurtosis /= u.size() * var * var;
+        kurtosis -= 3;
+        D::result_type x_mean = (d.a() + d.b()) / 2;
+        D::result_type x_var = sqr(d.b() - d.a()) / 12;
+        D::result_type x_skew = 0;
+        D::result_type x_kurtosis = -6./5;
+        assert(std::abs(mean - x_mean) / x_mean < 0.01);
+        assert(std::abs(var - x_var) / x_var < 0.01);
+        assert(std::abs(skew - x_skew) < 0.01);
+        assert(std::abs(kurtosis - x_kurtosis) / x_kurtosis < 0.01);
     }
 }
diff --git a/test/numerics/rand/rand.dis/rand.dist.uni/rand.dist.uni.real/eval_param.pass.cpp b/test/numerics/rand/rand.dis/rand.dist.uni/rand.dist.uni.real/eval_param.pass.cpp
index 2b33ea4..79e2924 100644
--- a/test/numerics/rand/rand.dis/rand.dist.uni/rand.dist.uni.real/eval_param.pass.cpp
+++ b/test/numerics/rand/rand.dis/rand.dist.uni/rand.dist.uni.real/eval_param.pass.cpp
@@ -16,20 +16,59 @@
 
 #include <random>
 #include <cassert>
+#include <vector>
+#include <numeric>
+
+template <class T>
+inline
+T
+sqr(T x)
+{
+    return x * x;
+}
 
 int main()
 {
     {
         typedef std::uniform_real_distribution<> D;
+        typedef std::minstd_rand G;
         typedef D::param_type P;
-        typedef std::minstd_rand0 G;
         G g;
-        D d(-6, 106);
+        D d(5.5, 25);
         P p(-10, 20);
-        for (int i = 0; i < 10000; ++i)
+        const int N = 100000;
+        std::vector<D::result_type> u;
+        for (int i = 0; i < N; ++i)
         {
-            int u = d(g, p);
-            assert(p.a() <= u && u < p.b());
+            D::result_type v = d(g, p);
+            assert(p.a() <= v && v < p.b());
+            u.push_back(v);
         }
+        D::result_type mean = std::accumulate(u.begin(), u.end(),
+                                              D::result_type(0)) / u.size();
+        D::result_type var = 0;
+        D::result_type skew = 0;
+        D::result_type kurtosis = 0;
+        for (int i = 0; i < u.size(); ++i)
+        {
+            D::result_type d = (u[i] - mean);
+            D::result_type d2 = sqr(d);
+            var += d2;
+            skew += d * d2;
+            kurtosis += d2 * d2;
+        }
+        var /= u.size();
+        D::result_type dev = std::sqrt(var);
+        skew /= u.size() * dev * var;
+        kurtosis /= u.size() * var * var;
+        kurtosis -= 3;
+        D::result_type x_mean = (p.a() + p.b()) / 2;
+        D::result_type x_var = sqr(p.b() - p.a()) / 12;
+        D::result_type x_skew = 0;
+        D::result_type x_kurtosis = -6./5;
+        assert(std::abs(mean - x_mean) / x_mean < 0.01);
+        assert(std::abs(var - x_var) / x_var < 0.01);
+        assert(std::abs(skew - x_skew) < 0.01);
+        assert(std::abs(kurtosis - x_kurtosis) / x_kurtosis < 0.01);
     }
 }