One approach is to enumerate all numeric floating-point encodings for a given floating-point type, then find the enumeration that is half-way between the enumerations of the two source operands. This answer assumes the use of a platform that uses IEEE-754 floating-point types and twos-complement integers.
The bit pattern of a floating-point operand itself can provide the desired enumeration. For example, for a double operand in IEEE-754 binary64 format, we copy the bits to a 64-bit unsigned integer, then reinterpret that as a signed 64-bit integer. However, before we re-interpret as a signed integer, we must take care of negative operands by reflecting them such negative zero is mapped onto positive zero, and negative infinity maps to a negative integer large in magnitude.
Now all possible numerical floating-point operands are contiguously enumerated in numerical order as a signed integer. If we do this for both source operands, then compute the average of the two enumerations, the resulting enumeration represents the desired midpoint. This midpoint is first transformed into an unsigned 64-bit integer, again reflecting negative operands. The unsigned 64-bit integer is copied to a double in IEEE-754 binary64 format.
The ISO-C program below implements the above algorithm. One caveat: Because IEEE-754 bindings are optional in C, not all compilers may handle -0.0 as a number distinct from 0.0; this may require at minimum compiling with the strictest floating-point setting of the compiler (e.g. -fp-model:strict). Alternative ways of constructing FP64_NEG_ZERO as an unsigned integer of same size as the floating-point type with only its most-significant bit set are possible, of course.
#include <stdio.h>
#include <stdlib.h>
#include <stdint.h>
#include <string.h>
#include <float.h>
#include <limits.h>
#include <math.h>
#define MORE_EFFICIENT (1)
#if MORE_EFFICIENT
/* Henry S. Warren Jr., "Hacker's Delight 2nd ed.", Addison-Wesley 2013, p. 19 */
/* arithmetic right shift by 1 bit */
int64_t asr1 (int64_t x)
{
return (x & ~1) / 2;
}
/* compute average of signed integers w/o overflow in intermediate computation*/
int64_t avg (int64_t x, int64_t y)
{
int64_t t = (x & y) + asr1 (x ^ y);
return t + (int64_t)(((uint64_t)t >> 63) & (x ^ y));
}
#else // MORE_EFFICIENT
/* compute average of signed integers w/o overflow in intermediate computation*/
int64_t avg (int64_t a, int64_t b)
{
int64_t r;
if ((a < 0) == (b < 0)) {
r = a / 2 + b / 2 + (a % 2 + b % 2) / 2;
} else {
r = (a + b) / 2;
}
return r;
}
#endif // MORE_EFFICIENT
/* reinterpret bit pattern of IEEE-754 'binary64' as unsigned 64-bit integer */
uint64_t double_as_uint64 (double a)
{
uint64_t r;
memcpy (&r, &a, sizeof r);
return r;
}
/* reinterpret bit pattern of unsigned 64-bit integer as IEEE-754 'binary64' */
double uint64_as_double (uint64_t a)
{
double r;
memcpy (&r, &a, sizeof r);
return r;
}
/* map an IEEE-754 'binary64' encoding to its enumeration position */
int64_t double_to_enumeration (double a)
{
const uint64_t FP64_NEG_ZERO = double_as_uint64 (-0.0);
uint64_t t = double_as_uint64 (a);
if ((int64_t)t < 0) {
t = FP64_NEG_ZERO - t;
}
return (int64_t)t;
}
/* map enumeration position to IEEE-754 'binary64' encoding */
double enumeration_to_double (int64_t a)
{
const uint64_t FP64_NEG_ZERO = double_as_uint64 (-0.0);
uint64_t t = a;
if (a < 0) {
t = FP64_NEG_ZERO - t;
}
return uint64_as_double (t);
}
/* find midpoint between two IEEE-754 'binary64' numbers a and b, such that
there is as equal a number of discrete encodings as possible between a
and the midpoint, and the midpoint and b.
*/
double midpoint (double a, double b)
{
int64_t ia = double_to_enumeration (a);
int64_t ib = double_to_enumeration (b);
int64_t midpoint_i = avg (ia, ib);
return enumeration_to_double (midpoint_i);
}
/* try test cases from question */
int main (void)
{
double lb [6] = {
-INFINITY,
1.7976931348623155e308,
-INFINITY,
0.0,
-5.0e-320,
-5.0e-315
};
double ub [6] = {
INFINITY,
INFINITY,
1.7976931348623155e308,
INFINITY,
1.0e-318,
1.000004e-318
};
printf ("midpoint of % 23.16e and % 23.16e is % 23.16e\n", lb[0], ub[0], midpoint (lb[0], ub[0]));
printf ("midpoint of % 23.16e and % 23.16e is % 23.16e\n", lb[1], ub[1], midpoint (lb[1], ub[1]));
printf ("midpoint of % 23.16e and % 23.16e is % 23.16e\n", lb[2], ub[2], midpoint (lb[2], ub[2]));
printf ("midpoint of % 23.16e and % 23.16e is % 23.16e\n", lb[3], ub[3], midpoint (lb[3], ub[3]));
printf ("midpoint of % 23.16e and % 23.16e is % 23.16e\n", lb[4], ub[4], midpoint (lb[4], ub[4]));
printf ("midpoint of % 23.16e and % 23.16e is % 23.16e\n", lb[5], ub[5], midpoint (lb[5], ub[5]));
return EXIT_SUCCESS;
}
The program's output should look similar to this:
midpoint of -1.#INF000000000000e+000 and 1.#INF000000000000e+000 is 0.0000000000000000e+000
midpoint of 1.7976931348623155e+308 and 1.#INF000000000000e+000 is 1.7976931348623157e+308
midpoint of -1.#INF000000000000e+000 and 1.7976931348623155e+308 is -4.9406564584124654e-324
midpoint of 0.0000000000000000e+000 and 1.#INF000000000000e+000 is 1.5000000000000000e+000
midpoint of -4.9999443359134150e-320 and 9.9999874849559983e-319 is 4.7499965256823284e-319
midpoint of -5.0000000022897320e-315 and 1.0000036891520582e-318 is -2.4994999993002900e-315
