1 // Boost.Polygon library detail/voronoi_robust_fpt.hpp header file
2
3 // Copyright Andrii Sydorchuk 2010-2012.
4 // Distributed under the Boost Software License, Version 1.0.
5 // (See accompanying file LICENSE_1_0.txt or copy at
6 // http://www.boost.org/LICENSE_1_0.txt)
7
8 // See http://www.boost.org for updates, documentation, and revision history.
9
10 #ifndef BOOST_POLYGON_DETAIL_VORONOI_ROBUST_FPT
11 #define BOOST_POLYGON_DETAIL_VORONOI_ROBUST_FPT
12
13 #include <algorithm>
14 #include <cmath>
15
16 // Geometry predicates with floating-point variables usually require
17 // high-precision predicates to retrieve the correct result.
18 // Epsilon robust predicates give the result within some epsilon relative
19 // error, but are a lot faster than high-precision predicates.
20 // To make algorithm robust and efficient epsilon robust predicates are
21 // used at the first step. In case of the undefined result high-precision
22 // arithmetic is used to produce required robustness. This approach
23 // requires exact computation of epsilon intervals within which epsilon
24 // robust predicates have undefined value.
25 // There are two ways to measure an error of floating-point calculations:
26 // relative error and ULPs (units in the last place).
27 // Let EPS be machine epsilon, then next inequalities have place:
28 // 1 EPS <= 1 ULP <= 2 EPS (1), 0.5 ULP <= 1 EPS <= 1 ULP (2).
29 // ULPs are good for measuring rounding errors and comparing values.
30 // Relative errors are good for computation of general relative
31 // error of formulas or expressions. So to calculate epsilon
32 // interval within which epsilon robust predicates have undefined result
33 // next schema is used:
34 // 1) Compute rounding errors of initial variables using ULPs;
35 // 2) Transform ULPs to epsilons using upper bound of the (1);
36 // 3) Compute relative error of the formula using epsilon arithmetic;
37 // 4) Transform epsilon to ULPs using upper bound of the (2);
38 // In case two values are inside undefined ULP range use high-precision
39 // arithmetic to produce the correct result, else output the result.
40 // Look at almost_equal function to see how two floating-point variables
41 // are checked to fit in the ULP range.
42 // If A has relative error of r(A) and B has relative error of r(B) then:
43 // 1) r(A + B) <= max(r(A), r(B)), for A * B >= 0;
44 // 2) r(A - B) <= B*r(A)+A*r(B)/(A-B), for A * B >= 0;
45 // 2) r(A * B) <= r(A) + r(B);
46 // 3) r(A / B) <= r(A) + r(B);
47 // In addition rounding error should be added, that is always equal to
48 // 0.5 ULP or at most 1 epsilon. As you might see from the above formulas
49 // subtraction relative error may be extremely large, that's why
50 // epsilon robust comparator class is used to store floating point values
51 // and compute subtraction as the final step of the evaluation.
52 // For further information about relative errors and ULPs try this link:
53 // http://docs.sun.com/source/806-3568/ncg_goldberg.html
54
55 namespace boost {
56 namespace polygon {
57 namespace detail {
58
59 template <typename T>
get_sqrt(const T & that)60 T get_sqrt(const T& that) {
61 return (std::sqrt)(that);
62 }
63
64 template <typename T>
is_pos(const T & that)65 bool is_pos(const T& that) {
66 return that > 0;
67 }
68
69 template <typename T>
is_neg(const T & that)70 bool is_neg(const T& that) {
71 return that < 0;
72 }
73
74 template <typename T>
is_zero(const T & that)75 bool is_zero(const T& that) {
76 return that == 0;
77 }
78
79 template <typename _fpt>
80 class robust_fpt {
81 public:
82 typedef _fpt floating_point_type;
83 typedef _fpt relative_error_type;
84
85 // Rounding error is at most 1 EPS.
86 enum {
87 ROUNDING_ERROR = 1
88 };
89
robust_fpt()90 robust_fpt() : fpv_(0.0), re_(0.0) {}
robust_fpt(floating_point_type fpv)91 explicit robust_fpt(floating_point_type fpv) :
92 fpv_(fpv), re_(0.0) {}
robust_fpt(floating_point_type fpv,relative_error_type error)93 robust_fpt(floating_point_type fpv, relative_error_type error) :
94 fpv_(fpv), re_(error) {}
95
fpv() const96 floating_point_type fpv() const { return fpv_; }
re() const97 relative_error_type re() const { return re_; }
ulp() const98 relative_error_type ulp() const { return re_; }
99
has_pos_value() const100 bool has_pos_value() const {
101 return is_pos(fpv_);
102 }
103
has_neg_value() const104 bool has_neg_value() const {
105 return is_neg(fpv_);
106 }
107
has_zero_value() const108 bool has_zero_value() const {
109 return is_zero(fpv_);
110 }
111
operator -() const112 robust_fpt operator-() const {
113 return robust_fpt(-fpv_, re_);
114 }
115
operator +=(const robust_fpt & that)116 robust_fpt& operator+=(const robust_fpt& that) {
117 floating_point_type fpv = this->fpv_ + that.fpv_;
118 if ((!is_neg(this->fpv_) && !is_neg(that.fpv_)) ||
119 (!is_pos(this->fpv_) && !is_pos(that.fpv_))) {
120 this->re_ = (std::max)(this->re_, that.re_) + ROUNDING_ERROR;
121 } else {
122 floating_point_type temp =
123 (this->fpv_ * this->re_ - that.fpv_ * that.re_) / fpv;
124 if (is_neg(temp))
125 temp = -temp;
126 this->re_ = temp + ROUNDING_ERROR;
127 }
128 this->fpv_ = fpv;
129 return *this;
130 }
131
operator -=(const robust_fpt & that)132 robust_fpt& operator-=(const robust_fpt& that) {
133 floating_point_type fpv = this->fpv_ - that.fpv_;
134 if ((!is_neg(this->fpv_) && !is_pos(that.fpv_)) ||
135 (!is_pos(this->fpv_) && !is_neg(that.fpv_))) {
136 this->re_ = (std::max)(this->re_, that.re_) + ROUNDING_ERROR;
137 } else {
138 floating_point_type temp =
139 (this->fpv_ * this->re_ + that.fpv_ * that.re_) / fpv;
140 if (is_neg(temp))
141 temp = -temp;
142 this->re_ = temp + ROUNDING_ERROR;
143 }
144 this->fpv_ = fpv;
145 return *this;
146 }
147
operator *=(const robust_fpt & that)148 robust_fpt& operator*=(const robust_fpt& that) {
149 this->re_ += that.re_ + ROUNDING_ERROR;
150 this->fpv_ *= that.fpv_;
151 return *this;
152 }
153
operator /=(const robust_fpt & that)154 robust_fpt& operator/=(const robust_fpt& that) {
155 this->re_ += that.re_ + ROUNDING_ERROR;
156 this->fpv_ /= that.fpv_;
157 return *this;
158 }
159
operator +(const robust_fpt & that) const160 robust_fpt operator+(const robust_fpt& that) const {
161 floating_point_type fpv = this->fpv_ + that.fpv_;
162 relative_error_type re;
163 if ((!is_neg(this->fpv_) && !is_neg(that.fpv_)) ||
164 (!is_pos(this->fpv_) && !is_pos(that.fpv_))) {
165 re = (std::max)(this->re_, that.re_) + ROUNDING_ERROR;
166 } else {
167 floating_point_type temp =
168 (this->fpv_ * this->re_ - that.fpv_ * that.re_) / fpv;
169 if (is_neg(temp))
170 temp = -temp;
171 re = temp + ROUNDING_ERROR;
172 }
173 return robust_fpt(fpv, re);
174 }
175
operator -(const robust_fpt & that) const176 robust_fpt operator-(const robust_fpt& that) const {
177 floating_point_type fpv = this->fpv_ - that.fpv_;
178 relative_error_type re;
179 if ((!is_neg(this->fpv_) && !is_pos(that.fpv_)) ||
180 (!is_pos(this->fpv_) && !is_neg(that.fpv_))) {
181 re = (std::max)(this->re_, that.re_) + ROUNDING_ERROR;
182 } else {
183 floating_point_type temp =
184 (this->fpv_ * this->re_ + that.fpv_ * that.re_) / fpv;
185 if (is_neg(temp))
186 temp = -temp;
187 re = temp + ROUNDING_ERROR;
188 }
189 return robust_fpt(fpv, re);
190 }
191
operator *(const robust_fpt & that) const192 robust_fpt operator*(const robust_fpt& that) const {
193 floating_point_type fpv = this->fpv_ * that.fpv_;
194 relative_error_type re = this->re_ + that.re_ + ROUNDING_ERROR;
195 return robust_fpt(fpv, re);
196 }
197
operator /(const robust_fpt & that) const198 robust_fpt operator/(const robust_fpt& that) const {
199 floating_point_type fpv = this->fpv_ / that.fpv_;
200 relative_error_type re = this->re_ + that.re_ + ROUNDING_ERROR;
201 return robust_fpt(fpv, re);
202 }
203
sqrt() const204 robust_fpt sqrt() const {
205 return robust_fpt(get_sqrt(fpv_),
206 re_ * static_cast<relative_error_type>(0.5) +
207 ROUNDING_ERROR);
208 }
209
210 private:
211 floating_point_type fpv_;
212 relative_error_type re_;
213 };
214
215 template <typename T>
get_sqrt(const robust_fpt<T> & that)216 robust_fpt<T> get_sqrt(const robust_fpt<T>& that) {
217 return that.sqrt();
218 }
219
220 template <typename T>
is_pos(const robust_fpt<T> & that)221 bool is_pos(const robust_fpt<T>& that) {
222 return that.has_pos_value();
223 }
224
225 template <typename T>
is_neg(const robust_fpt<T> & that)226 bool is_neg(const robust_fpt<T>& that) {
227 return that.has_neg_value();
228 }
229
230 template <typename T>
is_zero(const robust_fpt<T> & that)231 bool is_zero(const robust_fpt<T>& that) {
232 return that.has_zero_value();
233 }
234
235 // robust_dif consists of two not negative values: value1 and value2.
236 // The resulting expression is equal to the value1 - value2.
237 // Subtraction of a positive value is equivalent to the addition to value2
238 // and subtraction of a negative value is equivalent to the addition to
239 // value1. The structure implicitly avoids difference computation.
240 template <typename T>
241 class robust_dif {
242 public:
robust_dif()243 robust_dif() :
244 positive_sum_(0),
245 negative_sum_(0) {}
246
robust_dif(const T & value)247 explicit robust_dif(const T& value) :
248 positive_sum_((value > 0)?value:0),
249 negative_sum_((value < 0)?-value:0) {}
250
robust_dif(const T & pos,const T & neg)251 robust_dif(const T& pos, const T& neg) :
252 positive_sum_(pos),
253 negative_sum_(neg) {}
254
dif() const255 T dif() const {
256 return positive_sum_ - negative_sum_;
257 }
258
pos() const259 T pos() const {
260 return positive_sum_;
261 }
262
neg() const263 T neg() const {
264 return negative_sum_;
265 }
266
operator -() const267 robust_dif<T> operator-() const {
268 return robust_dif(negative_sum_, positive_sum_);
269 }
270
operator +=(const T & val)271 robust_dif<T>& operator+=(const T& val) {
272 if (!is_neg(val))
273 positive_sum_ += val;
274 else
275 negative_sum_ -= val;
276 return *this;
277 }
278
operator +=(const robust_dif<T> & that)279 robust_dif<T>& operator+=(const robust_dif<T>& that) {
280 positive_sum_ += that.positive_sum_;
281 negative_sum_ += that.negative_sum_;
282 return *this;
283 }
284
operator -=(const T & val)285 robust_dif<T>& operator-=(const T& val) {
286 if (!is_neg(val))
287 negative_sum_ += val;
288 else
289 positive_sum_ -= val;
290 return *this;
291 }
292
operator -=(const robust_dif<T> & that)293 robust_dif<T>& operator-=(const robust_dif<T>& that) {
294 positive_sum_ += that.negative_sum_;
295 negative_sum_ += that.positive_sum_;
296 return *this;
297 }
298
operator *=(const T & val)299 robust_dif<T>& operator*=(const T& val) {
300 if (!is_neg(val)) {
301 positive_sum_ *= val;
302 negative_sum_ *= val;
303 } else {
304 positive_sum_ *= -val;
305 negative_sum_ *= -val;
306 swap();
307 }
308 return *this;
309 }
310
operator *=(const robust_dif<T> & that)311 robust_dif<T>& operator*=(const robust_dif<T>& that) {
312 T positive_sum = this->positive_sum_ * that.positive_sum_ +
313 this->negative_sum_ * that.negative_sum_;
314 T negative_sum = this->positive_sum_ * that.negative_sum_ +
315 this->negative_sum_ * that.positive_sum_;
316 positive_sum_ = positive_sum;
317 negative_sum_ = negative_sum;
318 return *this;
319 }
320
operator /=(const T & val)321 robust_dif<T>& operator/=(const T& val) {
322 if (!is_neg(val)) {
323 positive_sum_ /= val;
324 negative_sum_ /= val;
325 } else {
326 positive_sum_ /= -val;
327 negative_sum_ /= -val;
328 swap();
329 }
330 return *this;
331 }
332
333 private:
swap()334 void swap() {
335 (std::swap)(positive_sum_, negative_sum_);
336 }
337
338 T positive_sum_;
339 T negative_sum_;
340 };
341
342 template<typename T>
operator +(const robust_dif<T> & lhs,const robust_dif<T> & rhs)343 robust_dif<T> operator+(const robust_dif<T>& lhs,
344 const robust_dif<T>& rhs) {
345 return robust_dif<T>(lhs.pos() + rhs.pos(), lhs.neg() + rhs.neg());
346 }
347
348 template<typename T>
operator +(const robust_dif<T> & lhs,const T & rhs)349 robust_dif<T> operator+(const robust_dif<T>& lhs, const T& rhs) {
350 if (!is_neg(rhs)) {
351 return robust_dif<T>(lhs.pos() + rhs, lhs.neg());
352 } else {
353 return robust_dif<T>(lhs.pos(), lhs.neg() - rhs);
354 }
355 }
356
357 template<typename T>
operator +(const T & lhs,const robust_dif<T> & rhs)358 robust_dif<T> operator+(const T& lhs, const robust_dif<T>& rhs) {
359 if (!is_neg(lhs)) {
360 return robust_dif<T>(lhs + rhs.pos(), rhs.neg());
361 } else {
362 return robust_dif<T>(rhs.pos(), rhs.neg() - lhs);
363 }
364 }
365
366 template<typename T>
operator -(const robust_dif<T> & lhs,const robust_dif<T> & rhs)367 robust_dif<T> operator-(const robust_dif<T>& lhs,
368 const robust_dif<T>& rhs) {
369 return robust_dif<T>(lhs.pos() + rhs.neg(), lhs.neg() + rhs.pos());
370 }
371
372 template<typename T>
operator -(const robust_dif<T> & lhs,const T & rhs)373 robust_dif<T> operator-(const robust_dif<T>& lhs, const T& rhs) {
374 if (!is_neg(rhs)) {
375 return robust_dif<T>(lhs.pos(), lhs.neg() + rhs);
376 } else {
377 return robust_dif<T>(lhs.pos() - rhs, lhs.neg());
378 }
379 }
380
381 template<typename T>
operator -(const T & lhs,const robust_dif<T> & rhs)382 robust_dif<T> operator-(const T& lhs, const robust_dif<T>& rhs) {
383 if (!is_neg(lhs)) {
384 return robust_dif<T>(lhs + rhs.neg(), rhs.pos());
385 } else {
386 return robust_dif<T>(rhs.neg(), rhs.pos() - lhs);
387 }
388 }
389
390 template<typename T>
operator *(const robust_dif<T> & lhs,const robust_dif<T> & rhs)391 robust_dif<T> operator*(const robust_dif<T>& lhs,
392 const robust_dif<T>& rhs) {
393 T res_pos = lhs.pos() * rhs.pos() + lhs.neg() * rhs.neg();
394 T res_neg = lhs.pos() * rhs.neg() + lhs.neg() * rhs.pos();
395 return robust_dif<T>(res_pos, res_neg);
396 }
397
398 template<typename T>
operator *(const robust_dif<T> & lhs,const T & val)399 robust_dif<T> operator*(const robust_dif<T>& lhs, const T& val) {
400 if (!is_neg(val)) {
401 return robust_dif<T>(lhs.pos() * val, lhs.neg() * val);
402 } else {
403 return robust_dif<T>(-lhs.neg() * val, -lhs.pos() * val);
404 }
405 }
406
407 template<typename T>
operator *(const T & val,const robust_dif<T> & rhs)408 robust_dif<T> operator*(const T& val, const robust_dif<T>& rhs) {
409 if (!is_neg(val)) {
410 return robust_dif<T>(val * rhs.pos(), val * rhs.neg());
411 } else {
412 return robust_dif<T>(-val * rhs.neg(), -val * rhs.pos());
413 }
414 }
415
416 template<typename T>
operator /(const robust_dif<T> & lhs,const T & val)417 robust_dif<T> operator/(const robust_dif<T>& lhs, const T& val) {
418 if (!is_neg(val)) {
419 return robust_dif<T>(lhs.pos() / val, lhs.neg() / val);
420 } else {
421 return robust_dif<T>(-lhs.neg() / val, -lhs.pos() / val);
422 }
423 }
424
425 // Used to compute expressions that operate with sqrts with predefined
426 // relative error. Evaluates expressions of the next type:
427 // sum(i = 1 .. n)(A[i] * sqrt(B[i])), 1 <= n <= 4.
428 template <typename _int, typename _fpt, typename _converter>
429 class robust_sqrt_expr {
430 public:
431 enum MAX_RELATIVE_ERROR {
432 MAX_RELATIVE_ERROR_EVAL1 = 4,
433 MAX_RELATIVE_ERROR_EVAL2 = 7,
434 MAX_RELATIVE_ERROR_EVAL3 = 16,
435 MAX_RELATIVE_ERROR_EVAL4 = 25
436 };
437
438 // Evaluates expression (re = 4 EPS):
439 // A[0] * sqrt(B[0]).
eval1(_int * A,_int * B)440 _fpt eval1(_int* A, _int* B) {
441 _fpt a = convert(A[0]);
442 _fpt b = convert(B[0]);
443 return a * get_sqrt(b);
444 }
445
446 // Evaluates expression (re = 7 EPS):
447 // A[0] * sqrt(B[0]) + A[1] * sqrt(B[1]).
eval2(_int * A,_int * B)448 _fpt eval2(_int* A, _int* B) {
449 _fpt a = eval1(A, B);
450 _fpt b = eval1(A + 1, B + 1);
451 if ((!is_neg(a) && !is_neg(b)) ||
452 (!is_pos(a) && !is_pos(b)))
453 return a + b;
454 return convert(A[0] * A[0] * B[0] - A[1] * A[1] * B[1]) / (a - b);
455 }
456
457 // Evaluates expression (re = 16 EPS):
458 // A[0] * sqrt(B[0]) + A[1] * sqrt(B[1]) + A[2] * sqrt(B[2]).
eval3(_int * A,_int * B)459 _fpt eval3(_int* A, _int* B) {
460 _fpt a = eval2(A, B);
461 _fpt b = eval1(A + 2, B + 2);
462 if ((!is_neg(a) && !is_neg(b)) ||
463 (!is_pos(a) && !is_pos(b)))
464 return a + b;
465 tA[3] = A[0] * A[0] * B[0] + A[1] * A[1] * B[1] - A[2] * A[2] * B[2];
466 tB[3] = 1;
467 tA[4] = A[0] * A[1] * 2;
468 tB[4] = B[0] * B[1];
469 return eval2(tA + 3, tB + 3) / (a - b);
470 }
471
472
473 // Evaluates expression (re = 25 EPS):
474 // A[0] * sqrt(B[0]) + A[1] * sqrt(B[1]) +
475 // A[2] * sqrt(B[2]) + A[3] * sqrt(B[3]).
eval4(_int * A,_int * B)476 _fpt eval4(_int* A, _int* B) {
477 _fpt a = eval2(A, B);
478 _fpt b = eval2(A + 2, B + 2);
479 if ((!is_neg(a) && !is_neg(b)) ||
480 (!is_pos(a) && !is_pos(b)))
481 return a + b;
482 tA[0] = A[0] * A[0] * B[0] + A[1] * A[1] * B[1] -
483 A[2] * A[2] * B[2] - A[3] * A[3] * B[3];
484 tB[0] = 1;
485 tA[1] = A[0] * A[1] * 2;
486 tB[1] = B[0] * B[1];
487 tA[2] = A[2] * A[3] * -2;
488 tB[2] = B[2] * B[3];
489 return eval3(tA, tB) / (a - b);
490 }
491
492 private:
493 _int tA[5];
494 _int tB[5];
495 _converter convert;
496 };
497 } // detail
498 } // polygon
499 } // boost
500
501 #endif // BOOST_POLYGON_DETAIL_VORONOI_ROBUST_FPT
502