Skip to main content
Code Review

Return to Question

replaced http://stackoverflow.com/ with https://stackoverflow.com/
Source Link

I found sorting large arrays by a comparator that looks at the floating-point difference problematic, especially -ffast-math, (http://stackoverflow.com/questions/24442725/is-floating-point-addition-commutative-in-c. 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:

#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{ httphttps://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;
}

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 Radix Sort, Sorting a float data.

I found sorting large arrays by a comparator that looks at the floating-point difference problematic, especially -ffast-math, (http://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:

#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{ http://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;
}

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.

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:

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

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.

Put it into a list; fixed a most glaring error pointed out by chux; charts.
Source Link
Neil
  • 1.1k
  • 8
  • 20

I found sorting large arrays by a comparator that looks at the floating-point difference problematic, especially -ffast-math, (http://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/ .sort:

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{ http://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;
 {
 const/* @chux unsigned -> uint32_t critical fix */
 const uint32_t ax_abs = ax.u & 0x7fffffff, ax_sign = ax.u & 0x80000000;
 const unsigneduint32_t bx_abs = bx.u & 0x7fffffff, bx_sign = bx.u & 0x80000000;
 const unsigneduint32_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. int32_t

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 .

IEEE 754 floating point.

I found sorting large arrays by a comparator that looks at the floating-point difference problematic, especially -ffast-math, (http://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{ http://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;
 {
 const unsigned ax_abs = ax.u & 0x7fffffff, ax_sign = ax.u & 0x80000000;
 const unsigned bx_abs = bx.u & 0x7fffffff, bx_sign = bx.u & 0x80000000;
 const unsigned 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?

I found sorting large arrays by a comparator that looks at the floating-point difference problematic, especially -ffast-math, (http://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:

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{ http://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. int32_t

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 .

IEEE 754 floating point.

Source Link
Neil
  • 1.1k
  • 8
  • 20

Sorting Floating Point Values

I found sorting large arrays by a comparator that looks at the floating-point difference problematic, especially -ffast-math, (http://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{ http://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;
 {
 const unsigned ax_abs = ax.u & 0x7fffffff, ax_sign = ax.u & 0x80000000;
 const unsigned bx_abs = bx.u & 0x7fffffff, bx_sign = bx.u & 0x80000000;
 const unsigned 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?

lang-c

AltStyle によって変換されたページ (->オリジナル) /