#pragma once

#include <psemek/random/uniform_real.hpp>
#include <psemek/math/constants.hpp>

#include <cmath>
#include <optional>

namespace psemek::random
{

	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);

		T mean() const { return mean_; }

		T stddev() const { return stddev_; }

	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 u1 = 0.f;
		while (u1 == 0.f)
			u1 = d(rng);
		T const u2 = d(rng);

		T const r = std::sqrt(- 2 * std::log(u1));
		T const th = 2 * math::pi * u2;

		T const z1 = mean_ + stddev_ * r * std::cos(th);
		T const z2 = mean_ + stddev_ * r * std::sin(th);

		cached_ = z2;
		return z1;
	}

}
