diff --git a/libcxx/include/CMakeLists.txt b/libcxx/include/CMakeLists.txt
--- a/libcxx/include/CMakeLists.txt
+++ b/libcxx/include/CMakeLists.txt
@@ -268,6 +268,7 @@
   __random/geometric_distribution.h
   __random/independent_bits_engine.h
   __random/is_seed_sequence.h
+  __random/is_valid_inttype.h
   __random/knuth_b.h
   __random/linear_congruential_engine.h
   __random/log2.h
diff --git a/libcxx/include/__random/binomial_distribution.h b/libcxx/include/__random/binomial_distribution.h
--- a/libcxx/include/__random/binomial_distribution.h
+++ b/libcxx/include/__random/binomial_distribution.h
@@ -10,6 +10,7 @@
 #define _LIBCPP___RANDOM_BINOMIAL_DISTRIBUTION_H
 
 #include <__config>
+#include <__random/is_valid_inttype.h>
 #include <__random/uniform_real_distribution.h>
 #include <cmath>
 #include <iosfwd>
@@ -26,6 +27,7 @@
 template<class _IntType = int>
 class _LIBCPP_TEMPLATE_VIS binomial_distribution
 {
+    static_assert(__random_is_valid_inttype<_IntType>::value, "IntType must be an integer type larger than char");
 public:
     // types
     typedef _IntType result_type;
diff --git a/libcxx/include/__random/discrete_distribution.h b/libcxx/include/__random/discrete_distribution.h
--- a/libcxx/include/__random/discrete_distribution.h
+++ b/libcxx/include/__random/discrete_distribution.h
@@ -11,6 +11,7 @@
 
 #include <__algorithm/upper_bound.h>
 #include <__config>
+#include <__random/is_valid_inttype.h>
 #include <__random/uniform_real_distribution.h>
 #include <cstddef>
 #include <iosfwd>
@@ -29,6 +30,7 @@
 template<class _IntType = int>
 class _LIBCPP_TEMPLATE_VIS discrete_distribution
 {
+    static_assert(__random_is_valid_inttype<_IntType>::value, "IntType must be an integer type larger than char");
 public:
     // types
     typedef _IntType result_type;
diff --git a/libcxx/include/__random/geometric_distribution.h b/libcxx/include/__random/geometric_distribution.h
--- a/libcxx/include/__random/geometric_distribution.h
+++ b/libcxx/include/__random/geometric_distribution.h
@@ -10,6 +10,7 @@
 #define _LIBCPP___RANDOM_GEOMETRIC_DISTRIBUTION_H
 
 #include <__config>
+#include <__random/is_valid_inttype.h>
 #include <__random/negative_binomial_distribution.h>
 #include <iosfwd>
 #include <limits>
@@ -26,6 +27,7 @@
 template<class _IntType = int>
 class _LIBCPP_TEMPLATE_VIS geometric_distribution
 {
+    static_assert(__random_is_valid_inttype<_IntType>::value, "IntType must be an integer type larger than char");
 public:
     // types
     typedef _IntType result_type;
diff --git a/libcxx/include/__random/is_valid_inttype.h b/libcxx/include/__random/is_valid_inttype.h
new file mode 100644
--- /dev/null
+++ b/libcxx/include/__random/is_valid_inttype.h
@@ -0,0 +1,44 @@
+//===----------------------------------------------------------------------===//
+//
+// Part of the LLVM Project, under the Apache License v2.0 with LLVM Exceptions.
+// See https://llvm.org/LICENSE.txt for license information.
+// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception
+//
+//===----------------------------------------------------------------------===//
+
+#ifndef _LIBCPP___RANDOM_IS_VALID_INTTYPE_H
+#define _LIBCPP___RANDOM_IS_VALID_INTTYPE_H
+
+#include <__config>
+#include <type_traits>
+
+#if !defined(_LIBCPP_HAS_NO_PRAGMA_SYSTEM_HEADER)
+#pragma GCC system_header
+#endif
+
+_LIBCPP_BEGIN_NAMESPACE_STD
+
+// [rand.req.genl]/1.5:
+// The effect of instantiating a template that has a template type parameter
+// named IntType is undefined unless the corresponding template argument is
+// cv-unqualified and is one of short, int, long, long long, unsigned short,
+// unsigned int, unsigned long, or unsigned long long.
+
+template<class> struct __random_is_valid_inttype : false_type {};
+template<> struct __random_is_valid_inttype<short> : true_type {};
+template<> struct __random_is_valid_inttype<int> : true_type {};
+template<> struct __random_is_valid_inttype<long> : true_type {};
+template<> struct __random_is_valid_inttype<long long> : true_type {};
+template<> struct __random_is_valid_inttype<unsigned short> : true_type {};
+template<> struct __random_is_valid_inttype<unsigned int> : true_type {};
+template<> struct __random_is_valid_inttype<unsigned long> : true_type {};
+template<> struct __random_is_valid_inttype<unsigned long long> : true_type {};
+
+#ifndef _LIBCPP_HAS_NO_INT128
+template<> struct __random_is_valid_inttype<__int128_t> : true_type {};
+template<> struct __random_is_valid_inttype<__uint128_t> : true_type {};
+#endif // _LIBCPP_HAS_NO_INT128
+
+_LIBCPP_END_NAMESPACE_STD
+
+#endif // _LIBCPP___RANDOM_IS_VALID_INTTYPE_H
diff --git a/libcxx/include/__random/negative_binomial_distribution.h b/libcxx/include/__random/negative_binomial_distribution.h
--- a/libcxx/include/__random/negative_binomial_distribution.h
+++ b/libcxx/include/__random/negative_binomial_distribution.h
@@ -12,6 +12,7 @@
 #include <__config>
 #include <__random/bernoulli_distribution.h>
 #include <__random/gamma_distribution.h>
+#include <__random/is_valid_inttype.h>
 #include <__random/poisson_distribution.h>
 #include <iosfwd>
 #include <limits>
@@ -28,6 +29,7 @@
 template<class _IntType = int>
 class _LIBCPP_TEMPLATE_VIS negative_binomial_distribution
 {
+    static_assert(__random_is_valid_inttype<_IntType>::value, "IntType must be an integer type larger than char");
 public:
     // types
     typedef _IntType result_type;
diff --git a/libcxx/include/__random/poisson_distribution.h b/libcxx/include/__random/poisson_distribution.h
--- a/libcxx/include/__random/poisson_distribution.h
+++ b/libcxx/include/__random/poisson_distribution.h
@@ -11,6 +11,7 @@
 
 #include <__config>
 #include <__random/exponential_distribution.h>
+#include <__random/is_valid_inttype.h>
 #include <__random/normal_distribution.h>
 #include <__random/uniform_real_distribution.h>
 #include <cmath>
@@ -29,6 +30,7 @@
 template<class _IntType = int>
 class _LIBCPP_TEMPLATE_VIS poisson_distribution
 {
+    static_assert(__random_is_valid_inttype<_IntType>::value, "IntType must be an integer type larger than char");
 public:
     // types
     typedef _IntType result_type;
diff --git a/libcxx/include/__random/uniform_int_distribution.h b/libcxx/include/__random/uniform_int_distribution.h
--- a/libcxx/include/__random/uniform_int_distribution.h
+++ b/libcxx/include/__random/uniform_int_distribution.h
@@ -11,6 +11,7 @@
 
 #include <__bits>
 #include <__config>
+#include <__random/is_valid_inttype.h>
 #include <__random/log2.h>
 #include <bit>
 #include <cstddef>
@@ -155,9 +156,10 @@
     return _Sp;
 }
 
-template<class _IntType = int> // __int128_t is also supported as an extension here
+template<class _IntType = int>
 class uniform_int_distribution
 {
+    static_assert(__random_is_valid_inttype<_IntType>::value, "IntType must be an integer type larger than char");
 public:
     // types
     typedef _IntType result_type;
diff --git a/libcxx/include/module.modulemap b/libcxx/include/module.modulemap
--- a/libcxx/include/module.modulemap
+++ b/libcxx/include/module.modulemap
@@ -715,6 +715,7 @@
       module geometric_distribution          { private header "__random/geometric_distribution.h" }
       module independent_bits_engine         { private header "__random/independent_bits_engine.h" }
       module is_seed_sequence                { private header "__random/is_seed_sequence.h" }
+      module is_valid_inttype                { private header "__random/is_valid_inttype.h" }
       module knuth_b                         { private header "__random/knuth_b.h" }
       module linear_congruential_engine      { private header "__random/linear_congruential_engine.h" }
       module log2                            { private header "__random/log2.h" }
diff --git a/libcxx/include/random b/libcxx/include/random
--- a/libcxx/include/random
+++ b/libcxx/include/random
@@ -1693,6 +1693,7 @@
 #include <__random/geometric_distribution.h>
 #include <__random/independent_bits_engine.h>
 #include <__random/is_seed_sequence.h>
+#include <__random/is_valid_inttype.h>
 #include <__random/knuth_b.h>
 #include <__random/linear_congruential_engine.h>
 #include <__random/log2.h>
diff --git a/libcxx/test/libcxx/diagnostics/detail.headers/random/is_valid_inttype.module.verify.cpp b/libcxx/test/libcxx/diagnostics/detail.headers/random/is_valid_inttype.module.verify.cpp
new file mode 100644
--- /dev/null
+++ b/libcxx/test/libcxx/diagnostics/detail.headers/random/is_valid_inttype.module.verify.cpp
@@ -0,0 +1,15 @@
+//===----------------------------------------------------------------------===//
+//
+// Part of the LLVM Project, under the Apache License v2.0 with LLVM Exceptions.
+// See https://llvm.org/LICENSE.txt for license information.
+// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception
+//
+//===----------------------------------------------------------------------===//
+
+// REQUIRES: modules-build
+
+// WARNING: This test was generated by 'generate_private_header_tests.py'
+// and should not be edited manually.
+
+// expected-error@*:* {{use of private header from outside its module: '__random/is_valid_inttype.h'}}
+#include <__random/is_valid_inttype.h>
diff --git a/libcxx/test/std/numerics/rand/rand.dis/rand.dist.uni/rand.dist.uni.int/eval.pass.cpp b/libcxx/test/std/numerics/rand/rand.dis/rand.dist.uni/rand.dist.uni.int/eval.pass.cpp
--- a/libcxx/test/std/numerics/rand/rand.dis/rand.dist.uni/rand.dist.uni.int/eval.pass.cpp
+++ b/libcxx/test/std/numerics/rand/rand.dis/rand.dist.uni/rand.dist.uni.int/eval.pass.cpp
@@ -17,443 +17,125 @@
 
 #include <random>
 #include <cassert>
-#include <vector>
-#include <numeric>
+#include <climits>
 #include <cstddef>
+#include <limits>
+#include <numeric>
+#include <vector>
 
 #include "test_macros.h"
 
 template <class T>
-inline
-T
-sqr(T x)
+T sqr(T x)
 {
     return x * x;
 }
 
-int main(int, char**)
+template <class ResultType, class EngineType>
+void test_statistics(ResultType a, ResultType b)
 {
-    {
-        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 (std::size_t i = 0; i < u.size(); ++i)
-        {
-            double dbl = (u[i] - mean);
-            double d2 = sqr(dbl);
-            var += d2;
-            skew += dbl * 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 (std::size_t i = 0; i < u.size(); ++i)
-        {
-            double dbl = (u[i] - mean);
-            double d2 = sqr(dbl);
-            var += d2;
-            skew += dbl * 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 (std::size_t i = 0; i < u.size(); ++i)
-        {
-            double dbl = (u[i] - mean);
-            double d2 = sqr(dbl);
-            var += d2;
-            skew += dbl * 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 (std::size_t i = 0; i < u.size(); ++i)
-        {
-            double dbl = (u[i] - mean);
-            double d2 = sqr(dbl);
-            var += d2;
-            skew += dbl * 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 (std::size_t i = 0; i < u.size(); ++i)
-        {
-            double dbl = (u[i] - mean);
-            double d2 = sqr(dbl);
-            var += d2;
-            skew += dbl * 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 (std::size_t i = 0; i < u.size(); ++i)
-        {
-            double dbl = (u[i] - mean);
-            double d2 = sqr(dbl);
-            var += d2;
-            skew += dbl * 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 (std::size_t i = 0; i < u.size(); ++i)
-        {
-            double dbl = (u[i] - mean);
-            double d2 = sqr(dbl);
-            var += d2;
-            skew += dbl * 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 (std::size_t i = 0; i < u.size(); ++i)
-        {
-            double dbl = (u[i] - mean);
-            double d2 = sqr(dbl);
-            var += d2;
-            skew += dbl * 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 (std::size_t i = 0; i < u.size(); ++i)
-        {
-            double dbl = (u[i] - mean);
-            double d2 = sqr(dbl);
-            var += d2;
-            skew += dbl * 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);
+    ASSERT_SAME_TYPE(typename std::uniform_int_distribution<ResultType>::result_type, ResultType);
+
+    EngineType g;
+    std::uniform_int_distribution<ResultType> dist(a, b);
+    assert(dist.a() == a);
+    assert(dist.b() == b);
+    std::vector<ResultType> u;
+    for (int i = 0; i < 10000; ++i) {
+        ResultType v = dist(g);
+        assert(a <= v && v <= b);
+        u.push_back(v);
     }
-    {
-        typedef std::uniform_int_distribution<> D;
-        typedef std::minstd_rand0 G;
-        G g;
-        D d(-6, 106);
-        for (int i = 0; i < 10000; ++i)
-        {
-            int u = d(g);
-            assert(-6 <= u && u <= 106);
-        }
+
+    // Quick check: The chance of getting *no* hits in any given tenth of the range
+    // is (0.9)^10000, or "ultra-astronomically low."
+    bool bottom_tenth = false;
+    bool top_tenth = false;
+    for (std::size_t i = 0; i < u.size(); ++i) {
+        bottom_tenth = bottom_tenth || (u[i] <= (a + (b / 10) - (a / 10)));
+        top_tenth = top_tenth || (u[i] >= (b - (b / 10) + (a / 10)));
     }
-    {
-        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 (std::size_t i = 0; i < u.size(); ++i)
-        {
-            double dbl = (u[i] - mean);
-            double d2 = sqr(dbl);
-            var += d2;
-            skew += dbl * 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);
+    assert(bottom_tenth);  // ...is populated
+    assert(top_tenth);  // ...is populated
+
+    // Now do some more involved statistical math.
+    double mean = std::accumulate(u.begin(), u.end(), 0.0) / u.size();
+    double var = 0;
+    double skew = 0;
+    double kurtosis = 0;
+    for (std::size_t i = 0; i < u.size(); ++i) {
+        double dbl = (u[i] - mean);
+        double d2 = dbl * dbl;
+        var += d2;
+        skew += dbl * d2;
+        kurtosis += d2 * d2;
     }
+    var /= u.size();
+    double dev = std::sqrt(var);
+    skew /= u.size() * dev * var;
+    kurtosis /= u.size() * var * var;
+
+    double expected_mean = double(a) + double(b)/2 - double(a)/2;
+    double expected_var = (sqr(double(b) - double(a) + 1) - 1) / 12;
+
+    double range = double(b) - double(a) + 1.0;
+    assert(range > range / 10);  // i.e., it's not infinity
+
+    assert(std::abs(mean - expected_mean) < range / 100);
+    assert(std::abs(var - expected_var) < expected_var / 50);
+    assert(-0.1 < skew && skew < 0.1);
+    assert(1.6 < kurtosis && kurtosis < 2.0);
+}
+
+template <class ResultType, class EngineType>
+void test_statistics()
+{
+    test_statistics<ResultType, EngineType>(0, std::numeric_limits<ResultType>::max());
+}
+
+int main(int, char**)
+{
+    test_statistics<int, std::minstd_rand0>();
+    test_statistics<int, std::minstd_rand>();
+    test_statistics<int, std::mt19937>();
+    test_statistics<int, std::mt19937_64>();
+    test_statistics<int, std::ranlux24_base>();
+    test_statistics<int, std::ranlux48_base>();
+    test_statistics<int, std::ranlux24>();
+    test_statistics<int, std::ranlux48>();
+    test_statistics<int, std::knuth_b>();
+    test_statistics<int, std::minstd_rand0>(-6, 106);
+    test_statistics<int, std::minstd_rand>(5, 100);
+
+    test_statistics<short, std::minstd_rand0>();
+    test_statistics<int, std::minstd_rand0>();
+    test_statistics<long, std::minstd_rand0>();
+    test_statistics<long long, std::minstd_rand0>();
+
+    test_statistics<unsigned short, std::minstd_rand0>();
+    test_statistics<unsigned int, std::minstd_rand0>();
+    test_statistics<unsigned long, std::minstd_rand0>();
+    test_statistics<unsigned long long, std::minstd_rand0>();
+
+    test_statistics<short, std::minstd_rand0>(SHRT_MIN, SHRT_MAX);
+
+    // http://eel.is/c++draft/rand.req#genl-1.5
+    // The effect of instantiating a template that has a parameter
+    // named IntType is undefined unless the corresponding template
+    // argument is cv-unqualified and is one of short, int, long,
+    // long long, unsigned short, unsigned int, unsigned long,
+    // or unsigned long long.
+    // (We support __int128 as an extension.)
+
+#ifndef _LIBCPP_HAS_NO_INT128
+    test_statistics<__int128_t, std::minstd_rand0>();
+    test_statistics<__uint128_t, std::minstd_rand0>();
+
+    test_statistics<__int128_t, std::minstd_rand0>(-100, 900);
+    test_statistics<__int128_t, std::minstd_rand0>(0, UINT64_MAX);
+    test_statistics<__int128_t, std::minstd_rand0>(std::numeric_limits<__int128_t>::min(), std::numeric_limits<__int128_t>::max());
+    test_statistics<__uint128_t, std::minstd_rand0>(0, UINT64_MAX);
+#endif
 
-  return 0;
+    return 0;
 }