I recently read a paper¹ about the difficulty of calculating the most exact midpoint between two floating point values, facing overflow, underflow and rounding errors. The intention is to assure that the midpoint(a,b) is

- actually within the interval [a,b]
- the floating point number in that interval that is closest to the real midpoint (i.e. (a+b)/2).

It feels like a function that should be part of Python's math module.

¹ https://hal.archives-ouvertes.fr/file/index/docid/576641/filename/computing-midpoint.pdf

The author proposes the following implementation (pages 20/21):

midpoint(a,b) = 
    a  if a == b else             # covers midpoint(inf, inf)
    0  if a == -b else            # covers midpoint(-inf, +inf)
    -realmax  if a == -inf else   # min. double value
    +realmax  if b == +inf else   # max. double value
    round_to_nearest_even((a - a/2) + b/2)

I guess nans should simply pass through as in other functions, i.e. midpoint(a,nan) == midpoint(nan,b) == nan.

The behaviour for [a,inf] is decidedly up for discussion. There are certainly cases where midpoint(a,+inf) would best return +inf, but +realmax as an actual finite value also seems reasonable. OTOH, it's easy for users to promote inf->realmax or realmax->inf or even inf->a themselves, just like it's easy to check for +/-inf before calling the function. It just takes a bit longer to do these checks on user side. There could also be a "mode" argument that makes it return one of: a or b (i.e. the finite bound), +/-realmax or +/-inf in the two half-infinity cases.

What do you think about this addition?