Sorting Floating Point Values

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 / EXIT_SUCCESS qsort /
#include / printf /
#include / assert /
#include / clock /
#include / INT_MAX /
#include / C99 floating point macros /
#include / 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 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,

-
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

-
, 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;

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

-
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);
}

-
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

StackExchange Code Review Q#161154, answer score: 4