CppUtils/modules/Math/Utility.mpp at master · MorganCaron/CppUtils · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
export module CppUtils.Math.Utility;

import std;
import CppUtils.Type.Concept;

export namespace CppUtils::Math
{
	// Plus utile en C++23 (std::abs constexpr)
	[[nodiscard]] inline constexpr auto absolute(auto number) -> decltype(auto)
	{
		if constexpr (std::unsigned_integral<decltype(number)>)
			return number;
		else if constexpr (std::floating_point<decltype(number)>)
			return number < 0 ? -number : number;
		else
			return number < 0 ?
				static_cast<std::make_unsigned_t<decltype(number)>>(-number) :
				static_cast<std::make_unsigned_t<decltype(number)>>(number);
	}

	static_assert(absolute(0) == 0u);
	static_assert(absolute(1) == 1u);
	static_assert(absolute(-1) == 1u);
	static_assert(absolute(0.0) == 0.0);
	static_assert(absolute(1.5) == 1.5);
	static_assert(absolute(-1.5) == 1.5);
	static_assert(absolute(0u) == 0u);
	static_assert(absolute(1u) == 1u);

	/**
	 * @brief Performs robust, magnitude-adaptive floating-point equality comparison.
	 *
	 * Direct comparison of floating-point numbers with `operator==` is notoriously unreliable
	 * due to IEEE 754 binary representation errors and rounding artifacts. For instance, the expression
	 * `(0.1 + 0.2) == 0.3` evaluates to `false` because in base-2 floating point, `0.1 + 0.2` results
	 * in `0.30000000000000004440892098500626...`.
	 *
	 * Furthermore, naive approaches like `std::abs(a - b) < fixed_epsilon` fail for two major reasons:
	 * 1. For large numbers, a fixed epsilon is far too tight (the gap between representable numbers
	 *    exceeds the epsilon), causing false negatives.
	 * 2. Near zero (subnormal numbers), relative scaling collapses, causing underflow issues.
	 *
	 * `isEqual` solves both problems by implementing an adaptive tolerance algorithm:
	 * `|lhs - rhs| <= std::max(minimumEpsilon, epsilon * std::max(|lhs|, |rhs|))`
	 *
	 * - **Adaptive scaling**: Scales tolerance proportionally to the operands' magnitude (`epsilon * max`),
	 *   maintaining accurate relative precision across small and arbitrarily large numbers.
	 * - **Zero-proximity fallback**: Protects comparisons near zero by falling back to `minimumEpsilon`
	 *   (`std::numeric_limits<Float>::min()` by default).
	 *
	 * @tparam Float A floating-point type satisfying `std::floating_point`.
	 * @param lhs Left-hand operand to compare.
	 * @param rhs Right-hand operand to compare.
	 * @param epsilon Relative tolerance factor (defaults to machine epsilon `std::numeric_limits<Float>::epsilon()`).
	 * @param minimumEpsilon Absolute tolerance floor for values close to zero (defaults to `std::numeric_limits<Float>::min()`).
	 * @return `true` if the difference between operands falls within the adaptive tolerance, `false` otherwise.
	 */
	template<std::floating_point Float>
	[[nodiscard]] inline constexpr auto isEqual(
		Float lhs,
		Float rhs,
		Float epsilon = std::numeric_limits<Float>::epsilon(),
		Float minimumEpsilon = std::numeric_limits<Float>::min()) noexcept -> bool
	{
		[[assume(epsilon > 0), assume(minimumEpsilon > 0)]];
		auto difference = absolute(lhs - rhs);
		auto maximum = std::ranges::max(absolute(lhs), absolute(rhs));
		return difference <= std::ranges::max(minimumEpsilon, epsilon * maximum);
	}

	static_assert(isEqual(0.0, 0.0));
	static_assert(isEqual(0.1 + 0.2, 0.3));
	static_assert(not isEqual(0.1 + 0.2, 0.29));
	static_assert(not isEqual(0.1 + 0.2, 0.31));
	static_assert(isEqual(1e-308, 0.0));
	static_assert(not isEqual(1e-307, 0.0));
	static_assert(not isEqual(std::numeric_limits<double>::epsilon(), 0.0));

	template<Type::Numeric T>
	[[nodiscard]] inline constexpr auto isBetween(T value, T low, T high) noexcept -> bool
	{
		return low <= value and value <= high;
	}

	static_assert(isBetween(5, 0, 10));
	static_assert(not isBetween(-5, 0, 10));
	static_assert(not isBetween(15, 0, 10));
	static_assert(isBetween(5.0, 0.0, 10.0));
	static_assert(not isBetween(-5.0, 0.0, 10.0));
	static_assert(not isBetween(15.0, 0.0, 10.0));
	static_assert(isBetween(5.0f, 0.0f, 10.0f));
	static_assert(not isBetween(-5.0f, 0.0f, 10.0f));
	static_assert(not isBetween(15.0f, 0.0f, 10.0f));
}