Floating-Point Error Bounds for Geometric Predicates

1. Overview

When evaluating geometric predicates in double-precision floating-point arithmetic, the computed result accumulates rounding error. If the true value of the determinant is close to zero, this error can flip the sign, yielding an incorrect predicate result. A static filter defines a threshold: if the absolute value of the computed determinant exceeds this threshold, the sign is reliable; otherwise, tangle falls back to exact arithmetic.

This document derives the error bound for both the orient2D and inCircle predicates, expressed in terms of \(M\), the maximum absolute value of the coordinate differences appearing in the computation. The derivation follows the standard forward error analysis framework (see Higham (2002)) applied to the specific evaluation order used in the code, and is consistent with the error coefficients derived by Shewchuk (1997).

2. Preliminaries

2.1 Floating-Point Error Model

Under IEEE 754 double precision with round-to-nearest, every basic arithmetic operation satisfies:

\[ \mathrm{fl}(a \oplus b) = (a \oplus b)(1 + \delta), \quad \left| \delta \right| \le \varepsilon = 2^{-53} = 1.1102230246251565 \times 10^{-16} \]
where \(\varepsilon\) is the unit roundoff (machine epsilon). For a chain of \(n\) such operations, the accumulated relative error is bounded by \(\gamma_n = n\varepsilon/(1 - n\varepsilon) \approx n\varepsilon\) for small \(n\varepsilon\).

2.2 The Permanent

For a sum of products \(\sum t_i\), the standard forward error bound is:

\[ \left|\mathrm{computed}\left(\sum t_i\right) - \mathrm{exact}\left(\sum t_i\right)\right| \le \gamma_n \times \sum |t_i| \]
The quantity \(\sum \left| t_i \right| \) is the permanent of the determinant expression: the same polynomial with all subtractions replaced by additions. Shewchuk's predicates.c computes the actual permanent from the input values; in our simplified approach, tangle bounds the permanent in terms of \(M = \max( \left| \text{coordinate differences} \right| )\).

3. orient2D

3.1 The Computation

Given three points a, b, c, tangle computes:

Ax = b.x - a.x,  Ay = b.y - a.y
Bx = c.x - a.x,  By = c.y - a.y
det = Ax * By - Bx * Ay

3.2 Bounding the Permanent

3.3 Error Coefficient

Shewchuk's careful operation count for the orient2D evaluation yields the error bound coefficient (from predicates.c):

ccwerrboundA = (3.0 + 16.0 * epsilon) * epsilon ≈ 3ε

This accounts for the subtractions producing Ax, Ay, Bx, By (each introducing one rounding error), the two multiplications, and the final subtraction.

3.4 Error Bound in Terms of M

Combining the error coefficient with the permanent bound:

\[ \left| \mathrm{err} \right| \le 3\varepsilon \times 2M^2 = 6\varepsilon \times M^2 \]
Numerically:

\[ 6\varepsilon = 6 \times 2^{-53} = 6.661338147750939 \times 10^{-16} \]
Therefore the safe threshold is approximately \(6.662 \times 10^{-16} \times M^2\).

3.5 Chosen Threshold

Tangle uses the theoretical bound directly as the threshold:

if(std::abs(det) > 6.662e-16 * M * M) return det;

No additional safety margin is needed. The bound \(P \le 2M^2\) is itself conservative: it requires all four coordinate differences to simultaneously have magnitude \(M\). On real inputs the actual permanent is typically much smaller, so the effective safety margin is substantial even at 1× the theoretical bound.

4. inCircle

4.1 The Computation

Given four points a, b, c, d, tangle computes:

Ax = a.x - d.x,  Ay = a.y - d.y
Bx = b.x - d.x,  By = b.y - d.y
Cx = c.x - d.x,  Cy = c.y - d.y
As = Ax*Ax + Ay*Ay
Bs = Bx*Bx + By*By
Cs = Cx*Cx + Cy*Cy
det = Ax*(By*Cs - Cy*Bs)
	- Ay*(Bx*Cs - Cx*Bs)
	+ As*(Bx*Cy - Cx*By)

4.2 The Determinant and Its Permanent

4.3 Bounding the Permanent in Terms of M

Since each coordinate difference is bounded by \(M\), the squared-norm terms satisfy \(As, Bs, Cs \le 2M^2\). Bounding each group:

Therefore:

\[ P \le 12M^4 \]

4.4 Error Coefficient

Shewchuk's operation count for the inCircle evaluation yields (from predicates.c):

iccerrboundA = (10.0 + 96.0 * epsilon) * epsilon ≈ 10ε

This accounts for the six subtractions producing coordinate differences, the multiplications forming the squared norms, the cofactor multiplications, and the final summation.

4.5 Error Bound in Terms of M

Combining the error coefficient with the permanent bound:

\[ \left| \mathrm{err} \right| \le 10\varepsilon \times 12M^4 = 120\varepsilon \times M^4 \]

Numerically:

\[ 120\varepsilon = 120 \times 2^{-53} = 1.3322676295501878 \times 10^{-14} \]

Therefore the safe threshold is approximately \(1.333 \times 10^{-14} \times M^4\).

4.6 Chosen Threshold

Tangle uses the theoretical bound directly as the threshold:

if(std::abs(det) > 1.333e-14 * M * M * M * M) return det;

As with orient2D, no additional safety margin is needed. The bound \(P \le 12M^4\) requires all six coordinate differences to simultaneously have magnitude \(M\), which essentially never occurs with real geometric inputs. The built-in conservatism of the M-based bound provides a substantial effective safety margin.

5. Comparison with Shewchuk's Approach

Shewchuk's predicates.c uses a different filter. Rather than bounding the permanent by a power of \(M\), it computes the actual permanent from the input coordinate differences and multiplies it by an error coefficient, ultimately creating a tighter but more expensive bound. In Shewchuk's implementation, the fallback is not directly to full exact arithmetic but to progressively more precise adaptive stages to minimize the time needed to evaluate the fallback.

Tangle's approach trades some unnecessary fallbacks to exact arithmetic for a less expensive bound. The built-in conservatism of the \(M\)-based bound -- typically an order of magnitude or more on real inputs -- makes this a robust choice for a two-stage (fast floating-point / exact fallback) implementation, even without an explicit additional safety factor.

6. Summary

Predicate	Theoretical Bound	Chosen Threshold
orient2D	\(6\varepsilon \times M^2 = 6.661338147750939 \times 10^{-16} \times M^2\)	\(6.662 \times 10^{-16} \times M^2\)
inCircle	\(120\varepsilon \times M^4 = 1.3322676295501878 \times 10^{-14} \times M^4\)	\(1.333 \times 10^{-14} \times M^4\)

In both cases, \(M\) is the maximum absolute value of the coordinate differences (not the raw coordinates). No explicit safety margin is applied because the use of \(M\) to bound the permanent already provides substantial conservatism on real inputs.

References

N. J. Higham, Accuracy and Stability of Numerical Algorithms, 2nd ed. , SIAM 2002. DOI: 10.1137/1.9780898718027

J. R. Shewchuk, "Adaptive Precision Floating-Point Arithmetic and Fast Robust Geometric Predicates", Carnegie Mellon University, Technical Report CMU-CS-96-140R, 01 Oct 1997.

Finite Element Method Magnetics : PredicateBounds