[C++][c++20] std::midpoint

Reference:
https://devblogs.microsoft.com/oldnewthing/20220207-00/?p=106223
https://ai.googleblog.com/2006/06/extra-extra-read-all-about-it-nearly.html
https://en.cppreference.com/w/cpp/numeric/midpoint
https://www.youtube.com/watch?v=sBtAGxBh-XI

Types for std::midpoint

1. integral
   5 of them:
   signed char, short int, int, long int, long long int https://en.cppreference.com/w/cpp/language/types
   signed overflow is UB; unsigned overflow is defined to wrap around.
   signed converted to unsigned is always safe.
   every signed integral type has a corresponding unsigned integral type and its name is std::make_unsigned_t<T>
2. pointer
3. floating point
   Denormalized numbers, INF, and NaN
   "at most one inexact operation"

No, don't do this:

unsigned average(unsigned a, unsigned b)
{
    return (a + b) / 2;
}

Better if we know which value is larger:

unsigned average(unsigned low, unsigned high)
{
    return low + (high - low) / 2;
}

Better:

unsigned average(unsigned low, unsigned high)
{
	return ((unsigned int)low + (unsigned int)high)) >> 1;
}

Better:

unsigned average(unsigned a, unsigned b)
{
    return (a / 2) + (b / 2) + (a & b & 1);
}

Better:

unsigned average(unsigned a, unsigned b)
{
    return (a & b) + (a ^ b) / 2;
}

If your compiler supports integers larger than the size of an unsigned, say because unsigned is a 32-bit value but the native register size is 64-bit, or because the compiler supports multiword arithmetic, then you can cast to the larger data type:

unsigned average(unsigned a, unsigned b)
{
    // Suppose "unsigned" is a 32-bit type and
    // "unsigned long long" is a 64-bit type.
    return ((unsigned long long)a + b) / 2;
}

code snipped from talk: CppCon 2019: Marshall Clow “std::midpoint? How Hard Could it Be?”

constexpr Integer midpoint(Integer a, Integer b) noexcept {
	using U = std::make_unsigned_t<Integer>;

	int sign = 1;
	U m = a;
	U M = b;

	if (a > b) {  // no branch generated.
		sign = -1;
		m = b;
		M = a;
	}

	return a + sign * Integer(U(M-m) >> 1); // C casting works while we are sure input is integral types.
}


template<typename T>
constexpr enable_if_t<is_pointer_v<T>, T>
midpoint_ptr(T a, T b) {
	return a + midpoint(ptrdiff_t{0}, b - a); // std::ptrdiff_t is the signed integer type of the result of subtracting two pointers.
}

template<typename T>
constexpr enable_if_t<is_floating_point_v<Fp>, Fp>
midpoint_fp(T a, T b) noexcept {
	Fp lo = numeric_limits<Fp>::min()*2;
	Fp hi = numeric_limits<Fp>::max()*2;

	return abs(a) <= hi && abs(b) <= hi ? // typical case
		(a + b) / 2: // alwas correctly rounded
		abs(a) < lo ? a + b/2: // not safe to halve a
		abs(b) < lo ? a/2 + b: // not safe to halve b
			a/2 + b/2; // otherwise correctly rounded
}