#pragma once

#include <psemek/math/sign.hpp>
#include <psemek/math/point.hpp>
#include <psemek/math/orientation.hpp>
#include <psemek/math/robust.hpp>

#ifdef PSEMEK_ROBUST_PREDICATES
#include <boost/multiprecision/gmp.hpp>
#endif

#include <type_traits>

namespace psemek::math
{

	template <typename T, std::size_t N, typename ... Points>
	sign_t in_circle(fast_predicate_tag, point<T, N> const & p0, Points const & ... points)
	{
		auto proj = [](vector<T, N> const & v)
		{
			vector<T, N + 1> u;
			u[N] = T{0};
			for (std::size_t i = 0; i < N; ++i)
			{
				u[i] = v[i];
				u[N] += v[i] * v[i];
			}
			return u;
		};

		auto result = orientation(fast, proj(points - p0) ...);
		if constexpr ((N % 2) == 0)
			return inverse(result);
		else
			return result;
	}

	template <typename T>
	std::enable_if_t<std::is_floating_point_v<T>, sign_t>
		in_circle(robust_predicate_tag, point<T, 2> const & p0, point<T, 2> const & p1, point<T, 2> const & p2, point<T, 2> const & p3)
	{
#ifdef PSEMEK_ROBUST_PREDICATES
		constexpr T error = std::numeric_limits<T>::epsilon() * T(29) / T(2);

		T const m01 = (p0[0] - p3[0]) * (p1[1] - p3[1]);
		T const m02 = (p0[0] - p3[0]) * (p2[1] - p3[1]);
		T const m10 = (p1[0] - p3[0]) * (p0[1] - p3[1]);
		T const m12 = (p1[0] - p3[0]) * (p2[1] - p3[1]);
		T const m20 = (p2[0] - p3[0]) * (p0[1] - p3[1]);
		T const m21 = (p2[0] - p3[0]) * (p1[1] - p3[1]);

		T const d = T(0)
			+ m01 * p2[0] * p2[0] + m01 * p2[1] * p2[1] - m01 * p3[0] * p3[0] - m01 * p3[1] * p3[1]
			- m02 * p1[0] * p1[0] - m02 * p1[1] * p1[1] + m02 * p3[0] * p3[0] + m02 * p3[1] * p3[1]
			- m10 * p2[0] * p2[0] - m10 * p2[1] * p2[1] + m10 * p3[0] * p3[0] + m10 * p3[1] * p3[1]
			+ m12 * p0[0] * p0[0] + m12 * p0[1] * p0[1] - m12 * p3[0] * p3[0] - m12 * p3[1] * p3[1]
			+ m20 * p1[0] * p1[0] + m20 * p1[1] * p1[1] - m20 * p3[0] * p3[0] - m20 * p3[1] * p3[1]
			- m21 * p0[0] * p0[0] - m21 * p0[1] * p0[1] + m21 * p3[0] * p3[0] + m21 * p3[1] * p3[1]
			;

		T const t = T(0)
			+ std::abs(m01 * p2[0] * p2[0]) + std::abs(m01 * p2[1] * p2[1]) + std::abs(m01 * p3[0] * p3[0]) + std::abs(m01 * p3[1] * p3[1])
			+ std::abs(m02 * p1[0] * p1[0]) + std::abs(m02 * p1[1] * p1[1]) + std::abs(m02 * p3[0] * p3[0]) + std::abs(m02 * p3[1] * p3[1])
			+ std::abs(m10 * p2[0] * p2[0]) + std::abs(m10 * p2[1] * p2[1]) + std::abs(m10 * p3[0] * p3[0]) + std::abs(m10 * p3[1] * p3[1])
			+ std::abs(m12 * p0[0] * p0[0]) + std::abs(m12 * p0[1] * p0[1]) + std::abs(m12 * p3[0] * p3[0]) + std::abs(m12 * p3[1] * p3[1])
			+ std::abs(m20 * p1[0] * p1[0]) + std::abs(m20 * p1[1] * p1[1]) + std::abs(m20 * p3[0] * p3[0]) + std::abs(m20 * p3[1] * p3[1])
			+ std::abs(m21 * p0[0] * p0[0]) + std::abs(m21 * p0[1] * p0[1]) + std::abs(m21 * p3[0] * p3[0]) + std::abs(m21 * p3[1] * p3[1])
			;

		if (d > t * error)
			return sign_t::positive;
		else if (d < - t * error)
			return sign_t::negative;
		else
		{
			using exact_type = boost::multiprecision::mpq_rational;

			return in_circle(cast<exact_type>(p0), cast<exact_type>(p1), cast<exact_type>(p2), cast<exact_type>(p3));
		}
#else
		return in_circle(fast_predicate_tag{}, p0, p1, p2, p3);
#endif
	}

	template <typename ... Args>
	sign_t in_circle(Args const & ... args)
	{
		return in_circle(default_robust_tag, args...);
	}

}
