Is my epsilon usage correct for comparing floating-point numbers in circle overlap calculations?

Question

I'm working on a C program that calculates the overlap area between two circles. The program needs to handle various cases: circles coinciding, one inside another, touching internally/externally, separated, or intersecting. I'm using epsilon for floating-point comparisons to handle precision issues, but I'm not sure if I'm using it correctly. Specifically:

1. Is my epsilon definition appropriate: const double eps = 1e6 * DBL_EPSILON;?
2. Am I using it correctly in comparisons like fabs(d) < eps and d < fabs(r1 - r2) - eps?

Problem description:

Input: Two circles with centers (x1, y1), (x2, y2) and r1, r2.

• Calculate distance d between centers.
• Determine relationship between circles and calculate overlap area

My concerns:

1. Is 1e6 * DBL_EPSILON too large or too small for my use case?

2. Should I be using relative epsilon instead of absolute epsilon for some comparisons?

3. Are there edge cases where my epsilon comparisons might fail?

Any guidance on proper epsilon usage for geometric calculations would be appreciated!

#include <math.h>
#include <stdio.h>
#include <float.h>

const double eps = 1e6 * DBL_EPSILON;

double middleSum(double x1, double x2, double y1, double y2) {
    return sqrt(pow(x2 - x1, 2) + pow(y2 - y1, 2));
}

double getOverlap(double x1, double y1, double r1, double x2, double y2, double r2) {
    double d = middleSum(x1, x2, y1, y2);
    double R1 = r1;
    double R2 = r2;

    if (fabs(d) < eps && fabs(r1 - r2) < eps) {
        double minR = (R1 < R2 ? R1 : R2);
        return PI * minR * minR;
    }

   return R1 * R1 * acos((d * d + R1 * R1 - R2 * R2) / (2 * d * R1)) + R2 * R2 * acos((d * d + R2 * R2 - R1 * R1) / (2 * d * R2)) - 0.5 * sqrt((-d + R1 + R2) * (d + R1 - R2) * (d - R1 + R2) * (d + R1 + R2));

}

int main(void) {
    // ... input reading ...
    
    double d = middleSum(x1, x2, y1, y2);

    // Check if circles coincide
    if (fabs(d) < eps && fabs(r1 - r2) < eps) {
        printf("Circles coincide\n");
    }
    // Check if one circle is inside another
    else if (d < fabs(r1 - r2) - eps) {
        printf("One circle inside another\n");
    }
    // Check for internal tangency
    else if (fabs(d - fabs(r1 - r2)) < eps) {
        printf("Internal tangency\n");
    }
    // Check for external tangency
    else if (fabs(d - (r1 + r2)) < eps) {
        printf("External tangency\n");
    }
    // Check if circles are separated
    else if (d > r1 + r2 + eps) {
        printf("Circles separated\n");
    }
    // Circles intersect
    else {
        printf("Circles intersect\n");
    }
    
    return 0;
}

Answer 1

1. Is 1e6 * DBL_EPSILON too large or too small for my use case?
   IMO, it is better as 1.0 * DBL_EPSILON or it is simply irrelevant.

2. Should I be using relative epsilon instead of absolute epsilon for some comparisons?
   It depends on the functional goals, yet usually with floating point relative epsilon is better than absolute epsilon. An even better approach involves using the difference in the total ordering. I. e., if you sequenced all double in a list, how far apart are the indexes of a and b in the list, that is, some integer.

3. Are there edge cases where my epsilon comparisons might fail?
   Yes, It is not clear that when 2 close values are compared, which path should be taken, if not a 3rd path.

---

It's floating point, not fixed point

Compares like fabs(d) < eps take the float out of floating point.

Say you have 2 circles, each defined by by 3 double: center_x, center_y, radius and code then reports: "circles coinciding, one inside another, touching internally/externally, separated, or intersecting".

Scale each of the 6 double by 2¹⁰, does one get the same report?

Scale each of the 6 double by 2^-10, does one get the same report?

Barring computational over/underflow, I'd expect the same reports. Yet fabs(d) < eps will always return true for tiny circles and too often false for large ones.

If an epsilon is needed here, using relative comparison is better.

---

Is epsilon truly needed?

This has not yet been demonstrated.

Consider x1 and the subsequent double: x2 = nextafter(x1, INFINITY);. Aside from the edge case, is x1 == x2 ever equal? Likewise for the circle compare, unless we have some +/- error information about the circle's input parameters, they should be considered exact and it is only our calculation's round-off errors should be considered.

// ... input reading ...

This unposted code deserves posting to understand clearly how input is taken and converted to double.

Text input like "12.3" could be considered to have an error as a double cannot encode that value exactly. Input like "1.25" or "0x1.23abcp0" could be considered exact when converted to a double.

If input is exact, the case of circle coincidence, it is simple enough to compare the 2 circles exactly as there is no computational error.

If input is not certainly exact, it should be considered to have an error of about DBL_EPSILON * the_value or even better as +/- 0.5 ULP.

Reduce errors

We can typically reduce some calculation errors as in:

double middleSum(double x1, double x2, double y1, double y2) {
    return sqrt(pow(x2 - x1, 2) + pow(y2 - y1, 2));
}
double middleSum_alt(double x1, double x2, double y1, double y2) {
    return hypot(x2 - x1, y2 - y1);
}

When higher precision floating point types are available, code could use those for intermediate calculations.

Account for computational errors

When a, b may have computational errors:

Rather the a 2-way branch of a < b, code needs a 3-way branch of

// Concept code
if (less_than_with_margin(a, b, a_error, b_error)) {
  // Do the a < b thing
else if (greater_than_or_equal_with_margin(a, b, a_error, b_error)) {
  // Do the !(a < b) thing
else
  // A 3rd path as a and b are so close to each other yet have errors such that 
  // determining a < b is not discernable.

less_than_with_margin() should return true only when a + a_error < b - b_error.

Minor: Unneeded code

fabs(d) can be replaced with d as the result of d = middleSum() is not less than 0.0.

---

Bottom line

Unless the problem has more info about the 1) nature of the input and 2) how to handle the grey zone cases, remove epsilon from all the categorizations.

If this is some academic assignment, I fear it is simply too under-defined for a good answer.