A Robust Predicate Library for Delaunay Triangulation on Mixed-Precision Hardware

1. Introduction

A triangulation of a point set is Delaunay when no point lies inside the circumscribing circle, or sphere, of any simplex, the empty sphere condition introduced by Delaunay (1934). Algorithms that construct such triangulations, whether incremental, divide-and-conquer or flip-based, take all their geometric decisions through two predicates (Guibas and Stolfi, 1985). The orientation test asks on which side of a line or plane a point lies. The in-circle or in-sphere test asks whether a point lies inside, on or outside the circle or sphere through the others. Both are the sign of a small determinant of the input coordinates.

Evaluated in floating-point arithmetic, such a determinant can come out with the wrong sign when its true value is close to zero. The consequences are out of proportion to the size of the error. Kettner et al. (2008) exhibit small inputs on which a textbook incremental algorithm, run with a naive floating-point orientation test, produces a wrong result or fails to terminate. The accepted cure is the exact geometric computation paradigm: the predicates must return the sign of the exact determinant, whatever arithmetic is used internally (Yap, 1997).

The best known implementation of that paradigm is the set of adaptive predicates of Shewchuk (1997), written for processors on which double precision is the native format. Triangulation codes now also run on graphics processors and other accelerators, where the fast format is narrower and the arithmetic environment is less uniform. This article proposes a design for a predicate library suited to such mixed-precision hardware. It describes a method and the reasoning behind it. No implementation is evaluated, and no performance claims are made.

2. Background

The design combines three bodies of established work: the analysis of rounding error in the predicates, the filtering techniques that avoid exact arithmetic in easy cases, and the error-free transformations from which exact arithmetic can be built out of ordinary floating-point operations.

2.1. Rounding and the IEEE formats

The IEEE 754 standard defines binary formats with significands of 11, 24 and 53 bits, known as binary16, binary32 and binary64, and requires that each basic operation return the correctly rounded result (IEEE, 2019). Correct rounding is what makes rigorous error bounds possible: the computed result of one operation differs from the exact one by at most half a unit in the last place, and bounds for a whole expression follow by the standard techniques described by Higham (2002). Goldberg (1991) remains a clear account of where the guarantees end, for example when a compiler keeps intermediate values in a wider register or reorders operations.

2.2. Filters

A filter evaluates the determinant in fast arithmetic together with a bound on the accumulated error, and accepts the sign when the magnitude of the computed value exceeds the bound. Static filters derive the bound before run time from assumed limits on the coordinates, as in the integer-coordinate analysis of Fortune and Van Wyk (1996). Semi-static filters scale a precomputed constant by the magnitude of the actual operands. Dynamic filters compute the bound alongside the value, usually with interval arithmetic (Brönnimann et al., 2001). Devillers and Pion (2003) show how these are layered for Delaunay predicates so that almost all calls are settled by the cheapest test. The general principle, cheap tests first and costly ones only on failure, is often called cascading.

2.3. Error-free transformations and expansions

The sum of two floating-point numbers can be written exactly as a floating-point result plus a floating-point error term, and the error term can be computed with a few further additions. One such algorithm is due to Knuth (Knuth, 1997), and a faster one that requires the operands to be ordered by magnitude, together with a splitting technique that gives the analogous result for products, is due to Dekker (1971). Chaining these transformations yields expansions, sums of floating-point components that represent a number exactly, an idea developed by Priest (1991) and made adaptive by Shewchuk (1997). Where a fused multiply-add is available, the error of a product costs a single extra operation (Muller et al., 2018).

3. What Changes on Mixed-Precision Hardware

Four assumptions behind the classical predicates need to be revisited. The first is the working format. Coordinates on an accelerator are commonly stored in binary32, sometimes in binary16 for compressed geometry, and double-precision operations may be markedly slower or absent. The second is the fused multiply-add. It changes the rounding of expressions that a compiler contracts silently, which can invalidate an error bound derived for separate multiplication and addition. The third is the treatment of subnormal numbers, which some devices flush to zero by default, breaking the exactness of error-free transformations near the underflow threshold. The fourth is control flow: data-parallel execution penalizes branches, so that a cascade in which different inputs leave at different stages has to be organized differently from a scalar routine.

Mixed precision also offers an opportunity that uniform double precision lacks. The product of two binary32 numbers has at most 48 significant bits and is therefore exactly representable in binary64, barring overflow and underflow. Narrow inputs evaluated in a wider format thus lose nothing in their first multiplications. This is the principle on which mixed-precision methods in numerical linear algebra are built (Higham and Mary, 2022), and it can be used in geometric predicates in the same way.

4. The Proposed Design

The library is organized as an input contract followed by four evaluation stages and a tie-breaking rule. Figure 1 shows the arrangement, and Table 1 states for each stage the arithmetic it uses, the condition under which it returns a certified sign, and the hardware property on which its correctness depends. A call descends to the next stage only when the current one cannot certify the sign.

4.1. Input contract

