diff --git a/src/base/test/CMakeLists.txt b/src/base/test/CMakeLists.txt
index c57fc5568f..d135209a20 100644
--- a/src/base/test/CMakeLists.txt
+++ b/src/base/test/CMakeLists.txt
@@ -18,3 +18,11 @@ wl_test(test_base
     base_time_string
     base_utf8
 )
+
+wl_binary(rng_stats
+  SRCS
+    rng_stats.cc
+  DEPENDS
+    base
+    base_random
+)
diff --git a/src/base/test/rng_stats.cc b/src/base/test/rng_stats.cc
new file mode 100644
index 0000000000..a18a567b67
--- /dev/null
+++ b/src/base/test/rng_stats.cc
@@ -0,0 +1,78 @@
+/*
+ * Copyright (C) 2022 by the Widelands Development Team
+ *
+ * This program is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU General Public License
+ * as published by the Free Software Foundation; either version 2
+ * of the License, or (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program; if not, see <https://www.gnu.org/licenses/>.
+ *
+ */
+
+#include <cmath>
+#include <cstdint>
+#include <vector>
+
+#include "base/log.h"
+#include "base/random.h"
+
+/**
+ * Check the random number distribution in range [0, n).
+ * @param n Interval upper bound.
+ */
+static void test(const uint32_t n) {
+	constexpr uint32_t kSamples = 100000000;
+	RNG rng;
+	std::vector<uint32_t> results(n, 0);
+
+	uint32_t max = 1;
+	for (uint32_t i = kSamples; i > 0; --i) {
+		uint32_t& ref = results.at(rng.rand() % n);
+		++ref;
+		max = std::max(max, ref);
+	}
+
+	double mean = 0;
+	double mean_val = 0;
+	for (uint32_t i = 0; i < n; ++i) {
+		mean += results.at(i);
+		mean_val += i * results.at(i);
+	}
+	mean /= results.size();
+	mean_val /= kSamples;
+
+	double sigma = 0;
+	for (uint32_t r : results) {
+		double d = r;
+		d -= mean;
+		sigma += d * d;
+	}
+	sigma = sqrt(sigma / n);
+
+	log_info("****************************************************************************");
+	log_info("           N = %u", n);
+	log_info("    mean_val = %f", mean_val);
+	log_info("        mean = %f", mean);
+	log_info("       sigma = %f", sigma);
+
+	constexpr size_t kMaxColumnW = 50;
+	for (uint32_t i = 0; i < n; ++i) {
+		std::string vis(kMaxColumnW * results.at(i) / max, '#');
+		log_info("[%10u] = %10u %s", i, results.at(i), vis.c_str());
+	}
+}
+
+/** Main function. */
+int main(int, char**) {
+	set_testcase_logging_dir();
+	for (uint32_t n : {6, 100, 0x10000}) {
+		test(n);
+	}
+}