Compiler Explorer

Source code

#include <algorithm>
#include <array>
#include <bit>
#include <cfloat>
#include <climits>
#include <cmath>
#include <cstdint>
#include <limits>
#include <stdexcept>
#include <version>

namespace in_range_ext {
	// integral excluding bool and character types
	template<class T>
	concept integer = std::is_integral_v<T> &&
		!std::is_same_v<std::remove_cv_t<T>, bool> &&
		!std::is_same_v<std::remove_cv_t<T>, char> &&
		!std::is_same_v<std::remove_cv_t<T>, char8_t> &&
		!std::is_same_v<std::remove_cv_t<T>, char16_t> &&
		!std::is_same_v<std::remove_cv_t<T>, char32_t> &&
		!std::is_same_v<std::remove_cv_t<T>, wchar_t>;

// bool in_range<integer>(floating_point f)
	// returns true if f is in the range of values that can be represented by the integer type
	template<integer I, std::floating_point F> constexpr bool in_range(F f);

// integer clamp(floating_point f, integer min, integer max)
	// returns f clamped to the integer range [min, max]
	// currently throws std::domain_error if isnan(f)
	template<integer I, std::floating_point F> constexpr I clamp(F f, I min, I max);

namespace detail {
		namespace constexpr_cmath {
			// need a few constexpr <cmath> functions, but these aren't widely available yet.
			// the stand-ins are slow, but only used to compute compile-time constants.
#if __cpp_lib_constexpr_cmath >= 202202L
			using std::copysign;
			using std::signbit;
			using std::fpclassify;
			using std::ilogb;
			using std::scalbn;
#else
			constexpr bool internal_test = false;

template<std::floating_point F>
			constexpr F copysign(F f, F s)
			{
				if (!internal_test && !std::is_constant_evaluated())
					return std::copysign(f, s);

#if defined __clang__ || defined __GNUC__
				return __builtin_copysign(f, s);
#elif defined _MSC_VER
				using U = std::conditional_t<std::is_same_v<std::remove_cv_t<F>, float>, uint32_t, uint64_t>;
				constexpr U signbit = std::numeric_limits::max() ^ (std::numeric_limits::max() >> 1);
				U sgn = std::bit_cast(s) & signbit;
				U mag = std::bit_cast(f) & ~signbit;
				return std::bit_cast<F>(sgn | mag);
#else
				#error I don't know what to do
#endif
			}

template<std::floating_point F>
			constexpr bool signbit(F f)
			{
				if (!internal_test && !std::is_constant_evaluated())
					return std::signbit(f);

return copysign(F(1.0), f) < 0;
			}

template<std::floating_point F>
			constexpr int fpclassify(F f)
			{
				if (!internal_test && !std::is_constant_evaluated())
					return std::fpclassify(f);

if (f == 0)
					return FP_ZERO;
				else if (-std::numeric_limits<F>::min() < f && f < +std::numeric_limits<F>::min())
					return FP_SUBNORMAL;
				else if (std::numeric_limits<F>::lowest() <= f && f <= std::numeric_limits<F>::max())
					return FP_NORMAL;
				else if (f == f)
					return FP_INFINITE;
				else
					return FP_NAN;
			}

template<std::floating_point F>
			constexpr int ilogb(F f)
			{
				if (!internal_test && !std::is_constant_evaluated())
					return std::ilogb(f);

if (f == 0)
					return FP_ILOGB0;
				else if (std::numeric_limits<F>::lowest() <= f && f <= std::numeric_limits<F>::max())
				{
					if (f < 0)
						f = -f;

constexpr F fradix = std::numeric_limits<F>::radix;
					int exp = 0;
					if (f < 1)
						do { --exp; f *= fradix; } while (f < 1);
					else
						while (f >= fradix) { ++exp; f /= fradix; }

return exp;
				}
				else if (f == f)
					return INT_MAX;
				else
					return FP_ILOGBNAN;
			}

template<std::floating_point F>
			constexpr F scalbn(F f, int exp)
			{
				if (!internal_test && !std::is_constant_evaluated())
					return std::scalbn(f, exp);

if (f == 0 || !(std::numeric_limits<F>::lowest() <= f && f <= std::numeric_limits<F>::max()))
					return f;

constexpr F fradix = std::numeric_limits<F>::radix;
				if (exp < 0)
					do { f /= fradix; ++exp; } while (exp < 0);
				else
					while (exp > 0) { f *= fradix; --exp; }

return f;
			}
#endif
			constexpr bool float_is_binary32 =
				std::numeric_limits<float>::radix == 2 &&
				std::numeric_limits<float>::digits == 24 &&
				std::numeric_limits<float>::is_iec559;
			constexpr bool double_is_binary64 =
				std::numeric_limits<double>::radix == 2 &&
				std::numeric_limits<double>::digits == 53 &&
				std::numeric_limits<double>::is_iec559;

// spot check ilogb for some binary32 values
			static_assert(ilogb(0.0f) == FP_ILOGB0);
			static_assert(!float_is_binary32 || ilogb(FLT_TRUE_MIN) == -149);
			static_assert(!float_is_binary32 || ilogb(FLT_TRUE_MIN * 2) == -148);
			static_assert(!float_is_binary32 || ilogb(FLT_MIN) == -126);
			static_assert(!float_is_binary32 || ilogb(0.5f) == -1);
			static_assert(!float_is_binary32 || ilogb(1.0f - FLT_EPSILON / 2) == -1);
			static_assert(!float_is_binary32 || ilogb(1.0f) == 0);
			static_assert(!float_is_binary32 || ilogb(2.0f - FLT_EPSILON) == 0);
			static_assert(!float_is_binary32 || ilogb(2.0f) == 1);
			static_assert(!float_is_binary32 || ilogb(FLT_MAX) == 127);
		}

// unpacked floating point representation for easier manipulation
		template<int radix, int digits>
		struct decomposed_float
		{
			int m_category{};
			bool m_signbit{};
			std::array<int, digits> m_digits{};
			int m_exponent{};
		};

// decomposed_float compare is probably the most complicated function here!
		template<int radix, int digits>
		constexpr bool operator<(const decomposed_float<radix, digits>& lhs, const decomposed_float<radix, digits>& rhs)
		{
			if (lhs.m_category == FP_NAN || rhs.m_category == FP_NAN) // isnan(lhs) || isnan(rhs)
				return false;
			else if (lhs.m_category == FP_INFINITE || rhs.m_category == FP_INFINITE) // isinf(lhs) || isinf(rhs)
			{
				if (lhs.m_category == FP_INFINITE) // isinf(lhs)
					return lhs.m_signbit && !(rhs.m_category == FP_INFINITE && rhs.m_signbit); // lhs == -inf && rhs > -inf
				else // !isinf(lhs) && isinf(rhs)
					return !rhs.m_signbit; // lhs < +inf && rhs == +inf
			}
			else if (lhs.m_category == FP_ZERO || rhs.m_category == FP_ZERO) // lhs == 0 || rhs == 0
			{
				if (lhs.m_category == FP_ZERO) // lhs == 0
					return rhs.m_category != FP_ZERO && !rhs.m_signbit; // lhs == 0 && rhs > 0
				else // lhs != 0 && rhs == 0
					return lhs.m_signbit; // lhs < 0 && rhs == 0
			}
			else
			{
				// both subnormal or normal
				if (lhs.m_signbit != rhs.m_signbit)
					return lhs.m_signbit; // lhs < 0 && rhs > 0
				else
				{
					if (lhs.m_exponent != rhs.m_exponent)
						return lhs.m_signbit ? lhs.m_exponent > rhs.m_exponent : lhs.m_exponent < rhs.m_exponent;
					else
					{
						for (int i = 0; i < digits; ++i)
							if (lhs.m_digits[i] != rhs.m_digits[i])
								return lhs.m_signbit ? lhs.m_digits[i] > rhs.m_digits[i] : lhs.m_digits[i] < rhs.m_digits[i];
						return false; // all digits equal
					}
				}
			}
		}

// Helper alias for decomposing std::floating_point types
		template<std::floating_point F>
		using decomposed = decomposed_float<std::numeric_limits<F>::radix, std::numeric_limits<F>::digits>;

// decompose floating_point value
		template<std::floating_point F>
		constexpr decomposed<F> decompose(F x)
		{
			decomposed<F> df{};
			df.m_category = constexpr_cmath::fpclassify(x);
			df.m_signbit = constexpr_cmath::signbit(x);
			df.m_exponent = constexpr_cmath::ilogb(x);
			if (df.m_category == FP_NORMAL || df.m_category == FP_SUBNORMAL)
			{
				if (df.m_signbit)
					x = -x;
				x = constexpr_cmath::scalbn(x, -df.m_exponent);
				for (int d = 0; d < std::numeric_limits<F>::digits; ++d)
				{
					int digit = static_cast<int>(x);
					df.m_digits[d] = digit;
					x -= static_cast<F>(digit);
					if (x == 0)
						break;
					x = constexpr_cmath::scalbn(x, 1);
				}
			}
			return df;
		}

// compose floating_point value from decomposed equivalent
		template<std::floating_point F>
		constexpr F compose(const decomposed<F>& df)
		{
			F f = 0;
			if (df.m_category == FP_ZERO)
				f = 0;
			else if (df.m_category == FP_SUBNORMAL || df.m_category == FP_NORMAL)
			{
				for (int d = 0; d < std::numeric_limits<F>::digits; ++d)
				{
					F digit = static_cast<F>(df.m_digits[d]);
					f += constexpr_cmath::scalbn(digit, -d);
				}
				f = constexpr_cmath::scalbn(f, df.m_exponent);
			}
			else if (df.m_category == FP_INFINITE)
				f = std::numeric_limits<F>::infinity();
			else if (df.m_category == FP_NAN)
				f = std::numeric_limits<F>::quiet_NaN();

if (df.m_signbit)
				f = constexpr_cmath::copysign(f, F(-1));

return f;
		}

// convert integer value to decomposed<F>
		template<std::floating_point F, integer I>
		constexpr decomposed<F> decompose_integer(I x)
		{
			// digit extraction done in full precision, truncated if needed at the end to fit F
			using U = std::make_unsigned_t;
			constexpr int ibits = std::numeric_limits::digits + std::numeric_limits::is_signed;
			constexpr int tradix = std::numeric_limits<F>::radix;
			constexpr int tbits_per_digit = std::bit_width(static_cast<unsigned>(tradix)) - 1;
			constexpr int tdigits = ibits / tbits_per_digit + (0 != ibits % tbits_per_digit);
			std::array<int, tdigits> digits{};

int tsigdigits = 0;
			U u = x < 0 ? U(0) - static_cast(x) : x;
			for (int i = 0; i < tdigits && u != 0; ++i, ++tsigdigits)
			{
				digits[tdigits - 1 - i] = u % tradix;
				u /= tradix;
			}
			std::rotate(digits.begin(), digits.begin() + (tdigits - tsigdigits), digits.end());

constexpr int fdigits = std::numeric_limits<F>::digits;
			decomposed<F> df{};
			df.m_category = x ? FP_NORMAL : FP_ZERO;
			df.m_signbit = x < 0;
			std::copy(digits.begin(), digits.begin() + std::min(fdigits, tsigdigits), df.m_digits.begin());
			df.m_exponent = std::max(tsigdigits - 1, 0);
			return df;
		}

static_assert(compose<float>(decompose(-0.0f)) == -0.0f && constexpr_cmath::signbit(compose<float>(decompose(-0.0f))));
		static_assert(compose<float>(decompose(+0.0f)) == +0.0f && !constexpr_cmath::signbit(compose<float>(detail::decompose(+0.0f))));
		static_assert(compose<float>(decompose(-FLT_TRUE_MIN)) == -FLT_TRUE_MIN);
		static_assert(compose<float>(decompose(+FLT_TRUE_MIN)) == +FLT_TRUE_MIN);
		static_assert(compose<float>(decompose(-FLT_MIN)) == -FLT_MIN);
		static_assert(compose<float>(decompose(+FLT_MIN)) == +FLT_MIN);
		static_assert(compose<float>(decompose(-0.1f)) == -0.1f);
		static_assert(compose<float>(decompose(+0.1f)) == +0.1f);
		static_assert(compose<float>(decompose(-1.0f)) == -1.0f);
		static_assert(compose<float>(decompose(+1.0f)) == +1.0f);
		static_assert(compose<float>(decompose(-10.0f)) == -10.0f);
		static_assert(compose<float>(decompose(+10.0f)) == +10.0f);
		static_assert(compose<float>(decompose(-FLT_MAX)) == -FLT_MAX);
		static_assert(compose<float>(decompose(+FLT_MAX)) == +FLT_MAX);
	}

// in_range<integer>(floating_point) implementation
	template<integer I, std::floating_point F>
	constexpr bool in_range(F f)
	{
		constexpr detail::decomposed<F> min_in_range_d = std::max(
			detail::decompose<F>(std::numeric_limits<F>::lowest()),
			detail::decompose_integer<F>(std::numeric_limits::lowest()));
		constexpr detail::decomposed<F> max_in_range_d = std::min(
			detail::decompose<F>(std::numeric_limits<F>::max()),
			detail::decompose_integer<F>(std::numeric_limits::max()));
		constexpr F min_in_range = detail::compose<F>(min_in_range_d);
		constexpr F max_in_range = detail::compose<F>(max_in_range_d);
		return min_in_range <= f && f <= max_in_range;
	}

#ifdef UINT32_MAX
	static_assert(!detail::constexpr_cmath::float_is_binary32 || !in_range<uint32_t>(float(UINT32_MAX)));
	static_assert(!detail::constexpr_cmath::double_is_binary64 || in_range<uint32_t>(double(UINT32_MAX)));
#endif
#ifdef INT32_MAX
	static_assert(!detail::constexpr_cmath::float_is_binary32 || in_range<int32_t>(float(INT32_MIN)));
	static_assert(!detail::constexpr_cmath::float_is_binary32 || !in_range<int32_t>(float(INT32_MAX)));
	static_assert(!detail::constexpr_cmath::double_is_binary64 || in_range<int32_t>(double(INT32_MIN)));
	static_assert(!detail::constexpr_cmath::double_is_binary64 || in_range<int32_t>(double(INT32_MAX)));
#endif

// clamp(floating_point, integer, integer) implementation
	template<integer I, std::floating_point F>
	constexpr I clamp(F f, I min, I max)
	{
		// What to do with NaN?
		if (!std::is_constant_evaluated())
			if (std::isnan(f))
				throw std::domain_error("isnan(f)");

I i = in_range(f) ? static_cast(f) :
			f < 0 ? std::numeric_limits::min() :
			std::numeric_limits::max();

return std::clamp(i, min, max);
	}

static_assert(clamp(float(INT_MAX), INT_MIN, INT_MAX) == INT_MAX);
	static_assert(clamp(float(INT_MIN), INT_MIN, INT_MAX) == INT_MIN);
	static_assert(clamp(FLT_MAX, -1, 1) == 1);
	static_assert(clamp(FLT_MIN, -1, 1) == 0);
	static_assert(clamp(-FLT_MIN, -1, 1) == 0);
	static_assert(clamp(-FLT_MAX, -1, 1) == -1);
}

// Non-constexpr / out-of-line versions to check code generation.
bool in_int_range(float f)
{
	return in_range_ext::in_range<int>(f);
}

int clamp_to_int(float f, int min, int max)
{
	return in_range_ext::clamp(f, min, max);
}

int main()
{
}