Voronoi Robust FPT

The Voronoi robust floating-point types are the set of classes and tools that allow to estimate relative error of the arithmetic expressions. It is assumed that the other Boost libraries may find this unit functionality extremely useful, as it can be used to implement robust and efficient arithmetic predicates or functors that compute values within known relative error.

Robust Fpt Type

The robust fpt type (robust floating-point type) - represents the IEEE-754 floating-point type wrapper that also contains information about the relative error of the underlying value. The implementation overloads 5 standard operations: +, -, *, /, sqrt and apart from the evaluating value of the expression also computes its relative error. Let's consider two values A and B; C - rounding error, re(X) - relative error of the X expression, then following rules apply:

re(A+B) <= max(re(A), re(B)) + C, if A * B >= 0;
re(A-B) <= (B * re(A) + A * re(B)) / |A - B| + C, if A * B < 0;
re(A*B) <= re(A) + re(B) + C;
re(A/B) <= re(A) + re(B) + C;
re(sqrt(A)) <= re(A) * 0.5 + C;

The constant C is equal to the rounding error, which for the above set of arithmetic operations in the IEEE-754 floating-point implementation should be equal to 1 machine epsilon.

Robust Difference Type

The robust difference type - represents expression wrapper that holds the positive and negative partial sums of the expression in a separate values in order to avoid the cancellation errors before evaluating the final difference. Following arithmetic operators are overloaded for the robust difference type: +, -, *, / (division operator is not overloaded for the case were both arguments have robust difference type).
Looking at the relative error formulas above one may notice a few facts:
1) all of the formulas evaluate upper bound of the relative error, the actual value could be a lot smaller;
2) relative error estimate for the expression depends on the order operations are evaluated;
3) relative error of  the difference of two positive numbers may be extremely large in case their values are close to each other (this is also known as the cancellation error).
To explain this a bit, consider the following expression (~ - stands for almost equal, << - many times larger than):

A - B + C, where A ~ B and C << B;

Computing the relative error of this expression from left to right will produce extremely large relative error:

re(A-B+C) = max(re(A-B), re(C)) = re(A-B) = (B * re(A) + A * re(B)) / 0 = INF;

While doing this from right to left will keep the relative error value small:

re(A-B+C) = re(C-B+A) = max(re(C-B), re(A)) = max(re(A), re(B));

While both estimates are valid (they define upper bound of the relative error), of course the second one is preferable.
Here is the place where robust difference type comes useful. Basically it splits expression onto positive and negative partial sums and evaluates the difference only when the result is required. And did I mention that positive and negative values might be of the robust fpt type, that's why the relative error is always known for the expression result.

Robust Sqrt Expression Structure

The robust square root expression structure allows to compute the result of the expression that contains square roots within defined relative error. As an example, consider the following expression:

A * sqrt(a) - B * sqrt(b), A * B > 0, a >= 0, b >= 0;

Computing this expressions directly may apply huge cancellation error, however it may be transformed to the next equivalent expression:

(A * A * a - B * B * b) / (A * sqrt(a) + B * sqrt(b));


The numerator and denominator of this expression could be computed directly as those won't lead to the cancellation errors.

In general case the robust sqrt expression structure allows to evaluate the following set of expressions:

sum(A[i] * sqrt(a[i]), i = 1 .. N), N <= 4;


This appears to be enough for the Boost.Polygon Voronoi.
 
Copyright: Copyright © Andrii Sydorchuk 2010-2012
License: Distributed under the Boost Software License, Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)