ceph/src/boost/libs/polygon/include/boost/polygon/detail/voronoi_robust_fpt.hpp

   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>
  60 T get_sqrt(const T& that) {
  61   return (std::sqrt)(that);
  62 }
  63
  64 template <typename T>
  65 bool is_pos(const T& that) {
  66   return that > 0;
  67 }
  68
  69 template <typename T>
  70 bool is_neg(const T& that) {
  71   return that < 0;
  72 }
  73
  74 template <typename T>
  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
  90   robust_fpt() : fpv_(0.0), re_(0.0) {}
  91   explicit robust_fpt(floating_point_type fpv) :
  92       fpv_(fpv), re_(0.0) {}
  93   robust_fpt(floating_point_type fpv, relative_error_type error) :
  94       fpv_(fpv), re_(error) {}
  95
  96   floating_point_type fpv() const { return fpv_; }
  97   relative_error_type re() const { return re_; }
  98   relative_error_type ulp() const { return re_; }
  99
 100   robust_fpt& operator=(const robust_fpt& that) {
 101     this->fpv_ = that.fpv_;
 102     this->re_ = that.re_;
 103     return *this;
 104   }
 105
 106   bool has_pos_value() const {
 107     return is_pos(fpv_);
 108   }
 109
 110   bool has_neg_value() const {
 111     return is_neg(fpv_);
 112   }
 113
 114   bool has_zero_value() const {
 115     return is_zero(fpv_);
 116   }
 117
 118   robust_fpt operator-() const {
 119     return robust_fpt(-fpv_, re_);
 120   }
 121
 122   robust_fpt& operator+=(const robust_fpt& that) {
 123     floating_point_type fpv = this->fpv_ + that.fpv_;
 124     if ((!is_neg(this->fpv_) && !is_neg(that.fpv_)) ||
 125         (!is_pos(this->fpv_) && !is_pos(that.fpv_))) {
 126       this->re_ = (std::max)(this->re_, that.re_) + ROUNDING_ERROR;
 127     } else {
 128       floating_point_type temp =
 129         (this->fpv_ * this->re_ - that.fpv_ * that.re_) / fpv;
 130       if (is_neg(temp))
 131         temp = -temp;
 132       this->re_ = temp + ROUNDING_ERROR;
 133     }
 134     this->fpv_ = fpv;
 135     return *this;
 136   }
 137
 138   robust_fpt& operator-=(const robust_fpt& that) {
 139     floating_point_type fpv = this->fpv_ - that.fpv_;
 140     if ((!is_neg(this->fpv_) && !is_pos(that.fpv_)) ||
 141         (!is_pos(this->fpv_) && !is_neg(that.fpv_))) {
 142        this->re_ = (std::max)(this->re_, that.re_) + ROUNDING_ERROR;
 143     } else {
 144       floating_point_type temp =
 145         (this->fpv_ * this->re_ + that.fpv_ * that.re_) / fpv;
 146       if (is_neg(temp))
 147          temp = -temp;
 148       this->re_ = temp + ROUNDING_ERROR;
 149     }
 150     this->fpv_ = fpv;
 151     return *this;
 152   }
 153
 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
 160   robust_fpt& operator/=(const robust_fpt& that) {
 161     this->re_ += that.re_ + ROUNDING_ERROR;
 162     this->fpv_ /= that.fpv_;
 163     return *this;
 164   }
 165
 166   robust_fpt operator+(const robust_fpt& that) const {
 167     floating_point_type fpv = this->fpv_ + that.fpv_;
 168     relative_error_type re;
 169     if ((!is_neg(this->fpv_) && !is_neg(that.fpv_)) ||
 170         (!is_pos(this->fpv_) && !is_pos(that.fpv_))) {
 171       re = (std::max)(this->re_, that.re_) + ROUNDING_ERROR;
 172     } else {
 173       floating_point_type temp =
 174         (this->fpv_ * this->re_ - that.fpv_ * that.re_) / fpv;
 175       if (is_neg(temp))
 176         temp = -temp;
 177       re = temp + ROUNDING_ERROR;
 178     }
 179     return robust_fpt(fpv, re);
 180   }
 181
 182   robust_fpt operator-(const robust_fpt& that) const {
 183     floating_point_type fpv = this->fpv_ - that.fpv_;
 184     relative_error_type re;
 185     if ((!is_neg(this->fpv_) && !is_pos(that.fpv_)) ||
 186         (!is_pos(this->fpv_) && !is_neg(that.fpv_))) {
 187       re = (std::max)(this->re_, that.re_) + ROUNDING_ERROR;
 188     } else {
 189       floating_point_type temp =
 190         (this->fpv_ * this->re_ + that.fpv_ * that.re_) / fpv;
 191       if (is_neg(temp))
 192         temp = -temp;
 193       re = temp + ROUNDING_ERROR;
 194     }
 195     return robust_fpt(fpv, re);
 196   }
 197
 198   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
 204   robust_fpt operator/(const robust_fpt& that) const {
 205     floating_point_type fpv = this->fpv_ / that.fpv_;
 206     relative_error_type re = this->re_ + that.re_ + ROUNDING_ERROR;
 207     return robust_fpt(fpv, re);
 208   }
 209
 210   robust_fpt sqrt() const {
 211     return robust_fpt(get_sqrt(fpv_),
 212                       re_ * static_cast<relative_error_type>(0.5) +
 213                       ROUNDING_ERROR);
 214   }
 215
 216  private:
 217   floating_point_type fpv_;
 218   relative_error_type re_;
 219 };
 220
 221 template <typename T>
 222 robust_fpt<T> get_sqrt(const robust_fpt<T>& that) {
 223   return that.sqrt();
 224 }
 225
 226 template <typename T>
 227 bool is_pos(const robust_fpt<T>& that) {
 228   return that.has_pos_value();
 229 }
 230
 231 template <typename T>
 232 bool is_neg(const robust_fpt<T>& that) {
 233   return that.has_neg_value();
 234 }
 235
 236 template <typename T>
 237 bool is_zero(const robust_fpt<T>& that) {
 238   return that.has_zero_value();
 239 }
 240
 241 // robust_dif consists of two not negative values: value1 and value2.
 242 // The resulting expression is equal to the value1 - value2.
 243 // Subtraction of a positive value is equivalent to the addition to value2
 244 // and subtraction of a negative value is equivalent to the addition to
 245 // value1. The structure implicitly avoids difference computation.
 246 template <typename T>
 247 class robust_dif {
 248  public:
 249   robust_dif() :
 250       positive_sum_(0),
 251       negative_sum_(0) {}
 252
 253   explicit robust_dif(const T& value) :
 254       positive_sum_((value > 0)?value:0),
 255       negative_sum_((value < 0)?-value:0) {}
 256
 257   robust_dif(const T& pos, const T& neg) :
 258       positive_sum_(pos),
 259       negative_sum_(neg) {}
 260
 261   T dif() const {
 262     return positive_sum_ - negative_sum_;
 263   }
 264
 265   T pos() const {
 266     return positive_sum_;
 267   }
 268
 269   T neg() const {
 270     return negative_sum_;
 271   }
 272
 273   robust_dif<T> operator-() const {
 274     return robust_dif(negative_sum_, positive_sum_);
 275   }
 276
 277   robust_dif<T>& operator+=(const T& val) {
 278     if (!is_neg(val))
 279       positive_sum_ += val;
 280     else
 281       negative_sum_ -= val;
 282     return *this;
 283   }
 284
 285   robust_dif<T>& operator+=(const robust_dif<T>& that) {
 286     positive_sum_ += that.positive_sum_;
 287     negative_sum_ += that.negative_sum_;
 288     return *this;
 289   }
 290
 291   robust_dif<T>& operator-=(const T& val) {
 292     if (!is_neg(val))
 293       negative_sum_ += val;
 294     else
 295       positive_sum_ -= val;
 296     return *this;
 297   }
 298
 299   robust_dif<T>& operator-=(const robust_dif<T>& that) {
 300     positive_sum_ += that.negative_sum_;
 301     negative_sum_ += that.positive_sum_;
 302     return *this;
 303   }
 304
 305   robust_dif<T>& operator*=(const T& val) {
 306     if (!is_neg(val)) {
 307       positive_sum_ *= val;
 308       negative_sum_ *= val;
 309     } else {
 310       positive_sum_ *= -val;
 311       negative_sum_ *= -val;
 312       swap();
 313     }
 314     return *this;
 315   }
 316
 317   robust_dif<T>& operator*=(const robust_dif<T>& that) {
 318     T positive_sum = this->positive_sum_ * that.positive_sum_ +
 319                      this->negative_sum_ * that.negative_sum_;
 320     T negative_sum = this->positive_sum_ * that.negative_sum_ +
 321                      this->negative_sum_ * that.positive_sum_;
 322     positive_sum_ = positive_sum;
 323     negative_sum_ = negative_sum;
 324     return *this;
 325   }
 326
 327   robust_dif<T>& operator/=(const T& val) {
 328     if (!is_neg(val)) {
 329       positive_sum_ /= val;
 330       negative_sum_ /= val;
 331     } else {
 332       positive_sum_ /= -val;
 333       negative_sum_ /= -val;
 334       swap();
 335     }
 336     return *this;
 337   }
 338
 339  private:
 340   void swap() {
 341     (std::swap)(positive_sum_, negative_sum_);
 342   }
 343
 344   T positive_sum_;
 345   T negative_sum_;
 346 };
 347
 348 template<typename T>
 349 robust_dif<T> operator+(const robust_dif<T>& lhs,
 350                         const robust_dif<T>& rhs) {
 351   return robust_dif<T>(lhs.pos() + rhs.pos(), lhs.neg() + rhs.neg());
 352 }
 353
 354 template<typename T>
 355 robust_dif<T> operator+(const robust_dif<T>& lhs, const T& rhs) {
 356   if (!is_neg(rhs)) {
 357     return robust_dif<T>(lhs.pos() + rhs, lhs.neg());
 358   } else {
 359     return robust_dif<T>(lhs.pos(), lhs.neg() - rhs);
 360   }
 361 }
 362
 363 template<typename T>
 364 robust_dif<T> operator+(const T& lhs, const robust_dif<T>& rhs) {
 365   if (!is_neg(lhs)) {
 366     return robust_dif<T>(lhs + rhs.pos(), rhs.neg());
 367   } else {
 368     return robust_dif<T>(rhs.pos(), rhs.neg() - lhs);
 369   }
 370 }
 371
 372 template<typename T>
 373 robust_dif<T> operator-(const robust_dif<T>& lhs,
 374                         const robust_dif<T>& rhs) {
 375   return robust_dif<T>(lhs.pos() + rhs.neg(), lhs.neg() + rhs.pos());
 376 }
 377
 378 template<typename T>
 379 robust_dif<T> operator-(const robust_dif<T>& lhs, const T& rhs) {
 380   if (!is_neg(rhs)) {
 381     return robust_dif<T>(lhs.pos(), lhs.neg() + rhs);
 382   } else {
 383     return robust_dif<T>(lhs.pos() - rhs, lhs.neg());
 384   }
 385 }
 386
 387 template<typename T>
 388 robust_dif<T> operator-(const T& lhs, const robust_dif<T>& rhs) {
 389   if (!is_neg(lhs)) {
 390     return robust_dif<T>(lhs + rhs.neg(), rhs.pos());
 391   } else {
 392     return robust_dif<T>(rhs.neg(), rhs.pos() - lhs);
 393   }
 394 }
 395
 396 template<typename T>
 397 robust_dif<T> operator*(const robust_dif<T>& lhs,
 398                         const robust_dif<T>& rhs) {
 399   T res_pos = lhs.pos() * rhs.pos() + lhs.neg() * rhs.neg();
 400   T res_neg = lhs.pos() * rhs.neg() + lhs.neg() * rhs.pos();
 401   return robust_dif<T>(res_pos, res_neg);
 402 }
 403
 404 template<typename T>
 405 robust_dif<T> operator*(const robust_dif<T>& lhs, const T& val) {
 406   if (!is_neg(val)) {
 407     return robust_dif<T>(lhs.pos() * val, lhs.neg() * val);
 408   } else {
 409     return robust_dif<T>(-lhs.neg() * val, -lhs.pos() * val);
 410   }
 411 }
 412
 413 template<typename T>
 414 robust_dif<T> operator*(const T& val, const robust_dif<T>& rhs) {
 415   if (!is_neg(val)) {
 416     return robust_dif<T>(val * rhs.pos(), val * rhs.neg());
 417   } else {
 418     return robust_dif<T>(-val * rhs.neg(), -val * rhs.pos());
 419   }
 420 }
 421
 422 template<typename T>
 423 robust_dif<T> operator/(const robust_dif<T>& lhs, const T& val) {
 424   if (!is_neg(val)) {
 425     return robust_dif<T>(lhs.pos() / val, lhs.neg() / val);
 426   } else {
 427     return robust_dif<T>(-lhs.neg() / val, -lhs.pos() / val);
 428   }
 429 }
 430
 431 // Used to compute expressions that operate with sqrts with predefined
 432 // relative error. Evaluates expressions of the next type:
 433 // sum(i = 1 .. n)(A[i] * sqrt(B[i])), 1 <= n <= 4.
 434 template <typename _int, typename _fpt, typename _converter>
 435 class robust_sqrt_expr {
 436  public:
 437   enum MAX_RELATIVE_ERROR {
 438     MAX_RELATIVE_ERROR_EVAL1 = 4,
 439     MAX_RELATIVE_ERROR_EVAL2 = 7,
 440     MAX_RELATIVE_ERROR_EVAL3 = 16,
 441     MAX_RELATIVE_ERROR_EVAL4 = 25
 442   };
 443
 444   // Evaluates expression (re = 4 EPS):
 445   // A[0] * sqrt(B[0]).
 446   _fpt eval1(_int* A, _int* B) {
 447     _fpt a = convert(A[0]);
 448     _fpt b = convert(B[0]);
 449     return a * get_sqrt(b);
 450   }
 451
 452   // Evaluates expression (re = 7 EPS):
 453   // A[0] * sqrt(B[0]) + A[1] * sqrt(B[1]).
 454   _fpt eval2(_int* A, _int* B) {
 455     _fpt a = eval1(A, B);
 456     _fpt b = eval1(A + 1, B + 1);
 457     if ((!is_neg(a) && !is_neg(b)) ||
 458         (!is_pos(a) && !is_pos(b)))
 459       return a + b;
 460     return convert(A[0] * A[0] * B[0] - A[1] * A[1] * B[1]) / (a - b);
 461   }
 462
 463   // Evaluates expression (re = 16 EPS):
 464   // A[0] * sqrt(B[0]) + A[1] * sqrt(B[1]) + A[2] * sqrt(B[2]).
 465   _fpt eval3(_int* A, _int* B) {
 466     _fpt a = eval2(A, B);
 467     _fpt b = eval1(A + 2, B + 2);
 468     if ((!is_neg(a) && !is_neg(b)) ||
 469         (!is_pos(a) && !is_pos(b)))
 470       return a + b;
 471     tA[3] = A[0] * A[0] * B[0] + A[1] * A[1] * B[1] - A[2] * A[2] * B[2];
 472     tB[3] = 1;
 473     tA[4] = A[0] * A[1] * 2;
 474     tB[4] = B[0] * B[1];
 475     return eval2(tA + 3, tB + 3) / (a - b);
 476   }
 477
 478
 479   // Evaluates expression (re = 25 EPS):
 480   // A[0] * sqrt(B[0]) + A[1] * sqrt(B[1]) +
 481   // A[2] * sqrt(B[2]) + A[3] * sqrt(B[3]).
 482   _fpt eval4(_int* A, _int* B) {
 483     _fpt a = eval2(A, B);
 484     _fpt b = eval2(A + 2, B + 2);
 485     if ((!is_neg(a) && !is_neg(b)) ||
 486         (!is_pos(a) && !is_pos(b)))
 487       return a + b;
 488     tA[0] = A[0] * A[0] * B[0] + A[1] * A[1] * B[1] -
 489             A[2] * A[2] * B[2] - A[3] * A[3] * B[3];
 490     tB[0] = 1;
 491     tA[1] = A[0] * A[1] * 2;
 492     tB[1] = B[0] * B[1];
 493     tA[2] = A[2] * A[3] * -2;
 494     tB[2] = B[2] * B[3];
 495     return eval3(tA, tB) / (a - b);
 496   }
 497
 498  private:
 499   _int tA[5];
 500   _int tB[5];
 501   _converter convert;
 502 };
 503 }  // detail
 504 }  // polygon
 505 }  // boost
 506
 507 #endif  // BOOST_POLYGON_DETAIL_VORONOI_ROBUST_FPT