[libc++] Explicitly reject `uniform_int_distribution<bool>` and `<char>`.

[• Quuxplusone]

Summary

`uniform_int_distribution<T>` is UB unless `T` is one of the non-character,
non-boolean integer types (`short` or larger). However, libc++ has never
enforced this. D114129 accidentally made `uniform_int_distribution<bool>`
into an error. Make it now *intentionally* an error; and likewise for the
character types and all user-defined class and enum types; but permit
`__[u]int128_t` to continue working.

Apply the same static_assert to all the integer distributions.

[• Quuxplusone]

libcxx/test/std/numerics/rand/rand.dis/rand.dist.uni/rand.dist.uni.int/eval.pass.cpp

69–76

Before you ask: yeah, these error bounds are crazy ad-hoc. But:

• the test isn't "flaky": in the sense that the engine g will always produce the same random sequence, so our code for dist should always cook that sequence into something that satisfies these arbitrary error bounds, and if it ever changes so that the error bounds fail, then we can just go edit this test.
• the old test was way more ad-hoc: like, what even is std::abs((mean - x_mean) / x_mean) < 0.01? Dividing by the mean of the distribution? We're just lucky they never tried a distribution with a mean of zero. 😀

[jloser]

libcxx/test/std/numerics/rand/rand.dis/rand.dist.uni/rand.dist.uni.int/eval.pass.cpp

69–76

One of those "oh boy" moments when looking at the existing trunk code (e.g. mean of zero comment).

[ldionne]

Let's just *not* support uniform_int_distribution<bool>, and make it into a static_assert.

Also, the test improvements should be made as a different patch. probably before this one cause it is needed to test with a wide variety of types like you do (which is a good improvement BTW).

[• Quuxplusone]

Let's just *not* support uniform_int_distribution<bool>, and make it into a static_assert.

Works for me, but (1) I'll feel better once we hear from @leonardchan, and (2) what do you want to do about char, signed char, char16_t, etc.? Also disallow them via static_assert, I'd guess?

[ldionne]

In D114920#3167060, @Quuxplusone wrote:

Let's just *not* support uniform_int_distribution<bool>, and make it into a static_assert.

Works for me, but (1) I'll feel better once we hear from @leonardchan, and (2) what do you want to do about char, signed char, char16_t, etc.? Also disallow them via static_assert, I'd guess?

Yeah, I'd do the same. If something is library UB and it's easy to diagnose, I would diagnose it.

[• Quuxplusone]

Change the intent of this PR from "re-enable" to "explicitly reject".

[• Quuxplusone]

libcxx/include/__random/is_valid_inttype.h

21–40

I'm not thrilled about this helper or its placement/naming, but it seemed like the best option.
Incidentally, this also gives people a sneaky libc++ internal detail that they can (unsupportedly) specialize in order to re-enable the removed functionality for specific user-defined types.
Should this helper be named __libcpp_*?

[ldionne]

Couple of comments.

Also, I'm going to try this patch out on a significant code base just to confirm it doesn't result in crazy breakage -- requesting changes until that's done.

Ref for myself: 86111407&86111417

libcxx/include/__random/is_valid_inttype.h

21–40

Yeah, let's use __libcpp_*. I don't like that name, but I think it's worth doing it if it reduces the likelihood of users messing around with our internal details.

37

This is where we could document that we support this as an extension, in a comment.

libcxx/test/std/numerics/rand/rand.dis/rand.dist.uni/rand.dist.uni.int/eval.pass.cpp

13

Any appetite for improving tests for other random distributions like this? It's *much* better than what we had before.

[ldionne]

Okay, so I have some initial results of my investigation. I found some uses of uniform_int_distribution with things like uint8_t, char, and unsigned char. There was definitely some breakage, however nothing that seemed entirely unmanageable.

I say we go for it, but we definitely need to add a release note. Also, we'll want to sync up with some folks that build their code bases from libc++ trunk (like I think @thakis @rnk @phosek) or they might encounter issues. They should proof their code bases before we land this (or if that can't be done in a reasonable time frame, then they can revert this commit internally).

One thing we could do to make things easier would be to have the static_asserts in a different commit, that way folks could just revert that one locally until they have solved their problems.

[• Quuxplusone]

@ldionne: Makes perfect sense to me. Would you care to suggest wording for a release note? Perhaps

(API Changes)
The integer distributions ``binomial_distribution``, ``discrete_distribution``,
``geometric_distribution``, ``negative_binomial_distribution``, ``poisson_distribution``,
and ``uniform_int_distribution`` now conform to the Standard by rejecting
template parameter types other than ``short``, ``int``, ``long``, ``long long``,
(as an extension) ``__int128_t``, and the unsigned versions thereof.
In particular, ``uniform_int_distribution<int8_t>`` is no longer supported.

Also, since "everything but the static_asserts" is really just "the test file," shall I go ahead and land the test-file change by itself, soonish?

[• Quuxplusone]

libcxx/include/__random/is_valid_inttype.h

21–40

Will rename to __libcpp_random_is_valid_inttype.

