Sorting Floating Point Values

Question

I found sorting large arrays by a comparator that looks at the floating-point difference problematic, especially -ffast-math, (https://stackoverflow.com/questions/24442725/is-floating-point-addition-commutative-in-c.) This is intended to compare 32-bit floats exactly using integer arithmetic. Most references cited concern themselves with radix sort:

• Nicholas Chapman, https://www.forwardscattering.org/post/34;
• Michael Herf, http://stereopsis.com/radix.html;
• Pierre Terdiman, http://codercorner.com/RadixSortRevisited.htm;
• https://hbfs.wordpress.com/2010/03/09/radix-sort-on-floating-point-numbers/;
• https://randomascii.wordpress.com/category/floating-point/.

Using a comparison sort, like qsort, it only has to decide on a total order of IEEE-754 numbers by comparing two at a time.

#include <stdlib.h> /* EXIT_SUCCESS qsort */
#include <stdio.h>  /* printf */
#include <assert.h> /* assert */
#include <time.h>   /* clock */
#include <limits.h> /* INT_MAX */
#include <math.h>   /* C99 floating point macros */
#include <stdint.h> /* C99 uint32_t */

struct Foo { float x; };

/** Compares float {x} values of {Foo} exactly. Assumes IEEE-754-ish 32-bit
 float has the same endianness as {uint32_t}. References:
 \cite{KimYoonKim2011FastSortFloatingPoint},
 Nicholas Chapman \url{ https://www.forwardscattering.org/post/34 },
 Michael Herf \url{ http://stereopsis.com/radix.html },
 Pierre Terdiman \url{ http://codercorner.com/RadixSortRevisited.htm },
 \url{ https://hbfs.wordpress.com/2010/03/09/radix-sort-on-floating-point-numbers/ },
 \url{ https://randomascii.wordpress.com/category/floating-point/ },
 \url{ https://stackoverflow.com/questions/10632237/any-c-compiler-where-evaluates-to-larger-than-one }.
 @implements <Foo>Comparator
 @return Greater than, equal to, or less than 0, if the {Foo.x} pointed to by
 {av} is greater than, equal to, or less than the {Foo.x} pointed to by {bv}. */
static int x_cmp(const void *av, const void *bv) {
    const struct Foo *const a = av, *const b = bv;
    union { float f; uint32_t u; } ax, bx;
    ax.f = a->x, bx.f = b->x;
    {
        /* @chux unsigned -> uint32_t critical fix */
        const uint32_t ax_abs = ax.u & 0x7fffffff, ax_sign = ax.u & 0x80000000;
        const uint32_t bx_abs = bx.u & 0x7fffffff, bx_sign = bx.u & 0x80000000;
        const uint32_t same = !(ax_sign ^ bx_sign);
        /* if(!same) return ax_sign^1; else return ax_sign^(ax_abs - bx_abs); */
        return ax_sign ^ (!same + same * (ax_abs - bx_abs));
    }
}

/** Tests {x_cmp}. */
int main(void) {
    struct Foo foos[32]; /* could be bigger */
    const size_t foos_size = sizeof foos / sizeof *foos;
    unsigned i;
    union { float f; uint32_t u; } x;

    assert(sizeof(float) == 4);
    srand((unsigned)clock()), rand();
    foos[0].x = 0.0f;
    foos[1].x = 0.0f;
    foos[2].x = -0.0f;
    foos[3].x = INFINITY; /* C99 */
    foos[4].x = -INFINITY; /* C99 */
    foos[5].x = NAN; /* C99 */
    foos[6].x = (x.u = 0x807fffff, x.f); /* subnormal */
    foos[7].x = (x.u = 1, x.f); /* subnormal */
    for(i = 8; i < 2 * foos_size / 3; i++) {
        foos[i].x = rand() / (0.5f * RAND_MAX / INT_MAX) - INT_MAX;
    } for( ; i < foos_size; i++) {
        foos[i].x = rand() / (0.5f * RAND_MAX / 100.0f) - 100.0f;
    }
    qsort(foos, foos_size, sizeof foos[0], &x_cmp);
    for(i = 1; i < foos_size; i++) {
        printf("%g <= %g\n", foos[i - 1].x, foos[i].x);
        /* isnan C99 */
        assert(isnan(foos[i - 1].x) || isnan(foos[i].x)
            || foos[i - 1].x <= foos[i].x);
    }
    return EXIT_SUCCESS;
}

I compiled it on GCC 4.2.1 MacOS, GCC 5.4.0 Linux, and MSVC15, (the order of the NaN macro switched, thus isnan,) I think it's valid? Could the comparator be done in less steps?

---

Edit:

int32_t: increasing injection mod the branch cut.

IEEE 754-ish float: (one way) to get rid of the singularity and make it monotonic is to invert the negative values and flip the sign bit on positive values, Radix Sort, Sorting a float data.

Answer 1

1. ax_sign ^ (!same + same * (ax_abs - bx_abs)); returns the incorrect signed result if int is not 32-bit. Certainly a problem if int is 16-bit and likely if int is 64-bit. If unsigned/int needs to be 32-bit, use (u)int32_t

2. , operator reduces clarity here. Suggest 2 lines of code

   // ax.f = a->x, bx.f = b->x;
   ax.f = a->x;
   bx.f = b->x;
   

3. I ran test cases and found OP's x_cmp() functionally correct for finite float over a 2,000,000,000 test cases.

4. As a test case, I tried OP's original x_cmp() versus the below and was at least 10% faster with the new code. Of course, that is just one platform comparison, yet aside from NaN issues, the below code is functionally similar to OP's and as a plus, is highly portable - unlike OP's. The point being that OP's compare method needs some reference point to justify the bit magic.

   static int x_cmp_ref(const void *av, const void *bv) {
     return (*(float*)av > *(float*)bv) - (*(float*)av < *(float*)bv);
   }
   

5. OP's has not stated the compare functionality of Not-a-number floats. A desirable aspect is that all NaN sort to one side, either all greater or all less than any other number, regardless of the NaN's "sign". x_cmp() considers sign first without regard to NaN-ness.

---

As a reference, I used the following to generate random float

float randf() {
  union {
      float f;
      unsigned char uc[sizeof (float)];
  } u;
  do {
    for (unsigned i=0; i<sizeof u.uc; i++) {
      u.uc[i] = (unsigned char) rand();
    }
  } while (!isfinite(u.f));
  return u.f;
}

---

At a later time, I may try to implement the idea in this comment