CODEWITHSUNDEEP
bit-manipulationieee-754
Simon
Simon

Reputation: 320

compute midpoint in floating point

Given two floating point numbers (IEEE single or double precision), I would like to find the number that lies half-way between them, but not in the sense of (x+y)/2 but with respect to actually representable numbers.

if both x and y are positive, the following works

float ieeeMidpoint(float x, float y)
{
    assert(x >= 0  && y >= 0);
    int xi = *(int*)&x;
    int yi = *(int*)&y;
    int zi = (xi+yi)/2;
    return *(float*)&zi;
}

The reason this works is that positive ieee floating point numbers (including subnormals and infinity) keep their order when doing a reinterpreting cast. (this is not true for the 80-bit extended format, but I don't need that anyway).

Now I am looking for an elegant way to do the same that includes the case when one or both of the numbers are negative. Of course it is easy to do with a bunch of if's, but I was wondering if there is some nice bit-magic, prefarably without any branching.

Upvotes: 0

Views: 232

Answers (1)

Simon
Simon

Reputation: 320

Figured it out myself. the order of negative number is reversed when doing the reinterpreting cast, so that is the only thing one needs to fix. This version is longer than I hoped it would be, but its only some bit-shuffling, so it should be fast.

float ieeeMidpoint(float x, float y)
{
    // check for NaN's (Note that subnormals and infinity work fine)
    assert(x ==x && y == y);

    // re-interpreting cast
    int xi = *(int*)&x;
    int yi = *(int*)&y;

    // reverse negative numbers
    // (would look cleaner with an 'if', but I like not branching)
    xi = xi ^ ((xi >> 31) & 0x3FFFFFFF);
    yi = yi ^ ((yi >> 31) & 0x3FFFFFFF);

    // compute average of xi,yi (overflow-safe)
    int zi = (xi >> 1) + (yi >> 1) + (xi & 1);

    // reverse negative numbers back
    zi = zi ^ ((zi >> 31) & 0x3FFFFFFF);

    // re-interpreting back to float
    return *(float*)&zi;
}

Upvotes: 1

Related Questions