Also yuck, libcxx/benchmarks/vector_operations.bench.cpp uses uniform_int_distribution<char> in order to get random example data for its vectors. Maybe we should just support that. It would match libstdc++, too.

uniform_int_distribution<char/int8_t/uint8_t>: GCC accepts, MSVC rejects with static_assert.
uniform_int_distribution<bool>: GCC and MSVC both reject, with ugly template error message.

[ldionne]

In D114920#3177676, @Quuxplusone wrote:

@ldionne: Makes perfect sense to me. Would you care to suggest wording for a release note? Perhaps

(API Changes)
The integer distributions ``binomial_distribution``, ``discrete_distribution``,
``geometric_distribution``, ``negative_binomial_distribution``, ``poisson_distribution``,
and ``uniform_int_distribution`` now conform to the Standard by rejecting
template parameter types other than ``short``, ``int``, ``long``, ``long long``,
(as an extension) ``__int128_t``, and the unsigned versions thereof.
In particular, ``uniform_int_distribution<int8_t>`` is no longer supported.

Yes, that sounds perfect to me.

Also, since "everything but the static_asserts" is really just "the test file," shall I go ahead and land the test-file change by itself, soonish?

Yes! You can do this as a NFC commit and repurpose this patch to be only the static asserts and the release note, and fixing the benchmarks.

libcxx/include/__random/is_valid_inttype.h

21–40

So basically we'd only exclude bool? IMO we should strive for strict conformance, however we should also write a half-page paper to support it for char and other small types cause it kind of makes sense.

[• Quuxplusone]

