HAL Id: inria-00071232
https://hal.inria.fr/inria-00071232v2
Submitted on 4 Dec 2008
HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entic research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.
L’archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diusion de documents
scientiques de niveau recherche, publiés ou non,
émanant des établissements d’enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.
Formally Certied Floating-Point Filters For
Homogeneous Geometric Predicates
Guillaume Melquiond, Sylvain Pion
To cite this version:
Guillaume Melquiond, Sylvain Pion. Formally Certied Floating-Point Filters For Homogeneous Geo-
metric Predicates. RAIRO - Theoretical Informatics and Applications (RAIRO: ITA), EDP Sciences,
2007, 41, pp.57-69. �10.1051/ita:2007005�. �inria-00071232v2�
Theoretical Informatics and Applications Will be set by the publisher
Informatique Th´eorique et Applications
FORMALLY CERTIFIED FLOATING-POINT FILTERS FOR
HOMOGENEOUS GEOMETRIC PREDICATES
Guillaume Melquiond
1
and Sylvain Pion
2
Abstract. Floating-point arithmetic provides a fast but inexact way
of computing geometric predicates. In order for these predicates to be
exact, it is important to rule out all the numerical situations where
floating-point computations could lead to wrong results. Taking into
account all the potential problems is a tedious work to do by hand. We
study in this paper a floating-point implementation of a filter for the
orientation-2 predicate, and how a formal and partially automatized
verification of this algorithm avoided many pitfalls. The presented
metho d is not limited to this particular predicate, it can easily be used
to produce correct semi-static floating-point filters for other geometric
predicates.
1991 Mathematics Subject Classification. 65G50,68Q60,65D18.
1. INTRODUCTION
Computational geometry algorithms, such as convex hull computations, Delau-
nay triangulations and arrangements computations, are notoriously sensitive to
numerical instability. The most important characteristic of these algorithms is
that they mix numerical com putations, with combinatorial ones. Their input is
for example a set of points in the plane or in space, given by their coordinates,
and their output contains a combinatorial part such as a graph: a simple linear
Keywords and phrases: Geometric predicates, semi-static filters, formal proofs, floating-point
1
Laboratoire de l’Informatique du Parall´elisme
UMR 5668 CNRS, ENS Lyon, INRIA, UCBL
46 all´ee d’Italie, 69364 Lyon C ede x 07, Franc e
E-mail: Guillaume.Melquiond@ens-lyon.fr
2
INRIA Sophia Antipolis,
2004 route des Lucioles, BP 93, 06902 Sophia Antipolis, France
E-mail: Sylvain.Pion@sophia.inria.fr
c
EDP Sciences 1999
2 TITLE WILL BE SET BY THE PUBLISHER
sequence for the 2-dimensional convex hull, or more complex ones in the case of
triangulations... Many such algorithms are gathered in the CGAL
1
library [3].
In order to derive a combinatorial, discrete, structure from a set of numerical
inputs, a set of functions are used, called geometric predicates, which evaluate
the relative p os itions of a few geometric objects, such as the orientation of three
points in the plane. From the implementation point of view, using floating-point
approximate arithmetic to evaluate these predicates has shown to be the source of
many non-robustness problems, because the incorrect geometry of these approxi-
mate predicates violates basic geometric theorems on which the algorithms rely [8].
One of the most appreciated solutions to this problem, due to its generality, is to
render these predicates exact, thus following the Exact Geometric Computation
paradigm [14].
The use of exact multi-precision arithmetic instead of floating-point provides
the necessary exactness, but it is not usable naively in practice, due to considerable
efficiency loss. Therefore an additional step has been developed in what is called
arithmetic filters, which are a way to filter out easy cases of the predicates using
certified floating-point arithmetic, calling the slow exact arithmetic as a last resort,
far less often. There are different kinds of filters, varying in efficiency and precision.
Static filters have been first developed [6], and rely on static forward analysis of
error propagation. They tend to be restricted in the input, but are fast. Many
variants exist [2, 5, 12]. Dynamic filters are easier to use and more general, but
slower [1].
Static filters variants use floating-point error propagation, which are risky and
error prone when done by hand. Some more or less automatic tools have been
developed to pro duce their code [2,7,9,11], but it remains an error-prone task due
to the complex nature of floating-point arithmetic.
In this paper, we detail the formal proof of a particular kind of static filters,
which is implemented in CGAL [3]. We focus on the 2-dimensional orientation
of points as predicate, since it is one of the most critical. We then show that
the techniques apply similarly to many more important other predicates. And we
conclude with a comparison to other existing methods.
2. ORIENTATION-2 PREDICATE
This predicate is one of the most often encountered geometric predicates. Given
three points p, q, and r in the plane, it answers if they are collinear, or if they are
clockwise or counter-clockwise oriented. If the Cartesian coordinates of the three
points are known, the answer is given by the sign of a 3×3-determinant of the
coordinates, or the sign of a 2×2-determinant of the vectors.
orient
2
(p, q, r) =
p
x
q
x
r
x
p
y
q
y
r
y
1 1 1
=
q
x
− p
x
r
x
− p
x
q
y
− p
y
r
y
− p
y
1
http://www.cgal.org/
TITLE WILL BE SET BY THE PUBLISHER 3
2.1. Naive implementation
Algorithm 1 Floating-point orientation-2 naive implementation
double pqx = qx - px, pqy = qy - py;
double prx = rx - px, pry = ry - py;
double det = pqx * pry - pqy * prx;
if (det > 0) return POSITIVE;
if (det < 0) return NEGATIVE;
return ZERO;
Algorithm 1 shows a naive implementation of the orientation-2 predicate. It
first computes floating-point approximations of the vectors pq and pr. Then it
simply computes an approximation of the 2×2-determinant and returns its sign.
If the computations were done on real numbers, this implementation would
lead to the correct result. Unfortunately floating-point numbers suffer from both
a limited precision and a limited range. These limitations can lead to a det value
with a different sign.
Let us consider the following example. When the result of a floating-point
computation is outside the range of representable numbers, it is replaced by an
infinite value. Let M be the biggest power-of-2 number that is still representable.
Let p, q, and r be
(p, q, r) =
−M M −7M/8
1 3 −17/16
Computing pqx will produce an unrepresentable value 2M. As a consequence
the variable will contain +∞. None of the other variables nor pqy * prx will
contain an infinite value. Hence the infinite value contained in pqx will propagate
till the final variable det and the predicate will answer POSITIVE. Unfortunately
the real value of the determinant is a negative number: −M/8.
This example shows that the limited range of floating-point numbers can lead to
the wrong sign. The limited precision can also lead to a wrong answer, especially
when the points are almost aligned. A small rounding error will spoil each floating-
point computations and the computed values may slowly drift away from the real
values.
2.2. Floating-point considerations
The IEEE-754 standard [13] covers both the data formats and the precise be-
havior of the arithmetic operators. In particular, it describes how rounding and
limited precision affect the operators.
A floating-point operation ˜c ← ˜a ⊕
˜
b can theoretically be decomposed in two
steps: the first one computes the exact real result c = ˜a +
˜
b and the second one
rounds this exact result c to the nearest floating-point number ˜c. The distance
between the exact result and the computed value is bounded by
0
· |c| and by
4 TITLE WILL BE SET BY THE PUBLISHER
0
· |˜c|, if the result is in the range of normal numbers, and η
0
otherwise. If |c| is
bounded by r, and if
0
r ≥ η
0
, then the computational error is bounded by
0
r.
Both
0
and η
0
depend on the floating-point format chosen to implement the
algorithm. In case of the IEEE-754 double precision floating-point arithmetic with
rounding to nearest, these constants are
0
= 2
−53
and η
0
= 2
−1075
.
Unfortunately, the standard does not prevent a processor to use an excess pre-
cision when computing in its default mode of operations. It does not describe how
programming languages and their compilers should comply with it either. As a
consequence, a double rounding phenomenon can appear. A value will first be
computed and rounded in the default extended precision of the processor. It will
then be stored in memory with a reduced precision. This storage induces a second
rounding.
The total rounding error of these two operations can then exceed the rounding
error that would have happened if a single rounding had been done. For example,
this phenomenon can happen in a program relying on double precision floating-
point computations compiled by GCC for the x86 processor family. De pending
on whether the result of a floating-point operation is only stored in a processor
register or goes through the memory, the value will be rounded one or two times.
The floating-point registers use 64-bit wide mantissas while the double type causes
the compiler to store the values in memory with 53-bit wide mantissas. In order
to take the double rounding phenomenon into account, the constants have to be
changed acc ordingly:
0
= 2
−53
(1 + 2
−11
) and η
0
= 2
−1075
(1 + 2
−12
).
2.3. Robust implementation
Algorithm 2 Floating-point orientation-2 filter
double pqx = qx - px, pqy = qy - py;
double prx = rx - px, pry = ry - py;
double maxx = max(abs(pqx), abs(prx));
double maxy = max(abs(pqy), abs(pry));
double eps = 8.8872057372592758e-16 * maxx * maxy;
if (maxx > maxy) swap(maxx, maxy);
if (maxx < 1e-146) { // underflows?
if (maxx == 0) return ZERO;
} else if (maxy < 1e153) { // no overflow?
double det = pqx * pry - pqy * prx;
if (det > eps) return POSITIVE;
if (det < -eps) return NEGATIVE;
}
// fall back to a more precise, slower method
