Huge pcg/random rewrite: get rid of <random>, use self-written distributions instead

libs/pcg/include/psemek/pcg/random/box.hpp

@@ -2,8 +2,7 @@
 
 #include <psemek/geom/point.hpp>
 #include <psemek/geom/box.hpp>
-
+#include <psemek/pcg/random/uniform_real.hpp>
-#include <random>
 
 namespace psemek::pcg
 {
@@ -11,23 +10,17 @@ namespace psemek::pcg
 	template <typename T, std::size_t N>
 	struct box_point_distribution
 	{
-		std::uniform_real_distribution<T> d[N];
+		uniform_real_distribution<T> d[N];
 
 		box_point_distribution(geom::box<T, N> const & b)
 		{
 			for (std::size_t i = 0; i < N; ++i)
 			{
-				d[i] = std::uniform_real_distribution<T>{b[i].min, b[i].max};
+				d[i] = uniform_real_distribution<T>{b[i]};
 			}
 		}
 
-		box_point_distribution()
+		box_point_distribution() = default;
-		{
-			for (std::size_t i = 0; i < N; ++i)
-			{
-				d[i] = std::uniform_real_distribution<T>{T{0}, T{1}};
-			}
-		}
 
 		template <typename RNG>
 		auto operator()(RNG && rng)

libs/pcg/include/psemek/pcg/random/normal.hpp

@@ -0,0 +1,68 @@
+#pragma once
+
+#include <psemek/pcg/random/uniform_real.hpp>
+#include <psemek/geom/constants.hpp>
+
+#include <cmath>
+#include <optional>
+
+namespace psemek::pcg
+{
+
+	template <typename T>
+	struct normal_distribution
+	{
+		static_assert(std::is_floating_point_v<T>);
+
+		using result_type = T;
+
+		normal_distribution()
+			: mean_{T(0)}
+			, stddev_{T(1)}
+		{}
+
+		normal_distribution(T mean, T stddev)
+			: mean_{mean}
+			, stddev_{stddev}
+		{}
+
+		normal_distribution(normal_distribution const &) = default;
+
+		template <typename RNG>
+		T operator()(RNG && rng);
+
+	private:
+		T mean_;
+		T stddev_;
+
+		std::optional<T> cached_;
+	};
+
+	template <typename T>
+	template <typename RNG>
+	T normal_distribution<T>::operator()(RNG && rng)
+	{
+		// Box-Muller algorithm
+
+		if (cached_)
+		{
+			T result = *cached_;
+			cached_ = std::nullopt;
+			return result;
+		}
+
+		uniform_real_distribution<T> d;
+		T const u1 = d(rng);
+		T const u2 = d(rng);
+
+		T const r = std::sqrt(- 2 * std::log(u1));
+		T const th = 2 * geom::pi * u2;
+
+		T const z1 = r * std::cos(th);
+		T const z2 = r * std::sin(th);
+
+		cached_ = z2;
+		return z1;
+	}
+
+}

libs/pcg/include/psemek/pcg/random/uniform_int.hpp

@@ -0,0 +1,93 @@
+#pragma once
+
+#include <psemek/geom/interval.hpp>
+
+#include <limits>
+#include <type_traits>
+
+namespace psemek::pcg
+{
+
+	template <typename T>
+	struct uniform_int_distribution
+	{
+		static_assert(std::is_integral_v<T>);
+
+		using result_type = T;
+
+		uniform_int_distribution()
+			: range_{std::numeric_limits<T>::min(), std::numeric_limits<T>::max()}
+		{}
+
+		uniform_int_distribution(T min, T max)
+			: range_{min, max}
+		{}
+
+		uniform_int_distribution(geom::interval<T> range)
+			: range_(range)
+		{}
+
+		uniform_int_distribution(uniform_int_distribution const &) = default;
+
+		template <typename RNG>
+		T operator()(RNG && rng);
+
+	private:
+		geom::interval<T> range_;
+
+		template <typename RNG>
+		static T generate(RNG & rng, geom::interval<T> const & range);
+	};
+
+	template <typename T>
+	template <typename RNG>
+	T uniform_int_distribution<T>::operator()(RNG && rng)
+	{
+		return generate(rng, range_);
+	}
+
+	template <typename T>
+	template <typename RNG>
+	T uniform_int_distribution<T>::generate(RNG & rng, geom::interval<T> const & range)
+	{
+		using ctype = std::common_type_t<typename RNG::result_type, std::make_unsigned_t<T>>;
+
+		ctype const rng_min = rng.min();
+		ctype const rng_max = rng.max();
+		ctype const rng_range = rng_max - rng_min;
+		ctype const gen_range = static_cast<ctype>(range.max) - static_cast<ctype>(range.min);
+
+		ctype result;
+
+		if (rng_range == gen_range)
+		{
+			result = static_cast<ctype>(rng()) - rng_min;
+		}
+		if (rng_range > gen_range)
+		{
+			ctype const gen_range_full = gen_range + 1;
+			ctype const buckets = rng_range / gen_range_full;
+			ctype const max = gen_range * buckets;
+
+			do
+			{
+				result = static_cast<ctype>(rng()) - rng_min;
+			} while (result >= max);
+			result /= buckets;
+		}
+		else // if (rng_range < gen_range)
+		{
+			ctype tmp;
+			do
+			{
+				const ctype rng_range_full = rng_range + 1;
+				tmp = rng_range_full * generate(rng, {0, gen_range / rng_range_full});
+				result = tmp + (static_cast<ctype>(rng()) - rng_min);
+			}
+			while (result > gen_range || result < tmp);
+		}
+
+		return result + range.min;
+	}
+
+}

libs/pcg/include/psemek/pcg/random/uniform_real.hpp

@@ -0,0 +1,74 @@
+#pragma once
+
+#include <psemek/geom/interval.hpp>
+
+#include <limits>
+#include <cmath>
+
+namespace psemek::pcg
+{
+
+	template <typename T>
+	struct uniform_real_distribution
+	{
+		static_assert(std::is_floating_point_v<T>);
+
+		using result_type = T;
+
+		uniform_real_distribution()
+			: range_{T(0), T(1)}
+		{}
+
+		uniform_real_distribution(T min, T max)
+			: range_{min, max}
+		{}
+
+		uniform_real_distribution(geom::interval<T> range)
+			: range_(range)
+		{}
+
+		uniform_real_distribution(uniform_real_distribution const &) = default;
+
+		template <typename RNG>
+		T operator()(RNG && rng);
+
+	private:
+		geom::interval<T> range_;
+
+		template <typename RNG>
+		T generate(RNG & rng);
+	};
+
+	template <typename T>
+	template <typename RNG>
+	T uniform_real_distribution<T>::operator()(RNG && rng)
+	{
+		return generate(rng) * range_.length() + range_.min;
+	}
+
+	template <typename T>
+	template <typename RNG>
+	T uniform_real_distribution<T>::generate(RNG & rng)
+	{
+		// see https://en.cppreference.com/w/cpp/numeric/random/generate_canonical
+
+		static constexpr std::size_t digits = std::numeric_limits<T>::digits;
+		static constexpr T const limit = std::pow<T>(2, digits);
+
+		T const R = static_cast<T>(rng.max()) - static_cast<T>(rng.min()) + 1;
+
+		T mult = T(1);
+
+		T sum = T(0);
+
+		while (mult < limit)
+		{
+			T inc = static_cast<T>(rng() - rng.min());
+			sum += inc * mult;
+			mult *= R;
+		}
+
+		return sum / mult;
+	}
+
+}

libs/pcg/include/psemek/pcg/random/uniform_sphere.hpp

@@ -0,0 +1,145 @@
+#pragma once
+
+#include <psemek/geom/vector.hpp>
+#include <psemek/geom/point.hpp>
+#include <psemek/geom/constants.hpp>
+#include <psemek/pcg/random/normal.hpp>
+#include <psemek/pcg/random/uniform_int.hpp>
+
+namespace psemek::pcg
+{
+
+	template <typename T, std::size_t N>
+	struct uniform_sphere_vector_distribution
+	{
+		using result_type = geom::vector<T, N>;
+
+		uniform_sphere_vector_distribution(T r = T{1})
+			: r_{r}
+		{}
+
+		uniform_sphere_vector_distribution(uniform_sphere_vector_distribution const &) = default;
+
+		template <typename RNG>
+		auto operator()(RNG && rng)
+		{
+			result_type v;
+			while (true)
+			{
+				for (std::size_t i = 0; i < N; ++i)
+					v[i] = d_(rng);
+
+				T l = length(v);
+				if (l != T{})
+					return v / l * r_;
+			}
+		}
+
+	private:
+		normal_distribution<T> d_;
+		T r_;
+	};
+
+	template <typename T>
+	struct uniform_sphere_vector_distribution<T, 1>
+	{
+		using result_type = geom::vector<T, 1>;
+
+		uniform_sphere_vector_distribution(T r = T{1})
+			: r_{r}
+		{}
+
+		uniform_sphere_vector_distribution(uniform_sphere_vector_distribution const &) = default;
+
+		template <typename RNG>
+		geom::vector<T, 1> operator()(RNG && rng)
+		{
+			uniform_int_distribution<int> d{0, 1};
+			if (d(rng) == 0)
+				return {r_};
+			return {-r_};
+		}
+
+	private:
+		T r_;
+	};
+
+	template <typename T>
+	struct uniform_sphere_vector_distribution<T, 2>
+	{
+		using result_type = geom::vector<T, 2>;
+
+		uniform_sphere_vector_distribution(T r = T{1})
+			: r_{r}
+		{}
+
+		uniform_sphere_vector_distribution(uniform_sphere_vector_distribution const &) = default;
+
+		template <typename RNG>
+		result_type operator()(RNG && rng)
+		{
+			uniform_real_distribution<T> d{0, 2 * geom::pi};
+			T a = d(rng);
+			return r_ * result_type{std::cos(a), std::sin(a)};
+		}
+
+	private:
+		T r_;
+	};
+
+	template <typename T>
+	struct uniform_sphere_vector_distribution<T, 3>
+	{
+		using result_type = geom::vector<T, 3>;
+
+		uniform_sphere_vector_distribution(T r = T{1})
+			: r_{r}
+		{}
+
+		template <typename RNG>
+		auto operator()(RNG && rng)
+		{
+			uniform_real_distribution<T> d{-1, 1};
+
+			while (true)
+			{
+				T const x = d(rng);
+				T const y = d(rng);
+
+				T const s = x * x + y * y;
+
+				if (s > 1) continue;
+
+				T const m = 2 * std::sqrt(1 - s);
+
+				return r_ * result_type{x * m, y * m, 2 * s - 1};
+			}
+		}
+
+	private:
+		T r_;
+	};
+
+	template <typename T, std::size_t N>
+	struct uniform_sphere_point_distribution
+	{
+		using result_type = geom::point<T, N>;
+
+		uniform_sphere_point_distribution(result_type origin, T r = T{1})
+			: origin_{origin}
+			, vector_d_{r}
+		{}
+
+		template <typename RNG>
+		auto operator()(RNG && rng)
+		{
+			return origin_ + vector_d_(rng);
+		}
+
+	private:
+		result_type origin_;
+		uniform_sphere_vector_distribution<T, N> vector_d_;
+	};
+
+
+}

libs/pcg/include/psemek/pcg/random/vector_sphere.hpp

@@ -1,101 +0,0 @@
-#pragma once
-
-#include <psemek/geom/vector.hpp>
-#include <psemek/geom/constants.hpp>
-
-#include <random>
-
-namespace psemek::pcg
-{
-
-	template <typename T, std::size_t N>
-	struct uniform_sphere_vector_distribution
-	{
-		std::normal_distribution<T> d;
-		T r;
-
-		uniform_sphere_vector_distribution(T r = T{1})
-			: r{r}
-		{}
-
-		template <typename RNG>
-		auto operator()(RNG && rng)
-		{
-			geom::vector<T, N> v;
-			while (true)
-			{
-				for (std::size_t i = 0; i < N; ++i)
-					v[i] = d(rng);
-
-				T l = length(v);
-				if (l != T{})
-					return v / l * r;
-			}
-		}
-	};
-
-	template <typename T>
-	struct uniform_sphere_vector_distribution<T, 1>
-	{
-		std::uniform_int_distribution<int> d{0, 1};
-		T r;
-
-		uniform_sphere_vector_distribution(T r = T{1})
-			: r{r}
-		{}
-
-		template <typename RNG>
-		geom::vector<T, 1> operator()(RNG && rng)
-		{
-			if (d(rng) == 0)
-				return {r};
-			return {-r};
-		}
-	};
-
-	template <typename T>
-	struct uniform_sphere_vector_distribution<T, 2>
-	{
-		std::uniform_real_distribution<T> d{0, 2 * geom::pi};
-		T r;
-
-		uniform_sphere_vector_distribution(T r = T{1})
-			: r{r}
-		{}
-
-		template <typename RNG>
-		geom::vector<T, 2> operator()(RNG && rng)
-		{
-			T a = d(rng);
-			return {r * std::cos(a), r * std::sin(a)};
-		}
-	};
-
-	template <typename T>
-	struct uniform_sphere_vector_distribution<T, 3>
-	{
-		std::uniform_real_distribution<T> d{-1, 1};
-		T r;
-
-		uniform_sphere_vector_distribution(T r = T{1})
-			: r{r}
-		{}
-
-		template <typename RNG>
-		auto operator()(RNG && rng)
-		{
-			geom::vector<T, 3> v;
-			while (true)
-			{
-				for (std::size_t i = 0; i < 3; ++i)
-					v[i] = d(rng);
-
-				T l = length(v);
-
-				if (l <= T{1})
-					return v / l * r;
-			}
-		}
-	};
-
-}

libs/pcg/include/psemek/pcg/white.hpp

@@ -3,7 +3,8 @@
 #include <psemek/geom/interval.hpp>
 #include <psemek/gfx/pixmap.hpp>
 
-#include <random>
+#include <psemek/pcg/random/uniform_int.hpp>
+#include <psemek/pcg/random/uniform_real.hpp>
 
 namespace psemek::pcg
 {
@@ -11,9 +12,9 @@ namespace psemek::pcg
 	template <typename T = float, typename RNG>
 	gfx::basic_pixmap<T> white(std::size_t width, std::size_t height, RNG && rng, geom::interval<T> const & range)
 	{
-		using dist = std::conditional_t<std::is_floating_point_v<T>, std::uniform_real_distribution<T>, std::uniform_int_distribution<T>>;
+		using dist = std::conditional_t<std::is_floating_point_v<T>, uniform_real_distribution<T>, uniform_int_distribution<T>>;
 
-		dist d{range.min, range.max};
+		dist d{range};
 
 		gfx::basic_pixmap<T> result(width, height);