Each point set passed to the library carries a descriptor stating the format of its coordinates and, optionally, a bounding box and a grid spacing. The descriptor selects the constants used by the filters. If the coordinates are known to be multiples of a grid spacing within a bounded box, they are scaled integers of known bit length, and the bit length required by each determinant is known in advance: roughly twice the coordinate length for the planar orientation test and four times for the in-circle test, plus a few guard bits. Inputs without such guarantees are treated as arbitrary values of their format. Not-a-number and infinite coordinates are rejected at this point and never inside a predicate.

4.2. Stage A: semi-static filter in the native format

The determinant is evaluated in the fast narrow format by branch-free code, after translating the points so that one of them is the origin. In parallel the same expression is evaluated on absolute values, and the result is multiplied by a constant derived from a forward error analysis of the expression in that format. The sign is accepted if the magnitude of the determinant exceeds this bound. The stage is written so that it can be applied to a whole batch of simplices at once, with the undecided ones compacted into a list for the next stage. Contraction into fused multiply-adds must be disabled for this code, or else the constant must be derived for the contracted form.

4.3. Stage B: filter in a wider format

Undecided calls are re-evaluated in binary64. For binary32 inputs the coordinate differences are formed in binary64, where they are exact whenever the operands are close enough in exponent for the result to fit in 53 bits, a condition the bounding box and grid in the input contract can guarantee. The first products are then nearly exact, and the error bound shrinks by many orders of magnitude compared with Stage A. On devices without binary64 the same stage is implemented with double-word arithmetic, in which each quantity is held as an unevaluated sum of two binary32 numbers and manipulated by error-free transformations (Muller et al., 2018), and the bound is adjusted accordingly.

4.4. Stage C: adaptive expansion arithmetic

The few calls that remain are close to degenerate. They are evaluated with expansions in the widest native format, following the adaptive scheme of Shewchuk (1997): the determinant is computed in successively more accurate approximations, each reusing the components already obtained, and evaluation stops as soon as the approximation exceeds its remaining error bound. Products use the fused multiply-add when the device provides a correctly rounded one and Dekker's splitting otherwise. The stage requires round-to-nearest and gradual underflow, and the library should verify both when it is initialized instead of assuming them.

4.5. Stage D and tie-breaking

If the adaptive stage runs to completion without certifying a nonzero sign, the full expansion has been computed and its sign is the exact answer. For grid inputs an alternative is multiword integer arithmetic of the length fixed by the input contract. A result of exactly zero means that the points are collinear, cocircular or the higher-dimensional equivalent. The library then applies a symbolic perturbation in the manner of the simulation of simplicity of Edelsbrunner and Mücke (1990), resolving the tie by comparing the global indices of the points, so that the triangulation algorithm never sees a degenerate case and the output is uniquely defined. Lévy (2016) describes a toolkit that generates filters, expansions and perturbation code from one formula, and the same generative approach is recommended here, since hand-written variants of one determinant are a known source of errors.

5. Verification Plan

A predicate library is useful only if its claim of exactness can be checked. The proposal is to test at three levels. Each error-free transformation is tested for exactness against multiword integer arithmetic on operands chosen to stress carries, cancellation and the underflow boundary. Each predicate is then compared with an exact rational evaluation on inputs that are deliberately near-degenerate: points on a common line or circle perturbed by one unit in the last place, in the style of the failure-inducing examples of Kettner et al. (2008). Finally the device itself is probed at start-up for rounding mode, subnormal support, fused multiply-add behaviour and contraction, and the library refuses to run, or selects a slower safe path, when a required property is absent.

6. Limitations and What Remains to Be Shown

The design set out here is a proposal and has not been validated. Its correctness rests on error bounds that have to be derived, and preferably machine-checked, for each predicate in each format, and a mistake in a single constant would defeat the purpose. Its efficiency rests on the expectation that the narrow filter settles the large majority of calls. That is plausible for well-distributed inputs, since the narrow format leaves fewer bits of margin but most configurations are far from degenerate, yet for inputs with many nearly cocircular points, such as samples of spheres or regular grids, the proportion of calls reaching the later stages may be high. Only measurement on real devices can settle this.

Two further limits should be stated. The library guarantees correct predicates on the coordinates as stored, and says nothing about the error committed when those coordinates were rounded from other data. And it covers predicates only. Constructions such as circumcentres, which Voronoi diagrams and mesh refinement require, cannot be made exact in fixed-size floating-point formats and need separate treatment.

7. Conclusion

Exact predicates for Delaunay triangulation need not be confined to double-precision processors. A cascade that starts in the narrow native format, widens where widening makes products exact, and falls back on expansion arithmetic and symbolic perturbation can in principle deliver exact signs on hardware whose fast arithmetic is single precision. The components are all established. What this article adds is their arrangement around an explicit input contract and an explicit list of the hardware properties each stage depends on. Implementation, formal verification of the error bounds and measurement on current accelerators are the necessary next steps.

Conflict of interest

To be declared by the author(s) before publication.

Funding

To be declared by the author(s) before publication.