From 1d8f611eb8ea1b1b821e9885fc1e214a30910670 Mon Sep 17 00:00:00 2001 From: lisyarus Date: Wed, 23 Sep 2020 11:12:23 +0300 Subject: [PATCH] Huge pcg/random rewrite: get rid of , use self-written distributions instead --- .../pcg/random/{point_box.hpp => box.hpp} | 15 +- libs/pcg/include/psemek/pcg/random/normal.hpp | 68 ++++++++ .../include/psemek/pcg/random/uniform_int.hpp | 93 +++++++++++ .../psemek/pcg/random/uniform_real.hpp | 74 +++++++++ .../psemek/pcg/random/uniform_sphere.hpp | 145 ++++++++++++++++++ .../psemek/pcg/random/vector_sphere.hpp | 101 ------------ libs/pcg/include/psemek/pcg/white.hpp | 7 +- 7 files changed, 388 insertions(+), 115 deletions(-) rename libs/pcg/include/psemek/pcg/random/{point_box.hpp => box.hpp} (62%) create mode 100644 libs/pcg/include/psemek/pcg/random/normal.hpp create mode 100644 libs/pcg/include/psemek/pcg/random/uniform_int.hpp create mode 100644 libs/pcg/include/psemek/pcg/random/uniform_real.hpp create mode 100644 libs/pcg/include/psemek/pcg/random/uniform_sphere.hpp delete mode 100644 libs/pcg/include/psemek/pcg/random/vector_sphere.hpp diff --git a/libs/pcg/include/psemek/pcg/random/point_box.hpp b/libs/pcg/include/psemek/pcg/random/box.hpp similarity index 62% rename from libs/pcg/include/psemek/pcg/random/point_box.hpp rename to libs/pcg/include/psemek/pcg/random/box.hpp index af90bd53..fbfd73c7 100644 --- a/libs/pcg/include/psemek/pcg/random/point_box.hpp +++ b/libs/pcg/include/psemek/pcg/random/box.hpp @@ -2,8 +2,7 @@ #include #include - -#include +#include namespace psemek::pcg { @@ -11,23 +10,17 @@ namespace psemek::pcg template struct box_point_distribution { - std::uniform_real_distribution d[N]; + uniform_real_distribution d[N]; box_point_distribution(geom::box const & b) { for (std::size_t i = 0; i < N; ++i) { - d[i] = std::uniform_real_distribution{b[i].min, b[i].max}; + d[i] = uniform_real_distribution{b[i]}; } } - box_point_distribution() - { - for (std::size_t i = 0; i < N; ++i) - { - d[i] = std::uniform_real_distribution{T{0}, T{1}}; - } - } + box_point_distribution() = default; template auto operator()(RNG && rng) diff --git a/libs/pcg/include/psemek/pcg/random/normal.hpp b/libs/pcg/include/psemek/pcg/random/normal.hpp new file mode 100644 index 00000000..585ebaee --- /dev/null +++ b/libs/pcg/include/psemek/pcg/random/normal.hpp @@ -0,0 +1,68 @@ +#pragma once + +#include +#include + +#include +#include + +namespace psemek::pcg +{ + + template + struct normal_distribution + { + static_assert(std::is_floating_point_v); + + 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 + T operator()(RNG && rng); + + private: + T mean_; + T stddev_; + + std::optional cached_; + }; + + template + template + T normal_distribution::operator()(RNG && rng) + { + // Box-Muller algorithm + + if (cached_) + { + T result = *cached_; + cached_ = std::nullopt; + return result; + } + + uniform_real_distribution 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; + } + +} diff --git a/libs/pcg/include/psemek/pcg/random/uniform_int.hpp b/libs/pcg/include/psemek/pcg/random/uniform_int.hpp new file mode 100644 index 00000000..b7c95393 --- /dev/null +++ b/libs/pcg/include/psemek/pcg/random/uniform_int.hpp @@ -0,0 +1,93 @@ +#pragma once + +#include + +#include +#include + +namespace psemek::pcg +{ + + template + struct uniform_int_distribution + { + static_assert(std::is_integral_v); + + using result_type = T; + + uniform_int_distribution() + : range_{std::numeric_limits::min(), std::numeric_limits::max()} + {} + + uniform_int_distribution(T min, T max) + : range_{min, max} + {} + + uniform_int_distribution(geom::interval range) + : range_(range) + {} + + uniform_int_distribution(uniform_int_distribution const &) = default; + + template + T operator()(RNG && rng); + + private: + geom::interval range_; + + template + static T generate(RNG & rng, geom::interval const & range); + }; + + template + template + T uniform_int_distribution::operator()(RNG && rng) + { + return generate(rng, range_); + } + + template + template + T uniform_int_distribution::generate(RNG & rng, geom::interval const & range) + { + using ctype = std::common_type_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(range.max) - static_cast(range.min); + + ctype result; + + if (rng_range == gen_range) + { + result = static_cast(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(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(rng()) - rng_min); + } + while (result > gen_range || result < tmp); + } + + return result + range.min; + } + +} diff --git a/libs/pcg/include/psemek/pcg/random/uniform_real.hpp b/libs/pcg/include/psemek/pcg/random/uniform_real.hpp new file mode 100644 index 00000000..ecb9157d --- /dev/null +++ b/libs/pcg/include/psemek/pcg/random/uniform_real.hpp @@ -0,0 +1,74 @@ +#pragma once + +#include + +#include +#include + +namespace psemek::pcg +{ + + template + struct uniform_real_distribution + { + static_assert(std::is_floating_point_v); + + 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 range) + : range_(range) + {} + + uniform_real_distribution(uniform_real_distribution const &) = default; + + template + T operator()(RNG && rng); + + private: + geom::interval range_; + + template + T generate(RNG & rng); + }; + + template + template + T uniform_real_distribution::operator()(RNG && rng) + { + return generate(rng) * range_.length() + range_.min; + } + + template + template + T uniform_real_distribution::generate(RNG & rng) + { + // see https://en.cppreference.com/w/cpp/numeric/random/generate_canonical + + static constexpr std::size_t digits = std::numeric_limits::digits; + static constexpr T const limit = std::pow(2, digits); + + T const R = static_cast(rng.max()) - static_cast(rng.min()) + 1; + + T mult = T(1); + + T sum = T(0); + + while (mult < limit) + { + T inc = static_cast(rng() - rng.min()); + sum += inc * mult; + mult *= R; + } + + return sum / mult; + } + +} diff --git a/libs/pcg/include/psemek/pcg/random/uniform_sphere.hpp b/libs/pcg/include/psemek/pcg/random/uniform_sphere.hpp new file mode 100644 index 00000000..663cb5c0 --- /dev/null +++ b/libs/pcg/include/psemek/pcg/random/uniform_sphere.hpp @@ -0,0 +1,145 @@ +#pragma once + +#include +#include +#include +#include +#include + +namespace psemek::pcg +{ + + template + struct uniform_sphere_vector_distribution + { + using result_type = geom::vector; + + uniform_sphere_vector_distribution(T r = T{1}) + : r_{r} + {} + + uniform_sphere_vector_distribution(uniform_sphere_vector_distribution const &) = default; + + template + 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 d_; + T r_; + }; + + template + struct uniform_sphere_vector_distribution + { + using result_type = geom::vector; + + uniform_sphere_vector_distribution(T r = T{1}) + : r_{r} + {} + + uniform_sphere_vector_distribution(uniform_sphere_vector_distribution const &) = default; + + template + geom::vector operator()(RNG && rng) + { + uniform_int_distribution d{0, 1}; + if (d(rng) == 0) + return {r_}; + return {-r_}; + } + + private: + T r_; + }; + + template + struct uniform_sphere_vector_distribution + { + using result_type = geom::vector; + + uniform_sphere_vector_distribution(T r = T{1}) + : r_{r} + {} + + uniform_sphere_vector_distribution(uniform_sphere_vector_distribution const &) = default; + + template + result_type operator()(RNG && rng) + { + uniform_real_distribution d{0, 2 * geom::pi}; + T a = d(rng); + return r_ * result_type{std::cos(a), std::sin(a)}; + } + + private: + T r_; + }; + + template + struct uniform_sphere_vector_distribution + { + using result_type = geom::vector; + + uniform_sphere_vector_distribution(T r = T{1}) + : r_{r} + {} + + template + auto operator()(RNG && rng) + { + uniform_real_distribution 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 + struct uniform_sphere_point_distribution + { + using result_type = geom::point; + + uniform_sphere_point_distribution(result_type origin, T r = T{1}) + : origin_{origin} + , vector_d_{r} + {} + + template + auto operator()(RNG && rng) + { + return origin_ + vector_d_(rng); + } + + private: + result_type origin_; + uniform_sphere_vector_distribution vector_d_; + }; + + +} diff --git a/libs/pcg/include/psemek/pcg/random/vector_sphere.hpp b/libs/pcg/include/psemek/pcg/random/vector_sphere.hpp deleted file mode 100644 index 927a68ed..00000000 --- a/libs/pcg/include/psemek/pcg/random/vector_sphere.hpp +++ /dev/null @@ -1,101 +0,0 @@ -#pragma once - -#include -#include - -#include - -namespace psemek::pcg -{ - - template - struct uniform_sphere_vector_distribution - { - std::normal_distribution d; - T r; - - uniform_sphere_vector_distribution(T r = T{1}) - : r{r} - {} - - template - auto operator()(RNG && rng) - { - geom::vector 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 - struct uniform_sphere_vector_distribution - { - std::uniform_int_distribution d{0, 1}; - T r; - - uniform_sphere_vector_distribution(T r = T{1}) - : r{r} - {} - - template - geom::vector operator()(RNG && rng) - { - if (d(rng) == 0) - return {r}; - return {-r}; - } - }; - - template - struct uniform_sphere_vector_distribution - { - std::uniform_real_distribution d{0, 2 * geom::pi}; - T r; - - uniform_sphere_vector_distribution(T r = T{1}) - : r{r} - {} - - template - geom::vector operator()(RNG && rng) - { - T a = d(rng); - return {r * std::cos(a), r * std::sin(a)}; - } - }; - - template - struct uniform_sphere_vector_distribution - { - std::uniform_real_distribution d{-1, 1}; - T r; - - uniform_sphere_vector_distribution(T r = T{1}) - : r{r} - {} - - template - auto operator()(RNG && rng) - { - geom::vector 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; - } - } - }; - -} diff --git a/libs/pcg/include/psemek/pcg/white.hpp b/libs/pcg/include/psemek/pcg/white.hpp index bc12aec1..fd660896 100644 --- a/libs/pcg/include/psemek/pcg/white.hpp +++ b/libs/pcg/include/psemek/pcg/white.hpp @@ -3,7 +3,8 @@ #include #include -#include +#include +#include namespace psemek::pcg { @@ -11,9 +12,9 @@ namespace psemek::pcg template gfx::basic_pixmap white(std::size_t width, std::size_t height, RNG && rng, geom::interval const & range) { - using dist = std::conditional_t, std::uniform_real_distribution, std::uniform_int_distribution>; + using dist = std::conditional_t, uniform_real_distribution, uniform_int_distribution>; - dist d{range.min, range.max}; + dist d{range}; gfx::basic_pixmap result(width, height);