Adaptive Precision FloatingPoint Arithmetic and Fast Robust Predicates for
Computational Geometry
Jonathan Richard Shewchuk
Computer Science Division
University of California at Berkeley
Berkeley, California 947201776
Created as part of the
Archimedes project (tools for parallel finite
element methods).
Supported in part by
NSF Grant CMS9318163 and an
NSERC 1967 Scholarship.
Many computational geometry applications use numerical tests known
as the orientation and incircle tests. The
orientation test determines whether a point lies to the left of, to the
right of, or on a line or plane defined by other points. The incircle
test determines whether a point lies inside, outside, or on a circle
defined by other points.
Each of these tests is performed by evaluating the sign of a determinant
(see the figure below).
The determinant is expressed in terms of the coordinates of the points.
If these coordinates are expressed as single or double precision floatingpoint
numbers, roundoff error may lead to an incorrect result when
the true determinant is near zero. In turn, this misinformation
can lead an application to fail or produce incorrect results.
One way to solve this problem is to use exact arithmetic. Unfortunately,
traditional libraries for arbitrary precision floatingpoint arithmetic
are quite slow, and can reduce the speed of an application by one or two
orders of magnitude.
To minimize this problem, I've produced algorithms and implementations
for performing the 2D and 3D orientation and incircle tests correctly and
reasonably quickly. From this page, you may access papers describing
the approach (short and long versions), and the code itself. My hope
is that working computational geometers will be able to simply plug
the code into their applications with little effort, and forget about
robustness problems (for those algorithms that only require
orientation and incircle tests).
My algorithms are faster than traditional arbitrary precision libraries
for two reasons. First, I use an unusual approach to exact arithmetic,
pioneered by Douglas Priest and sped by me, that has its greatest advantage
over traditional techniques when manipulating numbers of relatively small
complexity (a few hundred or
thousand bits). Second, the algorithms are adaptive in the sense
that they do only as much work as necessary to guarantee a correct result.
The sign of a small determinant can typically be determined quickly unless the
determinant is close to zero. See the papers for details.
Papers
The long version.

Jonathan Richard Shewchuk, Adaptive Precision FloatingPoint
Arithmetic and Fast Robust Geometric Predicates,
Discrete & Computational Geometry 18:305363, 1997.
Also Technical Report CMUCS96140, School of Computer Science,
Carnegie Mellon University, Pittsburgh, Pennsylvania, May 1996.
Abstract (with BibTeX citation),
PostScript (676k, 53 pages).
The short version.
Covers the most important points,
but leaves out all the proofs and many details.
Presents a slower version of the addition algorithm because
there wasn't room to explain the fastest one.

Robust Adaptive FloatingPoint Geometric
Predicates, Proceedings of the Twelfth Annual Symposium on
Computational Geometry, ACM, May 1996.
Abstract (with BibTeX citation),
PostScript (310k, 10 pages).
Software

Here's
C code for the 2D and 3D orientation and incircle tests,
and for arbitrary precision floatingpoint addition and multiplication.
Last revised May 18, 1996.
This code is in the public domain.
I'm planning to add subtraction routines when I have a day free,
but haven't gotten around to it yet.
Warning: This code will not work correctly on a processor that
uses extended precision internal floatingpoint registers, such as
the Intel 80486 and Pentium, unless the processor is configured to
round internally to double precision. See the
Predicates on PCs page for details.

Triangle is a 2D Delaunay triangulator and quality mesh generator that
uses the 2D orientation and incircle tests for robustness. Triangle
is copyrighted; it is not in the public domain. However, it is
freely available from the Triangle page.
Inspiration

Douglas M. Priest, Algorithms for
Arbitrary Precision Floating Point Arithmetic, Tenth Symposium on
Computer Arithmetic, pp. 132143, IEEE Computer Society Press, 1991.
Abstract (with BibTeX citation),
Copyright notice (please read),
PostScript (240k).

Steven Fortune and Christopher J. Van Wyk, Efficient Exact Arithmetic
for Computational Geometry, Proceedings of the Ninth Annual Symposium
on Computational Geometry, pp. 163172, Association for Computing Machinery,
May 1993.

Steven Fortune, Numerical Stability of Algorithms for 2D Delaunay
Triangulations, International Journal of Computational
Geometry & Applications 5(12):193213, MarchJune 1995.
Also check out the corresponding papers and software offered by
Avnaim, Boissonnat, Devillers, Preparata, and Yvinec. Their
technique may be faster than mine (I haven't timed theirs),
but is meant for points whose
coordinates are integers.
Jonathan Shewchuk