Add a release note.
Fix the benchmarks.
Scope-creep to [libc++] Explicitly reject URNG types with signed result_types (fixes https://github.com/llvm/llvm-project/issues/48965 ) while I'm here.

[ldionne]

Can we add libc++ specific tests to check that we properly diagnose those? One big .verify.cpp test for all the distributions would be great. LGTM apart from that.

[• Quuxplusone]

Rebase. Split out "Explicitly reject URNG types with signed result_types" into its own PR.
If CI on this part is green, I'll probably just land it on Monday, since it's been approved since forever.

[• Quuxplusone]

Add missing pragmas.

libcxx/include/CMakeLists.txt

set(files
   __random/exponential_distribution.h
   __random/extreme_value_distribution.h
   __random/fisher_f_distribution.h
   __random/gamma_distribution.h
   __random/generate_canonical.h
   __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
   __random/lognormal_distribution.h
   __random/mersenne_twister_engine.h
   __random/negative_binomial_distribution.h
   __random/normal_distribution.h
   __random/piecewise_constant_distribution.h

libcxx/include/__random/binomial_distribution.h

 //===----------------------------------------------------------------------===//
 //
 // 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_BINOMIAL_DISTRIBUTION_H
 #define _LIBCPP___RANDOM_BINOMIAL_DISTRIBUTION_H
 
 #include <__config>
+#include <__random/is_valid_inttype.h>
 #include <__random/uniform_real_distribution.h>
 #include <cmath>
 #include <iosfwd>
 
 #if !defined(_LIBCPP_HAS_NO_PRAGMA_SYSTEM_HEADER)
 #pragma GCC system_header
 #endif
 
 _LIBCPP_PUSH_MACROS
 #include <__undef_macros>
 
 _LIBCPP_BEGIN_NAMESPACE_STD
 
 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;
 
     class _LIBCPP_TEMPLATE_VIS param_type
     {
         result_type __t_;
         double __p_;

libcxx/include/__random/discrete_distribution.h

 //===----------------------------------------------------------------------===//
 //
 // 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_DISCRETE_DISTRIBUTION_H
 #define _LIBCPP___RANDOM_DISCRETE_DISTRIBUTION_H
 
 #include <__algorithm/upper_bound.h>
 #include <__config>
+#include <__random/is_valid_inttype.h>
 #include <__random/uniform_real_distribution.h>
 #include <cstddef>
 #include <iosfwd>
 #include <numeric>
 #include <vector>
 
 #if !defined(_LIBCPP_HAS_NO_PRAGMA_SYSTEM_HEADER)
 #pragma GCC system_header
 #endif
 
 _LIBCPP_PUSH_MACROS
 #include <__undef_macros>
 
 _LIBCPP_BEGIN_NAMESPACE_STD
 
 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;
 
     class _LIBCPP_TEMPLATE_VIS param_type
     {
         vector<double> __p_;
     public:

libcxx/include/__random/geometric_distribution.h

 //===----------------------------------------------------------------------===//
 //
 // 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_GEOMETRIC_DISTRIBUTION_H
 #define _LIBCPP___RANDOM_GEOMETRIC_DISTRIBUTION_H
 
 #include <__config>
+#include <__random/is_valid_inttype.h>
 #include <__random/negative_binomial_distribution.h>
 #include <iosfwd>
 #include <limits>
 
 #if !defined(_LIBCPP_HAS_NO_PRAGMA_SYSTEM_HEADER)
 #pragma GCC system_header
 #endif
 
 _LIBCPP_PUSH_MACROS
 #include <__undef_macros>
 
 _LIBCPP_BEGIN_NAMESPACE_STD
 
 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;
 
     class _LIBCPP_TEMPLATE_VIS param_type
     {
         double __p_;
     public:

libcxx/include/__random/is_valid_inttype.h

+//===----------------------------------------------------------------------===//
+//
+// 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

libcxx/include/__random/negative_binomial_distribution.h

 //===----------------------------------------------------------------------===//
 //
 // 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_NEGATIVE_BINOMIAL_DISTRIBUTION_H
 #define _LIBCPP___RANDOM_NEGATIVE_BINOMIAL_DISTRIBUTION_H
 
 #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>
 
 #if !defined(_LIBCPP_HAS_NO_PRAGMA_SYSTEM_HEADER)
 #pragma GCC system_header
 #endif
 
 _LIBCPP_PUSH_MACROS
 #include <__undef_macros>
 
 _LIBCPP_BEGIN_NAMESPACE_STD
 
 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;
 
     class _LIBCPP_TEMPLATE_VIS param_type
     {
         result_type __k_;
         double __p_;

libcxx/include/__random/poisson_distribution.h

 //===----------------------------------------------------------------------===//
 //
 // 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_POISSON_DISTRIBUTION_H
 #define _LIBCPP___RANDOM_POISSON_DISTRIBUTION_H
 
 #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>
 #include <iosfwd>
 #include <limits>
 
 #if !defined(_LIBCPP_HAS_NO_PRAGMA_SYSTEM_HEADER)
 #pragma GCC system_header
 #endif
 
 _LIBCPP_PUSH_MACROS
 #include <__undef_macros>
 
 _LIBCPP_BEGIN_NAMESPACE_STD
 
 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;
 
     class _LIBCPP_TEMPLATE_VIS param_type
     {
         double __mean_;
         double __s_;

libcxx/include/__random/uniform_int_distribution.h

 //===----------------------------------------------------------------------===//
 //
 // 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_UNIFORM_INT_DISTRIBUTION_H
 #define _LIBCPP___RANDOM_UNIFORM_INT_DISTRIBUTION_H
 
 #include <__bits>
 #include <__config>
+#include <__random/is_valid_inttype.h>
 #include <__random/log2.h>
 #include <bit>
 #include <cstddef>
 #include <cstdint>
 #include <iosfwd>
 #include <limits>
 #include <type_traits>
 
for (size_t __k = __n0_; __k < __n_; ++__k)
             _Sp <<= __w0_ + 1;
         else
             _Sp = 0;
         _Sp += __u & __mask1_;
     }
     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;
 
     class param_type
     {
         result_type __a_;
         result_type __b_;

libcxx/include/module.modulemap

module __random {
       module exponential_distribution        { private header "__random/exponential_distribution.h" }
       module extreme_value_distribution      { private header "__random/extreme_value_distribution.h" }
       module fisher_f_distribution           { private header "__random/fisher_f_distribution.h" }
       module gamma_distribution              { private header "__random/gamma_distribution.h" }
       module generate_canonical              { private header "__random/generate_canonical.h" }
       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" }
       module lognormal_distribution          { private header "__random/lognormal_distribution.h" }
       module mersenne_twister_engine         { private header "__random/mersenne_twister_engine.h" }
       module negative_binomial_distribution  { private header "__random/negative_binomial_distribution.h" }
       module normal_distribution             { private header "__random/normal_distribution.h" }
       module piecewise_constant_distribution { private header "__random/piecewise_constant_distribution.h" }

libcxx/include/random

 #include <__random/exponential_distribution.h>
 #include <__random/extreme_value_distribution.h>
 #include <__random/fisher_f_distribution.h>
 #include <__random/gamma_distribution.h>
 #include <__random/generate_canonical.h>
 #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>
 #include <__random/lognormal_distribution.h>
 #include <__random/mersenne_twister_engine.h>
 #include <__random/negative_binomial_distribution.h>
 #include <__random/normal_distribution.h>
 #include <__random/piecewise_constant_distribution.h>

libcxx/test/libcxx/diagnostics/detail.headers/random/is_valid_inttype.module.verify.cpp

+//===----------------------------------------------------------------------===//
+//
+// 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>

libcxx/test/std/numerics/rand/rand.dis/rand.dist.uni/rand.dist.uni.int/eval.pass.cpp

 //===----------------------------------------------------------------------===//
 //
 // 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: long_tests
 
 // <random>
 
 // template<class _IntType = int>
 // class uniform_int_distribution
 
 // template<class _URNG> result_type operator()(_URNG& g);
 
 #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**)
-{
-    {
-        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);
-    }
-    {
-        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);
-        }
-    }
-    {
-        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);
-    }
+template <class ResultType, class EngineType>
+void test_statistics(ResultType a, ResultType b)
+    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);
+
+    // 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)));
+    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();
+    for (std::size_t i = 0; i < u.size(); ++i) {
+        double d2 = dbl * dbl;
+
+    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;
 }