# Bisecting Intervals of floating point numbers containing 0 and infinity fairly

It is seldom considered that floating points are not evenly distributed in the real number line. I've been working with interval arithmetic and noticed when bisecting $$[a,b]$$ on the real number line into $$[a,(a+b)/2]$$ and $$[(a+b)/2,b]$$ starting from $$[-1.3552646593344489*10^{308} ,0]$$ until $$a,b$$ are neighboring floats takes way more steps when you always pick the interval containing $$-10$$ and way fewer steps if you always pick the interval containing $$-10^{305}$$. This is because taking the average of two floating point (fp) numbers with differing exponents divides a floating interval in sub intervals which contain a very different number of floats.

Given two interval boundaries $$(a,b)$$ how do i calculate a fp number $$m$$ such that between $$a$$ and $$m$$ there are nearly as many fp numbers as between $$m$$ and $$b$$ given that $$a,b$$ can have different signs or have values from $$\{-\infty,\infty\}$$?

A design for 64 bit floating point numbers would be welcome. Please document your design decisions so i can reimplement those for other precisions. Note that when dividing interval boxes the unevenness gets amplified exponentially with the number of dimensions so i am to trade speed for getting a better $$m$$.

### Implementation that works if the signs of $$a,b$$ are the same

I am already aware of Implementations like _middle(x1::Float64, x2::Float64) but it doesn't calculate the correct points if signs of $$a$$ and $$b$$ differ or if atleast one of them is infinite. Infinity should be treated as a floating point number one step larger then the largest finite floating point number.

### Test cases

• The midpoint of $$[-\infty,\infty]$$ should be 0.
• The midpoint of $$[1.7976931348623155*10^{308},\infty]$$ (the second largest finite float and prevfloat(prevfloat(Inf)) ) should be 1.7976931348623157e308 (the largest finite float and prevfloat(Inf)).
• The midpoint of $$[-\infty,1.7976931348623155*10^{308}]$$ should be -5.0e-324 (the largest negative float and prevfloat(0.0)).
• The midpoint of $$[0.0, \infty]$$ should be 1.5.
• The midpoint of $$[-5.0*10^{-320}, 1.0*10^{-318}]$$ should be $$4.75*10^{-319}$$.
• The midpoint of $$[-5.0*10^{-315}, 1.000004*10^{-318}]$$ should be $$-2.4995*10^{-315}$$.

### Red Herring

The publication How do you compute the midpoint of an interval? compares different numerical approaches to divide in a $$\frac{a+b}{2}$$ way which doesn't have the property i require.

• @D.W. Is my inquery more clear know? Jan 3, 2022 at 21:57
• That's a big improvement, thank you. Something seems missing: you write "...as a nextfloat precise interval and does so..." but I can't make any sense of that sentence. Are some words missing? I can't parse the grammar. Can you break it down into smaller sentences? Also, what's the meaning of "nextfloat precise"/"nextfloat precision"?
– D.W.
Jan 3, 2022 at 22:02
• Are you saying that you know of a solution for this problem that works except it doesn't handle certain special cases? Which special cases does it not handle, and have you tried working through each of those special cases to identify what is the desired/correct outcome for each of them? If so, can you identify the special cases for which you do know what the correct output should be, and for which ones you don't know what the correct output should be? I imagine if you know what is the desired output for each special case, you should be able to write code that checks for each such case.
– D.W.
Jan 3, 2022 at 22:06
• @D.W. i provided some test cases and pointed out what the implementation doesn't handle. Jan 4, 2022 at 15:21

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

• If you elaborate on ((uint64_t)1) << (CHAR_BIT * sizeof (a) - 1) and whether -0.0 would serve the same purpose i think i can accept this answer. Mar 26, 2022 at 14:21
• @worldsmithhelper Unfortunately, in my experience, one cannot rely on compilers to handle -0.0 correctly. Otherwise we could write double_as_uint64 (-0.0). The alternative is to do the work ourselves and write 0x8000000000000000ULL, assuming unsigned long long int is a 64-bit type. Is it obvious what needs to be changed when switching to float and uint32_t? I thought maybe not, and chose the portable expression that lets you search and replace double with float and 64 with 32 and be done. What do you consider preferable? Mar 26, 2022 at 18:23