I posted an earlier version of this on Tumblr. It had bugs (namely, it ignored signedness). Here’s an improved version.

#include <math.h>
/* Test (min / max) for closeness to 1 (i.e., what fraction of
* 'max' is 'min'?). And since (min / max >= cutoff) is equivalent
* to (min >= cutoff * max), we can avoid division. (I'm reminded
* of a young linear interpolation...)
*
* The best way to craft the cutoff constant is to subtract epsilon
* from 1.0, with epsilon defined in terms of the smallest value
* with exponent 0 (i.e., the ULP of 1.0). Here, we fake it.
*
* NOTE: .99999... approaches 1. Thou shall not Proliferate thy Nines
* lest thou incur the Confusion of the Internal strtod() of the Dread
* Compiler (Blessed be Her Grammar).
*
* NOTE: The approximation threshold scales. As the values approach
* zero, so does the threshold. In other words, a positive can never
* approximately equal a negative.
*
* NOTE: This has not been thoroughly tested. Approximate function
* is approximate.
*/
int
approx(double a, double b)
{
#define CUTOFF .9999999 /* 15-17 significant digits; stop halfway */
double da = fabs(a);
double db = fabs(b);
int samesigns = (da == a) == (db == b); /* both pos or neg */
if (da > db) { double tmp = da; da = db; db = tmp; }
return samesigns && da >= (db * CUTOFF);
}
#ifdef TEST
#include <stdio.h>
void
test_approx(void)
{
struct {
int rv;
double a;
double b;
} tests[] = {
{ 0, 1, 1.1 },
{ 1, 1, 1.00000005 },
{ 0, -1, 1.00000005 },
{ 0, 1, -1.00000005 },
{ 1, -1, -1.00000005 },
{ 0, .000000001, .000000005 },
{ 0, -.000000001, .000000005 },
{ 0, .000000001, -.000000005 },
{ 0, -.000000001, -.000000005 },
{ 0, 1, 2 },
{ 0, 0.5, 0.6 },
{ 0, 0.000005, 0.000006 },
{ 0, 1e100, 2e100 },
{ 1, .9999999999e0, 1e0 },
{ 1, .9999999999e1, 1e1 },
{ 1, .9999999999e10, 1e10 },
{ 1, .9999999999e20, 1e20 },
{ 1, .9999999999e30, 1e30 },
{ 1, .9999999999e31, 1e31 },
{ 1, .9999999999e32, 1e32 },
{ 1, .9999999999e33, 1e33 },
{ 1, .9999999999e50, 1e50 },
{ 1, .9999999999e62, 1e62 },
{ 1, .9999999999e63, 1e63 },
{ 1, .9999999999e64, 1e64 },
{ 1, .9999999999e65, 1e65 },
{ 1, .9999999999e100, 1e100 },
{ 1, .9999999999e200, 1e200 },
{ 1, .9999999999e300, 1e300 },
{ 1, .9999999999e-1, 1e-1 },
{ 1, .9999999999e-2, 1e-2 },
{ 1, .9999999999e-10, 1e-10 },
{ 1, .9999999999e-20, 1e-20 },
{ 1, .9999999999e-30, 1e-30 },
{ 1, .9999999999e-40, 1e-40 },
{ 1, .9999999999e-100, 1e-100 },
{ 1, .9999999999e-200, 1e-200 },
{ 1, .9999999999e-300, 1e-300 },
{ 0, 2e0, 1e0 },
{ 0, 2e1, 1e1 },
{ 0, 2e10, 1e10 },
{ 0, 2e20, 1e20 },
{ 0, 2e30, 1e30 },
{ 0, 2e31, 1e31 },
{ 0, 2e32, 1e32 },
{ 0, 2e33, 1e33 },
{ 0, 2e50, 1e50 },
{ 0, 2e62, 1e62 },
{ 0, 2e63, 1e63 },
{ 0, 2e64, 1e64 },
{ 0, 2e65, 1e65 },
{ 0, 2e100, 1e100 },
{ 0, 2e200, 1e200 },
{ 0, 2e300, 1e300 },
{ 0, 2e-1, 1e-1 },
{ 0, 2e-2, 1e-2 },
{ 0, 2e-10, 1e-10 },
{ 0, 2e-20, 1e-20 },
{ 0, 2e-30, 1e-30 },
{ 0, 2e-40, 1e-40 },
{ 0, 2e-100, 1e-100 },
{ 0, 2e-200, 1e-200 },
{ 0, 2e-300, 1e-300 },
{ -1, 0.0, 0.0 }
};
size_t i;
int rv;
for (i = 0; tests[i].rv >= 0; ++i) {
if (tests[i].a == tests[i].b)
printf("approx failed: malformed test %d: %15.15g == %15.15g\n",
(int)i, tests[i].a, tests[i].b);
rv = approx(tests[i].a, tests[i].b);
if (rv != tests[i].rv)
printf("approx failed: %15.15g %15.15g -> %d should be %d (test %d) # %15.15g %15.15g\n",
tests[i].a, tests[i].b, rv, tests[i].rv, (int)i,
tests[i].a * CUTOFF, tests[i].b * CUTOFF);
}
}
int
main(void)
{
test_approx();
printf("math done\n");
return 0;
}
#endif