view release on metacpan or  search on metacpan

src/boost/polygon/detail/voronoi_robust_fpt.hpp  view on Meta::CPAN

// Boost.Polygon library detail/voronoi_robust_fpt.hpp header file

//          Copyright Andrii Sydorchuk 2010-2012.
// 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)

// See http://www.boost.org for updates, documentation, and revision history.

#ifndef BOOST_POLYGON_DETAIL_VORONOI_ROBUST_FPT
#define BOOST_POLYGON_DETAIL_VORONOI_ROBUST_FPT

#include <cmath>

// Geometry predicates with floating-point variables usually require
// high-precision predicates to retrieve the correct result.
// Epsilon robust predicates give the result within some epsilon relative
// error, but are a lot faster than high-precision predicates.
// To make algorithm robust and efficient epsilon robust predicates are
// used at the first step. In case of the undefined result high-precision
// arithmetic is used to produce required robustness. This approach
// requires exact computation of epsilon intervals within which epsilon
// robust predicates have undefined value.
// There are two ways to measure an error of floating-point calculations:
// relative error and ULPs (units in the last place).
// Let EPS be machine epsilon, then next inequalities have place:
// 1 EPS <= 1 ULP <= 2 EPS (1), 0.5 ULP <= 1 EPS <= 1 ULP (2).
// ULPs are good for measuring rounding errors and comparing values.
// Relative errors are good for computation of general relative
// error of formulas or expressions. So to calculate epsilon
// interval within which epsilon robust predicates have undefined result
// next schema is used:
//     1) Compute rounding errors of initial variables using ULPs;
//     2) Transform ULPs to epsilons using upper bound of the (1);
//     3) Compute relative error of the formula using epsilon arithmetic;
//     4) Transform epsilon to ULPs using upper bound of the (2);
// In case two values are inside undefined ULP range use high-precision
// arithmetic to produce the correct result, else output the result.
// Look at almost_equal function to see how two floating-point variables
// are checked to fit in the ULP range.
// If A has relative error of r(A) and B has relative error of r(B) then:
//     1) r(A + B) <= max(r(A), r(B)), for A * B >= 0;
//     2) r(A - B) <= B*r(A)+A*r(B)/(A-B), for A * B >= 0;
//     2) r(A * B) <= r(A) + r(B);
//     3) r(A / B) <= r(A) + r(B);
// In addition rounding error should be added, that is always equal to
// 0.5 ULP or at most 1 epsilon. As you might see from the above formulas
// subtraction relative error may be extremely large, that's why
// epsilon robust comparator class is used to store floating point values
// and compute subtraction as the final step of the evaluation.
// For further information about relative errors and ULPs try this link:
// http://docs.sun.com/source/806-3568/ncg_goldberg.html

namespace boost {
namespace polygon {
namespace detail {

template <typename T>
T get_sqrt(const T& that) {
  return (std::sqrt)(that);
}

template <typename T>
bool is_pos(const T& that) {
  return that > 0;
}

template <typename T>
bool is_neg(const T& that) {
  return that < 0;
}

template <typename T>
bool is_zero(const T& that) {
  return that == 0;
}

template <typename _fpt>
class robust_fpt {
 public:
  typedef _fpt floating_point_type;
  typedef _fpt relative_error_type;

  // Rounding error is at most 1 EPS.
  enum {
    ROUNDING_ERROR = 1
  };

  robust_fpt() : fpv_(0.0), re_(0.0) {}
  explicit robust_fpt(floating_point_type fpv) :
      fpv_(fpv), re_(0.0) {}
  robust_fpt(floating_point_type fpv, relative_error_type error) :
      fpv_(fpv), re_(error) {}

  floating_point_type fpv() const { return fpv_; }
  relative_error_type re() const { return re_; }
  relative_error_type ulp() const { return re_; }

  robust_fpt& operator=(const robust_fpt& that) {
    this->fpv_ = that.fpv_;
    this->re_ = that.re_;
    return *this;
  }

  bool has_pos_value() const {
    return is_pos(fpv_);
  }