Exact

Exact numbers and exact predicates

What is the problem ?

The geometric algorithms in Geogram, such as 2D triangulations, 3D triangulations, mesh intersection and mesh CSG depend on on more "elementary questions", such as whether a point is above or below a plane (in a certain sense).

Such "elementary questions", called predicates, are functions that take as an argument a (small) number of points (or simple geometric objects) and that returns a set of discrete values.

For instance, consider four points ${\bf p}_1,{\bf p}_2,{\bf p}_3,{\bf p}_4$. One may want to know the position of ${\bf p}_4$ relative to the supporting plane of ${\bf p}_1,{\bf p}_2,{\bf p}_3$, that can be one of ABOVE, BELOW, ON_PLANE.

These predicates are the "nevralgic" point of mesh intersection methods: if at one moment the algorithm "thinks" that ${\bf p}_4$ is above the supporting plane of ${\bf p}_1,{\bf p}_2,{\bf p}_3$, it is important that at another moment the predicate does not say that ${\bf p}_4$ is below the same plane. How could this happen ?

In fact, these predicates can be expressed as the sign of a polynomial in the coordinates of the points, and due to the limited precision of floating point numbers, the output of the predicate can be different from the exact mathematical result, especially around zero, and it can depend on the order of the points. Try this: write the following little C program:

#include <stdio.h>

void test_add3(double x, double y, double z) {
   printf("%f + %f + %f = %f\n",x,y,z,x+y+z);
}

int main() {
   double x1 = 1e30;
   double x2 = -1e30;
   double x3 = 1e-6;
   test_add3(x1,x2,x3);
   test_add3(x1,x3,x2);
   return 0;
}

The first test computes $x_1+x_2+x_3$ and the second one $x_1+x_3+x_2$, mathematically they shoud return the same result, but we got limited precision: in the first expression, everything goes well, we first compute $x_1$ + $x_2$, the result is exactly zero, then add $x_3$, and the final result is $1e-6$. In the second expression, it is another story, when ones computes $x_1+x_3$, that is, $1e30$ + $1e-6$, double precision numbers cannot store all the digits, then the result is truncated (and in this case, it means ignoring $1e-6$. Finally, when we subtract $x_3$, we get exactly zero (which is wrong !).

Why should we care ? If you think of what may happen in the middle of a geometric algorithm, since the result of an operation depends on the order of the operands, you may obtain a different result when you ask for the position of ${\bf p}_4$ relative to the supporting plane of ${\bf p}_1, {\bf p}_2, {\bf p}_3$ or when you ask for the position of ${\bf p}_4$ relative to ${\bf p}_3, {\bf p}_2, {\bf p}_1$ ! It can have catastrophic consequences, such as generating an incorrect mesh.

So what can we do ?

The idea is to implement a C++ class expansion_nt that represents numbers exactly. How is it possible in a computer ? First thing, we will implement a very limited set of operations:

• addition: expansion_nt + expansion_nt $rightarrow$ expansion_nt
• subtraction: expansion_nt - expansion_nt $rightarrow$ expansion_nt
• product: expansion_nt * expansion_nt $rightarrow$ expansion_nt
• sign: expansion_nt $rightarrow$ NEGATIVE,ZERO,POSITIVE

See also...

• Programmer's reference: Exact numbers

• Want to learn more about how it works ? See this article and this one

The optional geogramplus expansion package

To gain more speed and more robustness in the extreme cases, read about GeogramPlus.