Subject: [Boost-commit] svn:boost r71681 - in sandbox/SOC/2010/sweepline: . boost/sweepline boost/sweepline/detail libs/sweepline/test
From: sydorchuk.andriy_at_[hidden]
Date: 2011-05-02 17:59:45
Author: asydorchuk
Date: 2011-05-02 17:59:42 EDT (Mon, 02 May 2011)
New Revision: 71681
URL: http://svn.boost.org/trac/boost/changeset/71681
Log:
-added mpt wrapper class;
-added robust fpt class;
-added sqrt expr evaluator class;
-implemented lazy logic for the ppp circle evaluator;
-updated benchmark.
Added:
sandbox/SOC/2010/sweepline/boost/sweepline/detail/mpt_wrapper.hpp (contents, props changed)
sandbox/SOC/2010/sweepline/boost/sweepline/detail/robust_fpt.hpp (contents, props changed)
sandbox/SOC/2010/sweepline/boost/sweepline/detail/sqrt_expr_evaluator.hpp (contents, props changed)
sandbox/SOC/2010/sweepline/libs/sweepline/test/robust_fpt_test.cpp (contents, props changed)
sandbox/SOC/2010/sweepline/libs/sweepline/test/sqrt_expr_evaluator_test.cpp
- copied, changed from r70966, /sandbox/SOC/2010/sweepline/libs/sweepline/test/mpz_arithmetic_test.cpp
Removed:
sandbox/SOC/2010/sweepline/boost/sweepline/detail/mpz_arithmetic.hpp
sandbox/SOC/2010/sweepline/libs/sweepline/test/mpz_arithmetic_test.cpp
Text files modified:
sandbox/SOC/2010/sweepline/Jamfile.v2 | 24
sandbox/SOC/2010/sweepline/boost/sweepline/detail/voronoi_formation.hpp | 1226 ++++++++++++++++-----------------------
sandbox/SOC/2010/sweepline/boost/sweepline/voronoi_sweepline.hpp | 8
sandbox/SOC/2010/sweepline/libs/sweepline/test/Jamfile.v2 | 37 -
sandbox/SOC/2010/sweepline/libs/sweepline/test/benchmark_test.cpp | 40
sandbox/SOC/2010/sweepline/libs/sweepline/test/sqrt_expr_evaluator_test.cpp | 40
sandbox/SOC/2010/sweepline/libs/sweepline/test/sweepline_test.cpp | 106 +-
7 files changed, 643 insertions(+), 838 deletions(-)
Modified: sandbox/SOC/2010/sweepline/Jamfile.v2
==============================================================================
--- sandbox/SOC/2010/sweepline/Jamfile.v2 (original)
+++ sandbox/SOC/2010/sweepline/Jamfile.v2 2011-05-02 17:59:42 EDT (Mon, 02 May 2011)
@@ -3,10 +3,34 @@
# (See accompanying file LICENSE_1_0.txt or copy at
# http://www.boost.org/LICENSE_1_0.txt)
+path-constant GMP_PATH : "c:/Users/Slevin/SWE/Sweepline/gmp" ;
+#path-constant GMP_PATH : "/usr/local" ;
+
+lib gmp
+ : # sources
+ : # requirements
+ <name>gmp <search>$(GMP_PATH)/lib
+ : #default-build
+ : # usage-requirements
+ ;
+
+lib gmpxx
+ : # sources
+ : # requirements
+ <name>gmpxx <search>$(GMP_PATH)/lib
+ : # default-build
+ : # usage-requirements
+ ;
+
+alias gmp_libs : gmp gmpxx : <link>static ;
+
+
project sweepline
: requirements
<include>.
<include>$(BOOST_ROOT)
+ <include>$(GMP_PATH)/include
+ <library>gmp_libs
:
build-dir bin.v2
;
Added: sandbox/SOC/2010/sweepline/boost/sweepline/detail/mpt_wrapper.hpp
==============================================================================
--- (empty file)
+++ sandbox/SOC/2010/sweepline/boost/sweepline/detail/mpt_wrapper.hpp 2011-05-02 17:59:42 EDT (Mon, 02 May 2011)
@@ -0,0 +1,221 @@
+// Boost sweepline/mpt_wrapper.hpp header file
+
+// Copyright Andrii Sydorchuk 2010.
+// 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_SWEEPLINE_MPT_WRAPPER
+#define BOOST_SWEEPLINE_MPT_WRAPPER
+
+#include <cmath>
+#include <string>
+
+namespace boost {
+namespace sweepline {
+namespace detail {
+
+ // Allows to compute expression without allocation of additional memory.
+ // Expression evaluation should contain less then N temporary variables
+ // (e.g. (a1 * b1) + (a2 * b2) - requires two temporary variables to hold
+ // a1 * b1 and a2 * b2). This class is not thread safe, use mpz_class or
+ // mpq_class instead (however there'll be at least 2 times slowdown).
+ template <typename mpt, int N>
+ class mpt_wrapper {
+ public:
+ mpt_wrapper() {}
+
+ explicit mpt_wrapper(int input) : m_(input) {}
+
+ explicit mpt_wrapper(double input) : m_(input) {}
+
+ mpt_wrapper(const mpt& input) : m_(input) {}
+
+ mpt_wrapper(const mpt_wrapper& w) : m_(w.m_) {}
+
+ template <typename mpt2, int N2>
+ mpt_wrapper(const mpt_wrapper<mpt2, N2>& w) : m_(w.m_) {}
+
+ mpt_wrapper& operator=(int that) {
+ m_ = that;
+ return *this;
+ }
+
+ mpt_wrapper& operator=(double that) {
+ m_ = that;
+ return *this;
+ }
+
+ mpt_wrapper& operator=(const mpt_wrapper& that) {
+ m_ = that.m_;
+ return *this;
+ }
+
+ template <typename mpt2, int N2>
+ mpt_wrapper& operator=(const mpt_wrapper<mpt2, N2>& that) {
+ m_ = that.get_mpt();
+ return *this;
+ }
+
+ double get_d() const {
+ return m_.get_d();
+ }
+
+ std::string get_str() const {
+ return m_.get_str();
+ }
+
+ const mpt& get_mpt() const {
+ return m_;
+ }
+
+ bool operator==(const mpt_wrapper& that) const {
+ return m_ == that.m_;
+ }
+
+ bool operator!=(const mpt_wrapper& that) const {
+ return m_ != that.m_;
+ }
+
+ bool operator<(const mpt_wrapper& that) const {
+ return m_ < that.m_;
+ }
+
+ bool operator<=(const mpt_wrapper& that) const {
+ return m_ <= that.m_;
+ }
+
+ bool operator>(const mpt_wrapper& that) const {
+ return m_ > that.m_;
+ }
+
+ bool operator>=(const mpt_wrapper& that) const {
+ return m_ >= that.m_;
+ }
+
+ bool operator==(int that) const {
+ if (that == 0)
+ return m_.__get_mp()->_mp_size == 0;
+ temp_[cur_] = that;
+ return m_ == temp_[cur_].m_;
+ }
+
+ bool operator!=(int that) const {
+ return !(*this == that);
+ }
+
+ bool operator<(int that) const {
+ if (that == 0)
+ return m_.__get_mp()->_mp_size < 0;
+ temp_[cur_] = that;
+ return m_ < temp_[cur_].m_;
+ }
+
+ bool operator<=(int that) const {
+ if (that == 0)
+ return m_.__get_mp()->_mp_size <= 0;
+ temp_[cur_] = that;
+ return m_ <= temp_[cur_].m_;
+ }
+
+ bool operator>(int that) const {
+ if (that == 0)
+ return m_.__get_mp()->_mp_size > 0;
+ temp_[cur_] = that;
+ return m_ > temp_[cur_].m_;
+ }
+
+ bool operator>=(int that) const {
+ if (that == 0)
+ return m_.__get_mp()->_mp_size >= 0;
+ temp_[cur_] = that;
+ return m_ >= temp_[cur_].m_;
+ }
+
+ mpt_wrapper& operator-() const {
+ temp_[cur_].m_ = -this->m_;
+ return temp_[next_cur()];
+ }
+
+ mpt_wrapper& operator+(const mpt_wrapper& that) const {
+ temp_[cur_].m_ = this->m_ + that.m_;
+ return temp_[next_cur()];
+ }
+
+ mpt_wrapper& operator-(const mpt_wrapper& that) const {
+ temp_[cur_].m_ = this->m_ - that.m_;
+ return temp_[next_cur()];
+ }
+
+ mpt_wrapper& operator*(const mpt_wrapper& that) const {
+ temp_[cur_].m_ = this->m_ * that.m_;
+ return temp_[next_cur()];
+ }
+
+ mpt_wrapper& operator/(const mpt_wrapper& that) const {
+ temp_[cur_].m_ = this->m_ / that.m_;
+ return temp_[next_cur()];
+ }
+
+ mpt_wrapper& operator*(double that) const {
+ temp_[cur_].m_ = that;
+ temp_[cur_].m_ *= this->m_;
+ return temp_[next_cur()];
+ }
+
+ mpt_wrapper& operator+=(const mpt_wrapper& that) {
+ this->m_ += that.m_;
+ return *this;
+ }
+
+ mpt_wrapper& operator-=(const mpt_wrapper& that) {
+ this->m_ -= that.m_;
+ return *this;
+ }
+
+ mpt_wrapper& operator*=(const mpt_wrapper& that) {
+ this->m_ *= that.m_;
+ return *this;
+ }
+
+ mpt_wrapper& operator/=(const mpt_wrapper& that) {
+ this->m_ /= that.m_;
+ return *this;
+ }
+
+ mpt_wrapper& get_sqrt() {
+ temp_[cur_].m_ = sqrt(m_);
+ return temp_[next_cur()];
+ }
+
+ private:
+ static int next_cur() {
+ int ret_val = cur_++;
+ if (cur_ == N)
+ cur_ = 0;
+ return ret_val;
+ }
+
+ mpt m_;
+ static int cur_;
+ static mpt_wrapper temp_[N];
+ };
+
+ template <typename mpt, int N>
+ int mpt_wrapper<mpt, N>::cur_ = 0;
+
+ template <typename mpt, int N>
+ mpt_wrapper<mpt, N> mpt_wrapper<mpt, N>::temp_[N];
+
+ template<int N>
+ mpt_wrapper<mpf_class, N>& sqrt(mpt_wrapper<mpf_class, N>& value) {
+ return value.get_sqrt();
+ }
+
+} // detail
+} // sweepline
+} // boost
+
+#endif
Deleted: sandbox/SOC/2010/sweepline/boost/sweepline/detail/mpz_arithmetic.hpp
==============================================================================
--- sandbox/SOC/2010/sweepline/boost/sweepline/detail/mpz_arithmetic.hpp 2011-05-02 17:59:42 EDT (Mon, 02 May 2011)
+++ (empty file)
@@ -1,257 +0,0 @@
-// Boost sweepline/mpz_arithmetic.hpp header file
-
-// Copyright Andrii Sydorchuk 2010.
-// 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_SWEEPLINE_MPZ_ARITHMETIC
-#define BOOST_SWEEPLINE_MPZ_ARITHMETIC
-
-#include <cmath>
-#include <string>
-
-namespace boost {
-namespace sweepline {
-namespace detail {
-
- // Allows to compute expression without allocation of additional memory.
- // Expression evaluation should contain less then N temporary variables
- // (e.g. (a1 * b1) + (a2 * b2) - requires two temporary variables to hold
- // a1 * b1 and a2 * b2). This class is not thread safe, use mpz_class or
- // mpq_class instead (however there'll be at least 2 times slowdown).
- template <typename mpt, int N>
- class mpt_wrapper {
- public:
- mpt_wrapper() {}
-
- explicit mpt_wrapper(int input) : m_(input) {}
-
- explicit mpt_wrapper(double input) : m_(input) {}
-
- mpt_wrapper(const mpt& input) : m_(input) {}
-
- mpt_wrapper(const mpt_wrapper& w) : m_(w.m_) {
- cur_ = 0;
- }
-
- mpt_wrapper& operator=(int that) {
- m_ = that;
- return *this;
- }
-
- mpt_wrapper& operator=(double that) {
- m_ = that;
- return *this;
- }
-
- mpt_wrapper& operator=(const mpt_wrapper& that) {
- m_ = that.m_;
- return *this;
- }
-
- double get_d() const {
- return m_.get_d();
- }
-
- std::string get_str() const {
- return m_.get_str();
- }
-
- bool operator==(const mpt_wrapper& that) const {
- return m_ == that.m_;
- }
-
- bool operator!=(const mpt_wrapper& that) const {
- return m_ != that.m_;
- }
-
- bool operator<(const mpt_wrapper& that) const {
- return m_ < that.m_;
- }
-
- bool operator<=(const mpt_wrapper& that) const {
- return m_ <= that.m_;
- }
-
- bool operator>(const mpt_wrapper& that) const {
- return m_ > that.m_;
- }
-
- bool operator>=(const mpt_wrapper& that) const {
- return m_ >= that.m_;
- }
-
- bool operator==(int that) const {
- temp_[cur_] = that;
- return m_ == temp_[cur_].m_;
- }
-
- bool operator!=(int that) const {
- temp_[cur_] = that;
- return m_ != temp_[cur_].m_;
- }
-
- bool operator<(int that) const {
- temp_[cur_] = that;
- return m_ < temp_[cur_].m_;
- }
-
- bool operator<=(int that) const {
- temp_[cur_] = that;
- return m_ <= temp_[cur_].m_;
- }
-
- bool operator>(int that) const {
- temp_[cur_] = that;
- return m_ > temp_[cur_].m_;
- }
-
- bool operator>=(int that) const {
- temp_[cur_] = that;
- return m_ >= temp_[cur_].m_;
- }
-
- mpt_wrapper& operator-() const {
- temp_[cur_].m_ = -this->m_;
- return temp_[next_cur()];
- }
-
- mpt_wrapper& operator+(const mpt_wrapper& that) const {
- temp_[cur_].m_ = this->m_ + that.m_;
- return temp_[next_cur()];
- }
-
- mpt_wrapper& operator-(const mpt_wrapper& that) const {
- temp_[cur_].m_ = this->m_ - that.m_;
- return temp_[next_cur()];
- }
-
- mpt_wrapper& operator*(const mpt_wrapper& that) const {
- temp_[cur_].m_ = this->m_ * that.m_;
- return temp_[next_cur()];
- }
-
- mpt_wrapper& operator*(double that) const {
- temp_[cur_].m_ = that;
- temp_[cur_].m_ *= this->m_;
- return temp_[next_cur()];
- }
-
- mpt_wrapper& operator+=(const mpt_wrapper& that) {
- this->m_ += that.m_;
- return *this;
- }
-
- mpt_wrapper& operator-=(const mpt_wrapper& that) {
- this->m_ -= that.m_;
- return *this;
- }
-
- mpt_wrapper& operator*=(const mpt_wrapper& that) {
- this->m_ *= that.m_;
- return *this;
- }
-
- private:
- static int next_cur() {
- int ret_val = cur_++;
- if (cur_ == N)
- cur_ = 0;
- return ret_val;
- }
-
- mpt m_;
- static int cur_;
- static mpt_wrapper temp_[N];
- };
-
- template <typename mpt, int N>
- int mpt_wrapper<mpt, N>::cur_ = 0;
-
- template <typename mpt, int N>
- mpt_wrapper<mpt, N> mpt_wrapper<mpt, N>::temp_[N];
-
- template <int N>
- struct sqr_expr_evaluator {
- template <typename mpt>
- static double eval(mpt *A, mpt *B);
- };
-
- // Evaluates expression:
- // A[0] * sqrt(B[0]) + A[1] * sqrt(B[1]) with
- // 7 * EPS relative error in the worst case.
- template <>
- struct sqr_expr_evaluator<2> {
- template <typename mpt>
- static double eval(mpt *A, mpt *B) {
-#ifndef THREAD_SAFETY
- static
-#endif
- mpt numerator;
- double lhs = A[0].get_d() * sqrt(B[0].get_d());
- double rhs = A[1].get_d() * sqrt(B[1].get_d());
- if ((lhs >= 0 && rhs >= 0) || (lhs <= 0 && rhs <= 0))
- return lhs + rhs;
- numerator = A[0] * A[0] * B[0] - A[1] * A[1] * B[1];
- return numerator.get_d() / (lhs - rhs);
- }
- };
-
- // Evaluates expression:
- // A[0] * sqrt(B[0]) + A[1] * sqrt(B[1]) + A[2] * sqrt(B[2])
- // with 16 * EPS relative error in the worst case.
- template <>
- struct sqr_expr_evaluator<3> {
- template <typename mpt>
- static double eval(mpt *A, mpt *B) {
-#ifndef THREAD_SAFETY
- static
-#endif
- mpt cA[2], cB[2];
- double lhs = sqr_expr_evaluator<2>::eval<mpt>(A, B);
- double rhs = A[2].get_d() * sqrt(B[2].get_d());
- if ((lhs >= 0 && rhs >= 0) || (lhs <= 0 && rhs <= 0))
- return lhs + rhs;
- cA[0] = A[0] * A[0] * B[0] + A[1] * A[1] * B[1];
- cA[0] -= A[2] * A[2] * B[2];
- cB[0] = 1;
- cA[1] = A[0] * A[1] * 2;
- cB[1] = B[0] * B[1];
- return sqr_expr_evaluator<2>::eval<mpt>(cA, cB) / (lhs - rhs);
- }
- };
-
- // Evaluates expression:
- // A[0] * sqrt(B[0]) + A[1] * sqrt(B[1]) + A[2] * sqrt(B[2]) + A[3] * sqrt(B[3])
- // with 25 * EPS relative error in the worst case.
- template <>
- struct sqr_expr_evaluator<4> {
- template <typename mpt>
- static double eval(mpt *A, mpt *B) {
-#ifndef THREAD_SAFETY
- static
-#endif
- mpt cA[3], cB[3];
- double lhs = sqr_expr_evaluator<2>::eval<mpt>(A, B);
- double rhs = sqr_expr_evaluator<2>::eval<mpt>(A + 2, B + 2);
- if ((lhs >= 0 && rhs >= 0) || (lhs <= 0 && rhs <= 0))
- return lhs + rhs;
- cA[0] = A[0] * A[0] * B[0] + A[1] * A[1] * B[1];
- cA[0] -= A[2] * A[2] * B[2] + A[3] * A[3] * B[3];
- cB[0] = 1;
- cA[1] = A[0] * A[1] * 2;
- cB[1] = B[0] * B[1];
- cA[2] = A[2] * A[3] * -2;
- cB[2] = B[2] * B[3];
- return sqr_expr_evaluator<3>::eval<mpt>(cA, cB) / (lhs - rhs);
- }
- };
-
-} // detail
-} // sweepline
-} // boost
-
-#endif
Added: sandbox/SOC/2010/sweepline/boost/sweepline/detail/robust_fpt.hpp
==============================================================================
--- (empty file)
+++ sandbox/SOC/2010/sweepline/boost/sweepline/detail/robust_fpt.hpp 2011-05-02 17:59:42 EDT (Mon, 02 May 2011)
@@ -0,0 +1,458 @@
+// Boost sweepline/mpt_wrapper.hpp header file
+
+// Copyright Andrii Sydorchuk 2010.
+// 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 VORONOI_SWEEPLINE_ROBUST_FPT
+#define VORONOI_SWEEPLINE_ROBUST_FPT
+
+namespace boost {
+namespace sweepline {
+namespace detail {
+
+ // Represents the result of the epsilon robust predicate.
+ // If the result is undefined some further processing is usually required.
+ enum kPredicateResult {
+ LESS = -1,
+ UNDEFINED = 0,
+ MORE = 1,
+ };
+
+ // If two floating-point numbers in the same format are ordered (x < y),
+ // then they are ordered the same way when their bits are reinterpreted as
+ // sign-magnitude integers. Values are considered to be almost equal if
+ // their integer reinterpretatoins differ in not more than maxUlps units.
+ template <typename T>
+ static bool almost_equal(T a, T b, unsigned int ulps);
+
+ template<>
+ bool almost_equal<float>(float a, float b, unsigned int maxUlps) {
+ unsigned int ll_a, ll_b;
+
+ // Reinterpret double bits as 32-bit signed integer.
+ memcpy(&ll_a, &a, sizeof(float));
+ memcpy(&ll_b, &b, sizeof(float));
+
+ if (ll_a < 0x80000000)
+ ll_a = 0x80000000 - ll_a;
+ if (ll_b < 0x80000000)
+ ll_b = 0x80000000 - ll_b;
+
+ if (ll_a > ll_b)
+ return ll_a - ll_b <= maxUlps;
+ return ll_b - ll_a <= maxUlps;
+ }
+
+ template<>
+ bool almost_equal<double>(double a, double b, unsigned int maxUlps) {
+ unsigned long long ll_a, ll_b;
+
+ // Reinterpret double bits as 64-bit signed integer.
+ memcpy(&ll_a, &a, sizeof(double));
+ memcpy(&ll_b, &b, sizeof(double));
+
+ // Positive 0.0 is integer zero. Negative 0.0 is 0x8000000000000000.
+ // Map negative zero to an integer zero representation - making it
+ // identical to positive zero - the smallest negative number is
+ // represented by negative one, and downwards from there.
+ if (ll_a < 0x8000000000000000ULL)
+ ll_a = 0x8000000000000000ULL - ll_a;
+ if (ll_b < 0x8000000000000000ULL)
+ ll_b = 0x8000000000000000ULL - ll_b;
+
+ // Compare 64-bit signed integer representations of input values.
+ // Difference in 1 Ulp is equivalent to a relative error of between
+ // 1/4,000,000,000,000,000 and 1/8,000,000,000,000,000.
+ if (ll_a > ll_b)
+ return ll_a - ll_b <= maxUlps;
+ return ll_b - ll_a <= maxUlps;
+ }
+
+ template <typename T>
+ double get_d(const T& value) {
+ return value.get_d();
+ }
+
+ template <>
+ double get_d(const float& value) {
+ return value;
+ }
+
+ template <>
+ double get_d(const double& value) {
+ return value;
+ }
+
+ template <typename _fpt>
+ class robust_fpt {
+ public:
+ typedef _fpt floating_point_type;
+ typedef double relative_error_type;
+
+ // Rounding error is at most 1 EPS.
+ static const relative_error_type ROUNDING_ERROR;
+
+ // Constructors.
+ robust_fpt() : fpv_(0.0), re_(0) {}
+ explicit robust_fpt(int fpv) : fpv_(fpv), re_(0) {}
+ explicit robust_fpt(floating_point_type fpv, bool rounded = true) : fpv_(fpv) {
+ re_ = rounded ? ROUNDING_ERROR : 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_; }
+
+ bool operator==(double that) const {
+ if (that == 0)
+ return this->fpv_ == that;
+ return almost_equal(this->fpv_, that,
+ static_cast<unsigned int>(this->ulp()));
+ }
+
+ bool operator!=(double that) const {
+ return !(*this == that);
+ }
+
+ bool operator<(double that) const {
+ if (*this == that)
+ return false;
+ return this->fpv_ < that;
+ }
+
+ bool operator<=(double that) const {
+ if (*this == that)
+ return true;
+ return this->fpv_ < that;
+ }
+
+ bool operator>(double that) const {
+ if (*this == that)
+ return false;
+ return this->fpv_ > that;
+ }
+
+ bool operator>=(double that) const {
+ if (*this == that)
+ return true;
+ return this->fpv_ > that;
+ }
+
+ bool operator==(const robust_fpt &that) const {
+ unsigned int ulp = static_cast<unsigned int>(ceil(this->re_ + that.re_));
+ return almost_equal(this->fpv_, that.fpv_, ulp);
+ }
+
+ bool operator!=(const robust_fpt &that) const {
+ return !(*this == that);
+ }
+
+ bool operator<(const robust_fpt &that) const {
+ if (*this == that)
+ return false;
+ return this->fpv_ < that.fpv_;
+ }
+
+ bool operator>(const robust_fpt &that) const {
+ return that < *this;
+ }
+
+ bool operator<=(const robust_fpt &that) const {
+ return !(that < *this);
+ }
+
+ bool operator>=(const robust_fpt &that) const {
+ return !(*this < that);
+ }
+
+ bool operator-() const {
+ return robust_fpt(-fpv_, re_);
+ }
+
+ robust_fpt& operator=(const robust_fpt &that) {
+ this->fpv_ = that.fpv_;
+ this->relative_error_ = that.re_;
+ return *this;
+ }
+
+ robust_fpt& operator+=(const robust_fpt &that) {
+ floating_point_type fpv = this->fpv_ + that.fpv_;
+ if ((this->fpv_ >= 0 && that.fpv_ >= 0) ||
+ (this->fpv_ <= 0 && that.fpv_ <= 0))
+ this->re_ = (std::max)(this->re_, that.re_) + ROUNDING_ERROR;
+ else {
+ floating_point_type temp = (this->fpv_ * this->re_ - that.fpv_ * that.re_) / fpv;
+ this->re_ = fabs(get_d(temp)) + ROUNDING_ERROR;
+ }
+ this->fpv_ = fpv;
+ return *this;
+ }
+
+ robust_fpt& operator-=(const robust_fpt &that) {
+ floating_point_type fpv = this->fpv_ - that.fpv_;
+ if ((this->fpv_ >= 0 && that.fpv_ <= 0) ||
+ (this->fpv_ <= 0 && that.fpv_ >= 0))
+ this->re_ = (std::max)(this->re_, that.re_) + ROUNDING_ERROR;
+ else {
+ floating_point_type temp = (this->fpv_ * this->re_ + that.fpv_ * that.re_) / fpv;
+ this->re_ = fabs(get_d(temp)) + ROUNDING_ERROR;
+ }
+ this->fpv_ = fpv;
+ return *this;
+ }
+
+ robust_fpt& operator*=(const robust_fpt &that) {
+ this->re_ += that.re_ + ROUNDING_ERROR;
+ this->fpv_ *= that.fpv_;
+ return *this;
+ }
+
+ robust_fpt& operator/=(const robust_fpt &that) {
+ this->re_ += that.re_ + ROUNDING_ERROR;
+ this->fpv_ /= that.fpv_;
+ return *this;
+ }
+
+ robust_fpt operator+(const robust_fpt &that) const {
+ floating_point_type fpv = this->fpv_ + that.fpv_;
+ relative_error_type re;
+ if ((this->fpv_ >= 0 && that.fpv_ >= 0) ||
+ (this->fpv_ <= 0 && that.fpv_ <= 0))
+ re = (std::max)(this->re_, that.re_) + ROUNDING_ERROR;
+ else {
+ floating_point_type temp = (this->fpv_ * this->re_ - that.fpv_ * that.re_) / fpv;
+ re = fabs(get_d(temp)) + ROUNDING_ERROR;
+ }
+ return robust_fpt(fpv, re);
+ }
+
+ robust_fpt operator-(const robust_fpt &that) const {
+ floating_point_type fpv = this->fpv_ - that.fpv_;
+ relative_error_type re;
+ if ((this->fpv_ >= 0 && that.fpv_ <= 0) ||
+ (this->fpv_ <= 0 && that.fpv_ >= 0))
+ re = (std::max)(this->re_, that.re_) + ROUNDING_ERROR;
+ else {
+ floating_point_type temp = (this->fpv_ * this->re_ + that.fpv_ * that.re_) / fpv;
+ re = fabs(get_d(temp)) + ROUNDING_ERROR;
+ }
+ return robust_fpt(fpv, re);
+ }
+
+ robust_fpt operator*(const robust_fpt &that) const {
+ floating_point_type fpv = this->fpv_ * that.fpv_;
+ relative_error_type re = this->re_ + that.re_ + ROUNDING_ERROR;
+ return robust_fpt(fpv, re);
+ }
+
+ robust_fpt operator/(const robust_fpt &that) const {
+ floating_point_type fpv = this->fpv_ / that.fpv_;
+ relative_error_type re = this->re_ + that.re_ + ROUNDING_ERROR;
+ return robust_fpt(fpv, re);
+ }
+
+ robust_fpt get_sqrt() const {
+ return robust_fpt(std::sqrt(fpv_), re_ * 0.5 + ROUNDING_ERROR);
+ }
+
+ private:
+ floating_point_type fpv_;
+ relative_error_type re_;
+ };
+
+ template <typename T>
+ const typename robust_fpt<T>::relative_error_type robust_fpt<T>::ROUNDING_ERROR = 1;
+
+ // Class used to make computations with epsilon relative error.
+ // ERC consists of two values: value1 and value2, both positive.
+ // The resulting expression is equal to the value1 - value2.
+ // The main idea is to represent any expression that consists of
+ // addition, substraction, multiplication and division operations
+ // to avoid using substraction. Substraction of a positive value
+ // is equivalent to the addition to value2 and substraction of
+ // a negative value is equivalent to the addition to value1.
+ // Cons: ERC gives error relative not to the resulting value,
+ // but relative to some expression instead. Example:
+ // center_x = 100, ERC's value1 = 10^20, value2 = 10^20,
+ // center_x = 1000, ERC's value3 = 10^21, value4 = 10^21,
+ // such two centers are considered equal(
+ // value1 + value4 = value2 + value3), while they are not.
+ // Pros: ERC is much faster then approaches based on the use
+ // of high-precision libraries. However this will give correct
+ // answer for the previous example.
+ // Solution: Use ERCs in case of defined comparison results and use
+ // high-precision libraries for undefined results.
+ template <typename T>
+ class epsilon_robust_comparator {
+ public:
+ epsilon_robust_comparator() :
+ positive_sum_(0),
+ negative_sum_(0) {}
+
+ epsilon_robust_comparator(const T &value) :
+ positive_sum_((value>0)?value:0),
+ negative_sum_((value<0)?-value:0) {}
+
+ epsilon_robust_comparator(const T &pos, const T &neg) :
+ positive_sum_(pos),
+ negative_sum_(neg) {}
+
+ T dif() const {
+ return positive_sum_ - negative_sum_;
+ }
+
+ T pos() const {
+ return positive_sum_;
+ }
+
+ T neg() const {
+ return negative_sum_;
+ }
+
+ // Equivalent to the unary minus.
+ void swap() {
+ (std::swap)(positive_sum_, negative_sum_);
+ }
+
+ bool abs() {
+ if (positive_sum_ < negative_sum_) {
+ swap();
+ return true;
+ }
+ return false;
+ }
+
+ epsilon_robust_comparator<T> &operator+=(const T &val) {
+ if (val >= 0)
+ positive_sum_ += val;
+ else
+ negative_sum_ -= val;
+ return *this;
+ }
+
+ epsilon_robust_comparator<T> &operator+=(
+ const epsilon_robust_comparator<T> &erc) {
+ positive_sum_ += erc.positive_sum_;
+ negative_sum_ += erc.negative_sum_;
+ return *this;
+ }
+
+ epsilon_robust_comparator<T> &operator-=(const T &val) {
+ if (val >= 0)
+ negative_sum_ += val;
+ else
+ positive_sum_ -= val;
+ return *this;
+ }
+
+ epsilon_robust_comparator<T> &operator-=(
+ const epsilon_robust_comparator<T> &erc) {
+ positive_sum_ += erc.negative_sum_;
+ negative_sum_ += erc.positive_sum_;
+ return *this;
+ }
+
+ epsilon_robust_comparator<T> &operator*=(const T &val) {
+ positive_sum_ *= fabs(val);
+ negative_sum_ *= fabs(val);
+ if (val < 0) {
+ swap();
+ }
+ return *this;
+ }
+
+ epsilon_robust_comparator<T> &operator/=(const T &val) {
+ positive_sum_ /= fabs(val);
+ negative_sum_ /= fabs(val);
+ if (val < 0) {
+ swap();
+ }
+ return *this;
+ }
+
+ // Compare predicate with value using ulp precision.
+ kPredicateResult compare(T value, int ulp = 0) const {
+ T lhs = positive_sum_ - ((value < 0) ? value : 0);
+ T rhs = negative_sum_ + ((value > 0) ? value : 0);
+ if (almost_equal(lhs, rhs, ulp))
+ return UNDEFINED;
+ return (lhs < rhs) ? LESS : MORE;
+ }
+
+ // Compare two predicats using ulp precision.
+ kPredicateResult compare(const epsilon_robust_comparator &rc,
+ int ulp = 0) const {
+ T lhs = positive_sum_ + rc.neg();
+ T rhs = negative_sum_ + rc.pos();
+ if (almost_equal(lhs, rhs, ulp))
+ return UNDEFINED;
+ return (lhs < rhs) ? LESS : MORE;
+ }
+
+ // Check whether comparison has undefined result.
+ bool compares_undefined(const epsilon_robust_comparator &rc,
+ int ulp) const {
+ return this->compare(rc, ulp) == UNDEFINED;
+ }
+
+ private:
+ T positive_sum_;
+ T negative_sum_;
+ };
+
+ template<typename T>
+ inline epsilon_robust_comparator<T> operator+(
+ const epsilon_robust_comparator<T>& lhs,
+ const epsilon_robust_comparator<T>& rhs) {
+ return epsilon_robust_comparator<T>(lhs.pos() + rhs.pos(),
+ lhs.neg() + rhs.neg());
+ }
+
+ template<typename T>
+ inline epsilon_robust_comparator<T> operator-(
+ const epsilon_robust_comparator<T>& lhs,
+ const epsilon_robust_comparator<T>& rhs) {
+ return epsilon_robust_comparator<T>(lhs.pos() - rhs.neg(),
+ lhs.neg() + rhs.pos());
+ }
+
+ template<typename T>
+ inline epsilon_robust_comparator<T> operator*(
+ const epsilon_robust_comparator<T>& lhs,
+ const epsilon_robust_comparator<T>& rhs) {
+ double res_pos = lhs.pos() * rhs.pos() + lhs.neg() * rhs.neg();
+ double res_neg = lhs.pos() * rhs.neg() + lhs.neg() * rhs.pos();
+ return epsilon_robust_comparator<T>(res_pos, res_neg);
+ }
+
+ template<typename T>
+ inline epsilon_robust_comparator<T> operator*(
+ const epsilon_robust_comparator<T>& lhs, const T& val) {
+ if (val >= 0)
+ return epsilon_robust_comparator<T>(lhs.pos() * val,
+ lhs.neg() * val);
+ return epsilon_robust_comparator<T>(-lhs.neg() * val,
+ -lhs.pos() * val);
+ }
+
+ template<typename T>
+ inline epsilon_robust_comparator<T> operator*(
+ const T& val, const epsilon_robust_comparator<T>& rhs) {
+ if (val >= 0)
+ return epsilon_robust_comparator<T>(val * rhs.pos(),
+ val * rhs.neg());
+ return epsilon_robust_comparator<T>(-val * rhs.neg(),
+ -val * rhs.pos());
+ }
+
+} // detail
+} // sweepline
+} // boost
+
+#endif
\ No newline at end of file
Added: sandbox/SOC/2010/sweepline/boost/sweepline/detail/sqrt_expr_evaluator.hpp
==============================================================================
--- (empty file)
+++ sandbox/SOC/2010/sweepline/boost/sweepline/detail/sqrt_expr_evaluator.hpp 2011-05-02 17:59:42 EDT (Mon, 02 May 2011)
@@ -0,0 +1,124 @@
+// Boost sweepline/sqr_expr_evaluator.hpp header file
+
+// Copyright Andrii Sydorchuk 2010.
+// 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_SWEEPLINE_SQRT_EXPR_EVALUATOR
+#define BOOST_SWEEPLINE_SQRT_EXPR_EVALUATOR
+
+namespace boost {
+namespace sweepline {
+namespace detail {
+
+ template <int N>
+ struct sqr_expr_evaluator {
+ template <typename mpt, typename mpf>
+ static mpf eval(mpt *A, mpt *B);
+ };
+
+ // Evaluates expression:
+ // A[0] * sqrt(B[0]).
+ template <>
+ struct sqr_expr_evaluator<1> {
+ template <typename mpt, typename mpf>
+ static mpf eval(mpt *A, mpt *B) {
+#ifndef THREAD_SAFETY
+ static
+#endif
+ mpf a, b, ret_val;
+ a = A[0];
+ b = B[0];
+ ret_val = a * sqrt(b);
+ return ret_val;
+ }
+ };
+
+ // Evaluates expression:
+ // A[0] * sqrt(B[0]) + A[1] * sqrt(B[1]) with
+ // 7 * EPS relative error in the worst case.
+ template <>
+ struct sqr_expr_evaluator<2> {
+ template <typename mpt, typename mpf>
+ static mpf eval(mpt *A, mpt *B) {
+#ifndef THREAD_SAFETY
+ static
+#endif
+ mpf lhs, rhs, numerator;
+ lhs = sqr_expr_evaluator<1>::eval<mpt, mpf>(A, B);
+ rhs = sqr_expr_evaluator<1>::eval<mpt, mpf>(A + 1, B + 1);
+ if ((lhs >= 0 && rhs >= 0) || (lhs <= 0 && rhs <= 0))
+ return lhs + rhs;
+ numerator = A[0] * A[0] * B[0] - A[1] * A[1] * B[1];
+ return numerator / (lhs - rhs);
+ }
+ };
+
+ // Evaluates expression:
+ // A[0] * sqrt(B[0]) + A[1] * sqrt(B[1]) + A[2] * sqrt(B[2])
+ // with 16 * EPS relative error in the worst case.
+ template <>
+ struct sqr_expr_evaluator<3> {
+ template <typename mpt, typename mpf>
+ static mpf eval(mpt *A, mpt *B) {
+#ifndef THREAD_SAFETY
+ static
+#endif
+ mpt cA[2], cB[2];
+#ifndef THREAD_SAFETY
+ static
+#endif
+ mpf lhs, rhs, numer;
+ lhs = sqr_expr_evaluator<2>::eval<mpt, mpf>(A, B);
+ rhs = sqr_expr_evaluator<1>::eval<mpt, mpf>(A + 2, B + 2);
+ if ((lhs >= 0 && rhs >= 0) || (lhs <= 0 && rhs <= 0))
+ return lhs + rhs;
+ cA[0] = A[0] * A[0] * B[0] + A[1] * A[1] * B[1];
+ cA[0] -= A[2] * A[2] * B[2];
+ cB[0] = 1;
+ cA[1] = A[0] * A[1] * 2;
+ cB[1] = B[0] * B[1];
+ numer = sqr_expr_evaluator<2>::eval<mpt, mpf>(cA, cB);
+ return numer / (lhs - rhs);
+ }
+ };
+
+ // Evaluates expression:
+ // A[0] * sqrt(B[0]) + A[1] * sqrt(B[1]) + A[2] * sqrt(B[2]) + A[3] * sqrt(B[3])
+ // with 25 * EPS relative error in the worst case.
+ template <>
+ struct sqr_expr_evaluator<4> {
+ template <typename mpt, typename mpf>
+ static mpf eval(mpt *A, mpt *B) {
+#ifndef THREAD_SAFETY
+ static
+#endif
+ mpt cA[3], cB[3];
+#ifndef THREAD_SAFETY
+ static
+#endif
+ mpf lhs, rhs, numer;
+ lhs = sqr_expr_evaluator<2>::eval<mpt, mpf>(A, B);
+ rhs = sqr_expr_evaluator<2>::eval<mpt, mpf>(A + 2, B + 2);
+ if ((lhs >= 0 && rhs >= 0) || (lhs <= 0 && rhs <= 0))
+ return lhs + rhs;
+ cA[0] = A[0] * A[0] * B[0] + A[1] * A[1] * B[1];
+ cA[0] -= A[2] * A[2] * B[2] + A[3] * A[3] * B[3];
+ cB[0] = 1;
+ cA[1] = A[0] * A[1] * 2;
+ cB[1] = B[0] * B[1];
+ cA[2] = A[2] * A[3] * -2;
+ cB[2] = B[2] * B[3];
+ numer = sqr_expr_evaluator<3>::eval<mpt, mpf>(cA, cB);
+ return numer / (lhs - rhs);
+ }
+ };
+
+} // detail
+} // sweepline
+} // boost
+
+#endif
\ No newline at end of file
Modified: sandbox/SOC/2010/sweepline/boost/sweepline/detail/voronoi_formation.hpp
==============================================================================
--- sandbox/SOC/2010/sweepline/boost/sweepline/detail/voronoi_formation.hpp (original)
+++ sandbox/SOC/2010/sweepline/boost/sweepline/detail/voronoi_formation.hpp 2011-05-02 17:59:42 EDT (Mon, 02 May 2011)
@@ -38,14 +38,6 @@
LEFT_ORIENTATION = 1,
};
- // Represents the result of the epsilon robust predicate.
- // If the result is undefined some further processing is usually required.
- enum kPredicateResult {
- LESS = -1,
- UNDEFINED = 0,
- MORE = 1,
- };
-
// Site event type.
// Occurs when the sweepline sweeps over one of the initial sites:
// 1) point site;
@@ -348,7 +340,7 @@
typedef typename std::map< Key, Value, NodeComparer >::iterator
beach_line_iterator;
- circle_event() : is_active_(true) {}
+ circle_event() : denom_(1), is_active_(true) {}
circle_event(coordinate_type c_x,
coordinate_type c_y,
@@ -452,16 +444,28 @@
return center_x_.dif() / denom_.dif();
}
+ void x(coordinate_type center_x) {
+ center_x_ = center_x;
+ }
+
// Evaluate y-coordinate of the center of the circle.
coordinate_type y() const {
return center_y_.dif() / denom_.dif();
}
+ void y(coordinate_type center_y) {
+ center_y_ = center_y;
+ }
+
// Evaluate x-coordinate of the rightmost point of the circle.
coordinate_type lower_x() const {
return lower_x_.dif() / denom_.dif();
}
+ void lower_x(coordinate_type lower_x) {
+ lower_x_ = lower_x;
+ }
+
point_2d_type center() const {
return make_point_2d(x(), y());
}
@@ -642,96 +646,6 @@
}
}
- // If two floating-point numbers in the same format are ordered (x < y),
- // then they are ordered the same way when their bits are reinterpreted as
- // sign-magnitude integers. Values are considered to be almost equal if
- // their integer reinterpretatoins differ in not more than maxUlps units.
- static inline bool almost_equal(double a, double b,
- unsigned int maxUlps) {
- polygon_ulong_long_type ll_a, ll_b;
-
- // Reinterpret double bits as 64-bit signed integer.
- memcpy(&ll_a, &a, sizeof(double));
- memcpy(&ll_b, &b, sizeof(double));
-
- // Positive 0.0 is integer zero. Negative 0.0 is 0x8000000000000000.
- // Map negative zero to an integer zero representation - making it
- // identical to positive zero - the smallest negative number is
- // represented by negative one, and downwards from there.
- if (ll_a < 0x8000000000000000ULL)
- ll_a = 0x8000000000000000ULL - ll_a;
- if (ll_b < 0x8000000000000000ULL)
- ll_b = 0x8000000000000000ULL - ll_b;
-
- // Compare 64-bit signed integer representations of input values.
- // Difference in 1 Ulp is equivalent to a relative error of between
- // 1/4,000,000,000,000,000 and 1/8,000,000,000,000,000.
- if (ll_a > ll_b)
- return ll_a - ll_b <= maxUlps;
- return ll_b - ll_a <= maxUlps;
- }
-
- // Robust orientation test. Works correctly for any input type that
- // can be casted without lose of data to the integer type within a range
- // [-2^32, 2^32-1].
- // Arguments: dif_x1_, dif_y1 - coordinates of the first vector.
- // dif_x2_, dif_y2 - coordinates of the second vector.
- // Returns orientation test result for input vectors.
- template <typename T>
- static kOrientation orientation_test(T dif_x1_, T dif_y1_,
- T dif_x2_, T dif_y2_) {
- polygon_ulong_long_type dif_x1, dif_y1, dif_x2, dif_y2;
- bool dif_x1_plus, dif_x2_plus, dif_y1_plus, dif_y2_plus;
- dif_x1_plus = convert_to_65_bit(dif_x1_, dif_x1);
- dif_y1_plus = convert_to_65_bit(dif_y1_, dif_y1);
- dif_x2_plus = convert_to_65_bit(dif_x2_, dif_x2);
- dif_y2_plus = convert_to_65_bit(dif_y2_, dif_y2);
-
- polygon_ulong_long_type expr_l = dif_x1 * dif_y2;
- polygon_ulong_long_type expr_r = dif_x2 * dif_y1;
-
- bool expr_l_plus = (dif_x1_plus == dif_y2_plus) ? true : false;
- bool expr_r_plus = (dif_x2_plus == dif_y1_plus) ? true : false;
-
- if (expr_l == 0)
- expr_l_plus = true;
- if (expr_r == 0)
- expr_r_plus = true;
-
- if ((expr_l_plus == expr_r_plus) && (expr_l == expr_r))
- return COLLINEAR;
-
- if (!expr_l_plus) {
- if (expr_r_plus)
- return RIGHT_ORIENTATION;
- else
- return (expr_l > expr_r) ?
- RIGHT_ORIENTATION : LEFT_ORIENTATION;
- } else {
- if (!expr_r_plus)
- return LEFT_ORIENTATION;
- else
- return (expr_l < expr_r) ?
- RIGHT_ORIENTATION : LEFT_ORIENTATION;
- }
- }
-
- // Robust orientation test. Works correctly for any input coordinate type
- // that can be casted without lose of data to integer type within a range
- // [-2^31, 2^31 - 1] - signed integer type.
- // Arguments: point1, point2 - represent the first vector;
- // point2, point3 - represent the second vector;
- // Returns orientation test result for input vectors.
- template <typename T>
- static inline kOrientation orientation_test(const point_2d<T> &point1,
- const point_2d<T> &point2,
- const point_2d<T> &point3) {
- return orientation_test(point1.x() - point2.x(),
- point1.y() - point2.y(),
- point2.x() - point3.x(),
- point2.y() - point3.y());
- }
-
// Value is a determinant of two vectors.
// Return orientation based on the sign of the determinant.
template <typename T>
@@ -785,186 +699,34 @@
}
}
- // Class used to make computations with epsilon relative error.
- // ERC consists of two values: value1 and value2, both positive.
- // The resulting expression is equal to the value1 - value2.
- // The main idea is to represent any expression that consists of
- // addition, substraction, multiplication and division operations
- // to avoid using substraction. Substraction of a positive value
- // is equivalent to the addition to value2 and substraction of
- // a negative value is equivalent to the addition to value1.
- // Cons: ERC gives error relative not to the resulting value,
- // but relative to some expression instead. Example:
- // center_x = 100, ERC's value1 = 10^20, value2 = 10^20,
- // center_x = 1000, ERC's value3 = 10^21, value4 = 10^21,
- // such two centers are considered equal(
- // value1 + value4 = value2 + value3), while they are not.
- // Pros: ERC is much faster then approaches based on the use
- // of high-precision libraries. However this will give correct
- // answer for the previous example.
- // Solution: Use ERCs in case of defined comparison results and use
- // high-precision libraries for undefined results.
+ // Robust orientation test. Works correctly for any input type that
+ // can be casted without lose of data to the integer type within a range
+ // [-2^32, 2^32-1].
+ // Arguments: dif_x1_, dif_y1 - coordinates of the first vector.
+ // dif_x2_, dif_y2 - coordinates of the second vector.
+ // Returns orientation test result for input vectors.
template <typename T>
- class epsilon_robust_comparator {
- public:
- epsilon_robust_comparator() :
- positive_sum_(0),
- negative_sum_(0) {}
-
- epsilon_robust_comparator(T value) :
- positive_sum_((value>0)?value:0),
- negative_sum_((value<0)?-value:0) {}
-
- epsilon_robust_comparator(T pos, T neg) :
- positive_sum_(pos),
- negative_sum_(neg) {}
-
- T dif() const {
- return positive_sum_ - negative_sum_;
- }
-
- T pos() const {
- return positive_sum_;
- }
-
- T neg() const {
- return negative_sum_;
- }
-
- // Equivalent to the unary minus.
- void swap() {
- (std::swap)(positive_sum_, negative_sum_);
- }
-
- bool abs() {
- if (positive_sum_ < negative_sum_) {
- swap();
- return true;
- }
- return false;
- }
-
- epsilon_robust_comparator<T> &operator+=(const T &val) {
- if (val >= 0)
- positive_sum_ += val;
- else
- negative_sum_ -= val;
- return *this;
- }
-
- epsilon_robust_comparator<T> &operator+=(
- const epsilon_robust_comparator<T> &erc) {
- positive_sum_ += erc.positive_sum_;
- negative_sum_ += erc.negative_sum_;
- return *this;
- }
-
- epsilon_robust_comparator<T> &operator-=(const T &val) {
- if (val >= 0)
- negative_sum_ += val;
- else
- positive_sum_ -= val;
- return *this;
- }
-
- epsilon_robust_comparator<T> &operator-=(
- const epsilon_robust_comparator<T> &erc) {
- positive_sum_ += erc.negative_sum_;
- negative_sum_ += erc.positive_sum_;
- return *this;
- }
-
- epsilon_robust_comparator<T> &operator*=(const T &val) {
- positive_sum_ *= fabs(val);
- negative_sum_ *= fabs(val);
- if (val < 0) {
- swap();
- }
- return *this;
- }
-
- epsilon_robust_comparator<T> &operator/=(const T &val) {
- positive_sum_ /= fabs(val);
- negative_sum_ /= fabs(val);
- if (val < 0) {
- swap();
- }
- return *this;
- }
-
- // Compare predicate with value using ulp precision.
- kPredicateResult compare(T value, int ulp = 0) const {
- T lhs = positive_sum_ - ((value < 0) ? value : 0);
- T rhs = negative_sum_ + ((value > 0) ? value : 0);
- if (almost_equal(lhs, rhs, ulp))
- return UNDEFINED;
- return (lhs < rhs) ? LESS : MORE;
- }
-
- // Compare two predicats using ulp precision.
- kPredicateResult compare(const epsilon_robust_comparator &rc,
- int ulp = 0) const {
- T lhs = positive_sum_ + rc.neg();
- T rhs = negative_sum_ + rc.pos();
- if (almost_equal(lhs, rhs, ulp))
- return UNDEFINED;
- return (lhs < rhs) ? LESS : MORE;
- }
-
- // Check whether comparison has undefined result.
- bool compares_undefined(const epsilon_robust_comparator &rc,
- int ulp) const {
- return this->compare(rc, ulp) == UNDEFINED;
- }
-
- private:
- T positive_sum_;
- T negative_sum_;
- };
+ static kOrientation orientation_test(T dif_x1_, T dif_y1_,
+ T dif_x2_, T dif_y2_) {
+ return orientation_test(
+ robust_cross_product(dif_x1_, dif_y1_, dif_x2_, dif_y2_));
+ }
- template<typename T>
- inline epsilon_robust_comparator<T> operator+(
- const epsilon_robust_comparator<T>& lhs,
- const epsilon_robust_comparator<T>& rhs) {
- return epsilon_robust_comparator<T>(lhs.pos() + rhs.pos(),
- lhs.neg() + rhs.neg());
- }
-
- template<typename T>
- inline epsilon_robust_comparator<T> operator-(
- const epsilon_robust_comparator<T>& lhs,
- const epsilon_robust_comparator<T>& rhs) {
- return epsilon_robust_comparator<T>(lhs.pos() - rhs.neg(),
- lhs.neg() + rhs.pos());
- }
-
- template<typename T>
- inline epsilon_robust_comparator<T> operator*(
- const epsilon_robust_comparator<T>& lhs,
- const epsilon_robust_comparator<T>& rhs) {
- double res_pos = lhs.pos() * rhs.pos() + lhs.neg() * rhs.neg();
- double res_neg = lhs.pos() * rhs.neg() + lhs.neg() * rhs.pos();
- return epsilon_robust_comparator<T>(res_pos, res_neg);
- }
-
- template<typename T>
- inline epsilon_robust_comparator<T> operator*(
- const epsilon_robust_comparator<T>& lhs, const T& val) {
- if (val >= 0)
- return epsilon_robust_comparator<T>(lhs.pos() * val,
- lhs.neg() * val);
- return epsilon_robust_comparator<T>(-lhs.neg() * val,
- -lhs.pos() * val);
- }
-
- template<typename T>
- inline epsilon_robust_comparator<T> operator*(
- const T& val, const epsilon_robust_comparator<T>& rhs) {
- if (val >= 0)
- return epsilon_robust_comparator<T>(val * rhs.pos(),
- val * rhs.neg());
- return epsilon_robust_comparator<T>(-val * rhs.neg(),
- -val * rhs.pos());
+ // Robust orientation test. Works correctly for any input coordinate type
+ // that can be casted without lose of data to integer type within a range
+ // [-2^31, 2^31 - 1] - signed integer type.
+ // Arguments: point1, point2 - represent the first vector;
+ // point2, point3 - represent the second vector;
+ // Returns orientation test result for input vectors.
+ template <typename T>
+ static inline kOrientation orientation_test(const point_2d<T> &point1,
+ const point_2d<T> &point2,
+ const point_2d<T> &point3) {
+ return orientation_test(
+ robust_cross_product(point1.x() - point2.x(),
+ point1.y() - point2.y(),
+ point2.x() - point3.x(),
+ point2.y() - point3.y()));
}
// Robust voronoi vertex data structure. Used during removing degenerate
@@ -1034,7 +796,7 @@
double b1 = segment_end.y() - segment_start.y();
double a3 = new_point.x() - segment_start.x();
double b3 = new_point.y() - segment_start.y();
- double k = sqrt(a1 * a1 + b1 * b1);
+ double k = std::sqrt(a1 * a1 + b1 * b1);
// Avoid substraction while computing k.
if (segment.is_inverse()) {
if (b1 >= 0.0) {
@@ -1296,385 +1058,416 @@
return dif_x * dif_x + dif_y * dif_y;
}
- //// Recompute parameters of the circle event using high-precision library.
- //template <typename T>
- //static bool create_circle_event_ppp_gmpxx(const site_event<T> &site1,
- // const site_event<T> &site2,
- // const site_event<T> &site3,
- // circle_event<T> &c_event) {
- // typedef mpt_wrapper<mpz_class, 8> mpt_type;
- // static mpt_type mpz_dif_x[3], mpz_dif_y[3], mpz_sum_x[3], mpz_sum_y[3],
- // mpz_numerator[3], mpz_c_x, mpz_c_y, mpz_sqr_r;
- // mpz_dif_x[0] = site1.x() - site2.x();
- // mpz_dif_x[1] = site2.x() - site3.x();
- // mpz_dif_x[2] = site1.x() - site3.x();
- // mpz_dif_y[0] = site1.y() - site2.y();
- // mpz_dif_y[1] = site2.y() - site3.y();
- // mpz_dif_y[2] = site1.y() - site3.y();
- //
- // // Evaluate orientation.
- // double orientation = (mpz_dif_x[0] * mpz_dif_y[1] - mpz_dif_x[1] * mpz_dif_y[0]).get_d();
-
- // // If use this function only to recompute parameters of the circle
- // // event, this check won't be required.
- // if (orientation_test(orientation) != RIGHT_ORIENTATION)
- // return false;
-
- // // Evaluate inverse orientation to avoid division in calculations.
- // // r(inv_orientation) = 2 * EPS.
- // double inv_orientation = 0.5 / orientation;
- //
- // mpz_sum_x[0] = site1.x() + site2.x();
- // mpz_sum_x[1] = site2.x() + site3.x();
- // mpz_sum_y[0] = site1.y() + site2.y();
- // mpz_sum_y[1] = site2.y() + site3.y();
- //
- // mpz_numerator[1] = mpz_dif_x[0] * mpz_sum_x[0] + mpz_dif_y[0] * mpz_sum_y[0];
- // mpz_numerator[2] = mpz_dif_x[1] * mpz_sum_x[1] + mpz_dif_y[1] * mpz_sum_y[1];
-
- // mpz_c_x = mpz_numerator[1] * mpz_dif_y[1] - mpz_numerator[2] * mpz_dif_y[0];
- // mpz_c_y = mpz_numerator[2] * mpz_dif_x[0] - mpz_numerator[1] * mpz_dif_x[1];
-
- // // Evaluate radius of the circle.
- // mpz_sqr_r = 1.0;
- // for (int i = 0; i < 3; ++i)
- // mpz_sqr_r *= mpz_dif_x[i] * mpz_dif_x[i] + mpz_dif_y[i] * mpz_dif_y[i];
-
- // // Evaluate coordinates of the center of the circle.
- // // r(c_x) = r(c_y) = 4 * EPS.
- // double c_x = mpz_c_x.get_d() * inv_orientation;
- // double c_y = mpz_c_y.get_d() * inv_orientation;
-
- // // r(r) = 1.5 * EPS < 2 * EPS.
- // double r = sqrt(mpz_sqr_r.get_d());
-
- // // If c_x >= 0 then lower_x = c_x + r,
- // // else lower_x = (c_x * c_x - r * r) / (c_x - r).
- // // To guarantee epsilon relative error.
- // if (c_x >= 0) {
- // // r(lower_x) = 5 * EPS.
- // c_event = circle_event<double>(c_x, c_y, c_x + r * fabs(inv_orientation));
- // } else {
- // mpz_numerator[0] = mpz_c_x * mpz_c_x - mpz_sqr_r;
- // // r(lower_x) = 8 * EPS.
- // double lower_x = mpz_numerator[0].get_d() * fabs(inv_orientation);
- // lower_x /= (mpz_c_x >= 0) ? (-mpz_c_x.get_d() - r) : (mpz_c_x.get_d() - r);
- // c_event = circle_event<double>(c_x, c_y, lower_x);
- // }
- // return true;
- //}
-
- //// Recompute parameters of the circle event using high-precision library.
- //template <typename T>
- //static bool create_circle_event_pps_gmpxx(const site_event<T> &site1,
- // const site_event<T> &site2,
- // const site_event<T> &site3,
- // int segment_index,
- // circle_event<T> &c_event) {
- // typedef mpt_wrapper<mpz_class, 8> mpt_type;
- // // TODO(asydorchuk): Checks are not required if we recompute event parameters.
- // if (segment_index != 2) {
- // kOrientation orient1 = orientation_test(site1.point0(),
- // site2.point0(), site3.point0(true));
- // kOrientation orient2 = orientation_test(site1.point0(),
- // site2.point0(), site3.point1(true));
- // if (segment_index == 1 && site1.x0() >= site2.x0()) {
- // if (orient1 != RIGHT_ORIENTATION)
- // return false;
- // } else if (segment_index == 3 && site2.x0() >= site1.x0()) {
- // if (orient2 != RIGHT_ORIENTATION)
- // return false;
- // } else if (orient1 != RIGHT_ORIENTATION && orient2 != RIGHT_ORIENTATION) {
- // return false;
- // }
- // } else {
- // if (site3.point0(true) == site1.point0() &&
- // site3.point1(true) == site2.point0())
- // return false;
- // }
-
- // static mpt_type line_a, line_b, segm_len, vec_x, vec_y, sum_x, sum_y, teta,
- // denom, A, B, sum_AB, dif[2], numer, cA[4], cB[4], det;
-
- // line_a = site3.point1(true).y() - site3.point0(true).y();
- // line_b = site3.point0(true).x() - site3.point1(true).x();
- // segm_len = line_a * line_a + line_b * line_b;
- // vec_x = site2.y() - site1.y();
- // vec_y = site1.x() - site2.x();
- // sum_x = site1.x() + site2.x();
- // sum_y = site1.y() + site2.y();
- // teta = line_a * vec_x + line_b * vec_y;
- // denom = vec_x * line_b - vec_y * line_a;
- //
- // dif[0] = site3.point1().y() - site1.y();
- // dif[1] = site1.x() - site3.point1().x();
- // A = line_a * dif[1] - line_b * dif[0];
- // dif[0] = site3.point1().y() - site2.y();
- // dif[1] = site2.x() - site3.point1().x();
- // B = line_a * dif[1] - line_b * dif[0];
- // sum_AB = A + B;
-
- // if (denom == 0) {
- // numer = teta * teta - sum_AB * sum_AB;
- // denom = teta * sum_AB;
- // cA[0] = denom * sum_x * 2 + numer * vec_x;
- // cB[0] = segm_len;
- // cA[1] = denom * sum_AB * 2 + numer * teta;
- // cB[1] = 1;
- // cA[2] = denom * sum_y * 2 + numer * vec_y;
- // double c_x = 0.25 * cA[0].get_d() / denom.get_d();
- // double c_y = 0.25 * cA[2].get_d() / denom.get_d();
- // double lower_x = 0.25 * sqr_expr_evaluator<2>::eval(cA, cB) /
- // (denom.get_d() * sqrt(segm_len.get_d()));
- // c_event = circle_event<double>(c_x, c_y, lower_x);
- // return true;
- // }
-
- // det = (teta * teta + denom * denom) * A * B * 4;
- // cA[0] = sum_x * denom * denom + teta * sum_AB * vec_x;
- // cB[0] = 1;
- // cA[1] = (segment_index == 2) ? -vec_x : vec_x;
- // cB[1] = det;
- // double c_x = 0.5 * sqr_expr_evaluator<2>::eval(cA, cB) / (denom.get_d() * denom.get_d());
- //
- // cA[2] = sum_y * denom * denom + teta * sum_AB * vec_y;
- // cB[2] = 1;
- // cA[3] = (segment_index == 2) ? -vec_y : vec_y;
- // cB[3] = det;
- // double c_y = 0.5 * sqr_expr_evaluator<2>::eval(&cA[2], &cB[2]) / (denom.get_d() * denom.get_d());
- //
- // cB[0] *= segm_len;
- // cB[1] *= segm_len;
- // cA[2] = sum_AB * (denom * denom + teta * teta);
- // cB[2] = 1;
- // cA[3] = (segment_index == 2) ? -teta : teta;
- // cB[3] = det;
- // double lower_x = 0.5 * sqr_expr_evaluator<4>::eval(cA, cB) /
- // (denom.get_d() * denom.get_d() * sqrt(segm_len.get_d()));
- //
- // c_event = circle_event<double>(c_x, c_y, lower_x);
- // return true;
- //}
-
- //// Evaluates A[3] + A[0] * sqrt(B[0]) + A[1] * sqrt(B[1]) +
- //// A[2] * sqrt(B[3] * (sqrt(B[0] * B[1]) + B[2]));
- //template <typename mpt>
- //static double sqr_expr_evaluator_pss(mpt *A, mpt *B) {
- // static mpt cA[4], cB[4];
- // if (A[3] == 0) {
- // double lh = sqr_expr_evaluator<2>::eval<mpt>(A, B);
- // cA[0] = 1;
- // cB[0] = B[0] * B[1];
- // cA[1] = B[2];
- // cB[1] = 1;
- // double rh = A[2].get_d() * sqrt(B[3].get_d() *
- // sqr_expr_evaluator<2>::eval<mpt>(cA, cB));
- // if (((lh >= 0.0) && (rh >= 0.0) || (lh <= 0.0) && (rh <= 0.0))) {
- // return lh + rh;
- // }
- // cA[0] = A[0] * A[0] * B[0] + A[1] * A[1] * B[1];
- // cA[0] -= A[2] * A[2] * B[3] * B[2];
- // cB[0] = 1;
- // cA[1] = A[0] * A[1] * 2 - A[2] * A[2] * B[3];
- // cB[1] = B[0] * B[1];
- // return sqr_expr_evaluator<2>::eval(cA, cB) / (lh - rh);
- // }
- // cA[0] = A[3];
- // cB[0] = 1;
- // cA[1] = A[0];
- // cB[1] = B[0];
- // cA[2] = A[1];
- // cB[2] = B[1];
- // double lh = sqr_expr_evaluator<3>::eval<mpt>(cA, cB);
- // cA[0] = 1;
- // cB[0] = B[0] * B[1];
- // cA[1] = B[2];
- // cB[1] = 1;
- // double rh = A[2].get_d() * sqrt(B[3].get_d() *
- // sqr_expr_evaluator<2>::eval<mpt>(cA, cB));
- // if ((lh >= 0.0) && (rh >= 0.0) || (lh <= 0.0) && (rh <= 0.0)) {
- // return lh + rh;
- // }
- // cA[0] = A[0] * A[0] * B[0] + A[1] * A[1] * B[1];
- // cA[0] += A[3] * A[3] - A[2] * A[2] * B[2] * B[3];
- // cB[0] = 1;
- // cA[1] = A[3] * A[0] * 2;
- // cB[1] = B[0];
- // cA[2] = A[3] * A[1] * 2;
- // cB[2] = B[1];
- // cA[3] = A[0] * A[1] * 2 - A[2] * A[2] * B[3];
- // cB[3] = B[0] * B[1];
- // return sqr_expr_evaluator<4>::eval(cA, cB) / (lh - rh);
- //}
-
- //// Recompute parameters of the circle event using high-precision library.
- //template <typename T>
- //static bool create_circle_event_pss_gmpxx(const site_event<T> &site1,
- // const site_event<T> &site2,
- // const site_event<T> &site3,
- // int point_index,
- // circle_event<T> &c_event) {
- // typedef mpt_wrapper<mpz_class, 8> mpt_type;
- // if (site2.index() == site3.index()) {
- // return false;
- // }
-
- // const point_2d<T> &segm_start1 = site2.point1(true);
- // const point_2d<T> &segm_end1 = site2.point0(true);
- // const point_2d<T> &segm_start2 = site3.point0(true);
- // const point_2d<T> &segm_end2 = site3.point1(true);
-
- // if (point_index == 2) {
- // if (!site2.is_inverse() && site3.is_inverse())
- // return false;
- // if (site2.is_inverse() == site3.is_inverse() &&
- // orientation_test(segm_end1, site1.point0(), segm_end2) != RIGHT_ORIENTATION)
- // return false;
- // }
-
- // static mpt_type a[2], b[2], c[2], cA[4], cB[4], or, dx, dy, ix, iy;
- // a[0] = segm_end1.x() - segm_start1.x();
- // b[0] = segm_end1.y() - segm_start1.y();
- // a[1] = segm_end2.x() - segm_start2.x();
- // b[1] = segm_end2.y() - segm_start2.y();
- // or = a[1] * b[0] - a[0] * b[1];
- // if (or == 0) {
- // double denom = (a[0] * a[0] + b[0] * b[0]).get_d() * 2.0;
-
- // c[0] = b[0] * (segm_start2.x() - segm_start1.x()) -
- // a[0] * (segm_start2.y() - segm_start1.y());
- // dx = a[0] * (site1.y() - segm_start1.y()) -
- // b[0] * (site1.x() - segm_start1.x());
- // dy = b[0] * (site1.x() - segm_start2.x()) -
- // a[0] * (site1.y() - segm_start2.y());
- // cA[0] = b[0] * ((point_index == 2) ? 2 : -2);
- // cB[0] = dx * dy;
- // cA[1] = a[0] * a[0] * (segm_start1.y() + segm_start2.y()) -
- // a[0] * b[0] * (segm_start1.x() + segm_start2.x() - 2.0 * site1.x()) +
- // b[0] * b[0] * (2.0 * site1.y());
- // cB[1] = 1;
- // double c_y = sqr_expr_evaluator<2>::eval(cA, cB);
-
- // cA[0] = a[0] * ((point_index == 2) ? 2 : -2);
- // cA[1] = b[0] * b[0] * (segm_start1.x() + segm_start2.x()) -
- // a[0] * b[0] * (segm_start1.y() + segm_start2.y() - 2.0 * site1.y()) +
- // a[0] * a[0] * (2.0 * site1.x());
- // double c_x = sqr_expr_evaluator<2>::eval(cA, cB);
-
- // cA[2] = (c[0] < 0) ? -c[0] : c[0];
- // cB[2] = a[0] * a[0] + b[0] * b[0];
- // double lower_x = sqr_expr_evaluator<3>::eval(cA, cB);
- // c_event = circle_event<T>(c_x / denom, c_y / denom, lower_x / denom);
- // return true;
- // }
- // c[0] = b[0] * segm_end1.x() - a[0] * segm_end1.y();
- // c[1] = a[1] * segm_end2.y() - b[1] * segm_end2.x();
- // ix = a[0] * c[1] + a[1] * c[0];
- // iy = b[0] * c[1] + b[1] * c[0];
- // dx = ix - or * site1.x();
- // dy = iy - or * site1.y();
- // if ((dx == 0) && (dy == 0)) {
- // double denom = or.get_d();
- // double c_x = ix.get_d() / denom;
- // double c_y = iy.get_d() / denom;
- // c_event = circle_event<T>(c_x, c_y, c_x);
- // return true;
- // }
-
- // double sign = ((point_index == 2) ? 1 : -1) * ((or < 0) ? 1 : -1);
- // cA[0] = a[1] * -dx + b[1] * -dy;
- // cA[1] = a[0] * -dx + b[0] * -dy;
- // cA[2] = sign;
- // cA[3] = 0;
- // cB[0] = a[0] * a[0] + b[0] * b[0];
- // cB[1] = a[1] * a[1] + b[1] * b[1];
- // cB[2] = a[0] * a[1] + b[0] * b[1];
- // cB[3] = (a[0] * dy - b[0] * dx) * (a[1] * dy - b[1] * dx) * -2;
- // double denom = sqr_expr_evaluator_pss<mpt_type>(cA, cB);
-
- // cA[0] = b[1] * (dx * dx + dy * dy) - iy * (dx * a[1] + dy * b[1]);
- // cA[1] = b[0] * (dx * dx + dy * dy) - iy * (dx * a[0] + dy * b[0]);
- // cA[2] = iy * sign;
- // double cy = sqr_expr_evaluator_pss<mpt_type>(cA, cB);
-
- // cA[0] = a[1] * (dx * dx + dy * dy) - ix * (dx * a[1] + dy * b[1]);
- // cA[1] = a[0] * (dx * dx + dy * dy) - ix * (dx * a[0] + dy * b[0]);
- // cA[2] = ix * sign;
- // double cx = sqr_expr_evaluator_pss<mpt_type>(cA, cB);
-
- // cA[3] = or * (dx * dx + dy * dy) * (denom < 0 ? -1 : 1);
- // double lower_x = sqr_expr_evaluator_pss<mpt_type>(cA, cB);
- // denom *= or.get_d();
- // c_event = circle_event<T>(cx / denom, cy / denom, lower_x / denom);
- // return true;
- //}
-
- //template <typename T>
- //static mpt_wrapper<mpz_class, 8> &mpt_cross(T a0, T b0, T a1, T b1) {
- // static mpt_wrapper<mpz_class, 8> temp[2];
- // temp[0] = a0;
- // temp[1] = b0;
- // temp[0] = temp[0] * b1;
- // temp[1] = temp[1] * a1;
- // temp[0] -= temp[1];
- // return temp[0];
- //}
-
- //// Recompute parameters of the circle event using high-precision library.
- //template <typename T>
- //static bool create_circle_event_sss_gmpxx(const site_event<T> &site1,
- // const site_event<T> &site2,
- // const site_event<T> &site3,
- // circle_event<T> &c_event) {
- // typedef mpt_wrapper<mpz_class, 8> mpt_type;
- // if (site1.index() == site2.index() ||
- // site2.index() == site3.index())
- // return false;
- // static mpt_type a[3], b[3], c[3], sqr_len[4], cross[4];
- // a[0] = site1.x1(true) - site1.x0(true);
- // a[1] = site2.x1(true) - site2.x0(true);
- // a[2] = site3.x1(true) - site3.x0(true);
- //
- // b[0] = site1.y1(true) - site1.y0(true);
- // b[1] = site2.y1(true) - site2.y0(true);
- // b[2] = site3.y1(true) - site3.y0(true);
-
- // c[0] = mpt_cross(site1.x0(true), site1.y0(true), site1.x1(true), site1.y1(true));
- // c[1] = mpt_cross(site2.x0(true), site2.y0(true), site2.x1(true), site2.y1(true));
- // c[2] = mpt_cross(site3.x0(true), site3.y0(true), site3.x1(true), site3.y1(true));
-
- // for (int i = 0; i < 3; ++i) {
- // int j = (i+1) % 3;
- // int k = (i+2) % 3;
- // cross[i] = a[j] * b[k] - a[k] * b[j];
- // sqr_len[i] = a[i] * a[i] + b[i] * b[i];
- // }
- // double denom = sqr_expr_evaluator<3>::eval(cross, sqr_len);
-
- // for (int i = 0; i < 3; ++i) {
- // int j = (i+1) % 3;
- // int k = (i+2) % 3;
- // cross[i] = b[j] * c[k] - b[k] * c[j];
- // sqr_len[i] = a[i] * a[i] + b[i] * b[i];
- // }
- // double c_y = sqr_expr_evaluator<3>::eval(cross, sqr_len) / denom;
-
- // cross[3] = 0;
- // for (int i = 0; i < 3; ++i) {
- // int j = (i+1) % 3;
- // int k = (i+2) % 3;
- // cross[i] = a[j] * c[k] - a[k] * c[j];
- // sqr_len[i] = a[i] * a[i] + b[i] * b[i];
- // cross[3] += cross[i] * b[i];
- // }
- // double c_x = sqr_expr_evaluator<3>::eval(cross, sqr_len) / denom;
- //
- // sqr_len[3] = 1;
- // double lower_x = sqr_expr_evaluator<4>::eval(cross, sqr_len) / denom;
- //
- // c_event = circle_event<double>(c_x, c_y, lower_x);
- // return true;
- //}
+ // Recompute parameters of the circle event using high-precision library.
+ template <typename T>
+ static bool create_circle_event_ppp_gmpxx(const site_event<T> &site1,
+ const site_event<T> &site2,
+ const site_event<T> &site3,
+ circle_event<T> &c_event,
+ bool recompute_c_x = true,
+ bool recompute_c_y = true,
+ bool recompute_lower_x = true) {
+ typedef mpt_wrapper<mpz_class, 8> mpt_type;
+ static mpt_type mpz_dif_x[3], mpz_dif_y[3], mpz_sum_x[3], mpz_sum_y[3],
+ mpz_numerator[3], mpz_c_x, mpz_c_y, mpz_sqr_r, denom,
+ cA[2], cB[2];
+ mpz_dif_x[0] = site1.x() - site2.x();
+ mpz_dif_x[1] = site2.x() - site3.x();
+ mpz_dif_x[2] = site1.x() - site3.x();
+ mpz_dif_y[0] = site1.y() - site2.y();
+ mpz_dif_y[1] = site2.y() - site3.y();
+ mpz_dif_y[2] = site1.y() - site3.y();
+
+ // Evaluate orientation.
+ denom = (mpz_dif_x[0] * mpz_dif_y[1] - mpz_dif_x[1] * mpz_dif_y[0]) * 2.0;
+
+ // Evaluate inverse orientation to avoid division in calculations.
+ // r(inv_orientation) = 2 * EPS.
+ if (denom.get_d() >= 0)
+ return false;
+
+ mpz_sum_x[0] = site1.x() + site2.x();
+ mpz_sum_x[1] = site2.x() + site3.x();
+ mpz_sum_y[0] = site1.y() + site2.y();
+ mpz_sum_y[1] = site2.y() + site3.y();
+
+ mpz_numerator[1] = mpz_dif_x[0] * mpz_sum_x[0] + mpz_dif_y[0] * mpz_sum_y[0];
+ mpz_numerator[2] = mpz_dif_x[1] * mpz_sum_x[1] + mpz_dif_y[1] * mpz_sum_y[1];
+
+ if (recompute_c_x || recompute_lower_x) {
+ mpz_c_x = mpz_numerator[1] * mpz_dif_y[1] - mpz_numerator[2] * mpz_dif_y[0];
+ c_event.x(mpz_c_x.get_d() / denom.get_d());
+
+ if (recompute_lower_x) {
+ // Evaluate radius of the circle.
+ mpz_sqr_r = 1.0;
+ for (int i = 0; i < 3; ++i)
+ mpz_sqr_r *= mpz_dif_x[i] * mpz_dif_x[i] + mpz_dif_y[i] * mpz_dif_y[i];
+
+ // r(r) = 1.5 * EPS < 2 * EPS.
+ double r = std::sqrt(mpz_sqr_r.get_d());
+
+ // If c_x >= 0 then lower_x = c_x + r,
+ // else lower_x = (c_x * c_x - r * r) / (c_x - r).
+ // To guarantee epsilon relative error.
+ if (c_event.x() >= 0) {
+ // r(lower_x) = 5 * EPS.
+ c_event.lower_x(c_event.x() + r / fabs(denom.get_d()));
+ } else {
+ mpz_numerator[0] = mpz_c_x * mpz_c_x - mpz_sqr_r;
+ // r(lower_x) = 8 * EPS.
+ double lower_x = mpz_numerator[0].get_d() / fabs(denom.get_d());
+ lower_x /= (mpz_c_x >= 0) ? (-mpz_c_x.get_d() - r) : (mpz_c_x.get_d() - r);
+ c_event.lower_x(lower_x);
+ }
+ }
+ }
+
+ if (recompute_c_y) {
+ mpz_c_y = mpz_numerator[2] * mpz_dif_x[0] - mpz_numerator[1] * mpz_dif_x[1];
+ c_event.y(mpz_c_y.get_d() / denom.get_d());
+ }
+
+ return true;
+ }
+
+ // Recompute parameters of the circle event using high-precision library.
+ template <typename T>
+ static bool create_circle_event_pps_gmpxx(const site_event<T> &site1,
+ const site_event<T> &site2,
+ const site_event<T> &site3,
+ int segment_index,
+ circle_event<T> &c_event) {
+ typedef mpt_wrapper<mpz_class, 8> mpt_type;
+ typedef mpt_wrapper<mpf_class, 2> mpf_type;
+ // TODO(asydorchuk): Checks are not required if we recompute event parameters.
+ if (segment_index != 2) {
+ kOrientation orient1 = orientation_test(site1.point0(),
+ site2.point0(), site3.point0(true));
+ kOrientation orient2 = orientation_test(site1.point0(),
+ site2.point0(), site3.point1(true));
+ if (segment_index == 1 && site1.x0() >= site2.x0()) {
+ if (orient1 != RIGHT_ORIENTATION)
+ return false;
+ } else if (segment_index == 3 && site2.x0() >= site1.x0()) {
+ if (orient2 != RIGHT_ORIENTATION)
+ return false;
+ } else if (orient1 != RIGHT_ORIENTATION && orient2 != RIGHT_ORIENTATION) {
+ return false;
+ }
+ } else {
+ if (site3.point0(true) == site1.point0() &&
+ site3.point1(true) == site2.point0())
+ return false;
+ }
+
+ static mpt_type line_a, line_b, segm_len, vec_x, vec_y, sum_x, sum_y, teta,
+ denom, A, B, sum_AB, dif[2], numer, cA[4], cB[4], det;
+
+ line_a = site3.point1(true).y() - site3.point0(true).y();
+ line_b = site3.point0(true).x() - site3.point1(true).x();
+ segm_len = line_a * line_a + line_b * line_b;
+ vec_x = site2.y() - site1.y();
+ vec_y = site1.x() - site2.x();
+ sum_x = site1.x() + site2.x();
+ sum_y = site1.y() + site2.y();
+ teta = line_a * vec_x + line_b * vec_y;
+ denom = vec_x * line_b - vec_y * line_a;
+
+ dif[0] = site3.point1().y() - site1.y();
+ dif[1] = site1.x() - site3.point1().x();
+ A = line_a * dif[1] - line_b * dif[0];
+ dif[0] = site3.point1().y() - site2.y();
+ dif[1] = site2.x() - site3.point1().x();
+ B = line_a * dif[1] - line_b * dif[0];
+ sum_AB = A + B;
+
+ if (denom == 0) {
+ numer = teta * teta - sum_AB * sum_AB;
+ denom = teta * sum_AB;
+ cA[0] = denom * sum_x * 2 + numer * vec_x;
+ cB[0] = segm_len;
+ cA[1] = denom * sum_AB * 2 + numer * teta;
+ cB[1] = 1;
+ cA[2] = denom * sum_y * 2 + numer * vec_y;
+ double c_x = 0.25 * cA[0].get_d() / denom.get_d();
+ double c_y = 0.25 * cA[2].get_d() / denom.get_d();
+ double lower_x = 0.25 * sqr_expr_evaluator<2>::eval<mpt_type, mpf_type>(cA, cB).get_d() /
+ (denom.get_d() * std::sqrt(segm_len.get_d()));
+ c_event = circle_event<double>(c_x, c_y, lower_x);
+ return true;
+ }
+
+ det = (teta * teta + denom * denom) * A * B * 4;
+ cA[0] = sum_x * denom * denom + teta * sum_AB * vec_x;
+ cB[0] = 1;
+ cA[1] = (segment_index == 2) ? -vec_x : vec_x;
+ cB[1] = det;
+ double c_x = 0.5 * sqr_expr_evaluator<2>::eval<mpt_type, mpf_type>(cA, cB).get_d() /
+ (denom.get_d() * denom.get_d());
+
+ cA[2] = sum_y * denom * denom + teta * sum_AB * vec_y;
+ cB[2] = 1;
+ cA[3] = (segment_index == 2) ? -vec_y : vec_y;
+ cB[3] = det;
+ double c_y = 0.5 * sqr_expr_evaluator<2>::eval<mpt_type, mpf_type>(&cA[2], &cB[2]).get_d() /
+ (denom.get_d() * denom.get_d());
+
+ cB[0] *= segm_len;
+ cB[1] *= segm_len;
+ cA[2] = sum_AB * (denom * denom + teta * teta);
+ cB[2] = 1;
+ cA[3] = (segment_index == 2) ? -teta : teta;
+ cB[3] = det;
+ double lower_x = 0.5 * sqr_expr_evaluator<4>::eval<mpt_type, mpf_type>(cA, cB).get_d() /
+ (denom.get_d() * denom.get_d() * std::sqrt(segm_len.get_d()));
+
+ c_event = circle_event<double>(c_x, c_y, lower_x);
+ return true;
+ }
+
+ // Evaluates A[3] + A[0] * sqrt(B[0]) + A[1] * sqrt(B[1]) +
+ // A[2] * sqrt(B[3] * (sqrt(B[0] * B[1]) + B[2]));
+ template <typename mpt, typename mpf>
+ static mpf sqr_expr_evaluator_pss(mpt *A, mpt *B) {
+ static mpt cA[4], cB[4];
+ static mpf lh, rh, numer;
+ if (A[3] == 0) {
+ lh = sqr_expr_evaluator<2>::eval<mpt, mpf>(A, B);
+ cA[0] = 1;
+ cB[0] = B[0] * B[1];
+ cA[1] = B[2];
+ cB[1] = 1;
+ rh = A[2].get_d() * std::sqrt(B[3].get_d() *
+ sqr_expr_evaluator<2>::eval<mpt, mpf>(cA, cB).get_d());
+ if (((lh >= 0) && (rh >= 0)) || ((lh <= 0) && (rh <= 0))) {
+ return lh + rh;
+ }
+ cA[0] = A[0] * A[0] * B[0] + A[1] * A[1] * B[1];
+ cA[0] -= A[2] * A[2] * B[3] * B[2];
+ cB[0] = 1;
+ cA[1] = A[0] * A[1] * 2 - A[2] * A[2] * B[3];
+ cB[1] = B[0] * B[1];
+ numer = sqr_expr_evaluator<2>::eval<mpt, mpf>(cA, cB);
+ return numer / (lh - rh);
+ }
+ cA[0] = A[3];
+ cB[0] = 1;
+ cA[1] = A[0];
+ cB[1] = B[0];
+ cA[2] = A[1];
+ cB[2] = B[1];
+ lh = sqr_expr_evaluator<3>::eval<mpt, mpf>(cA, cB);
+ cA[0] = 1;
+ cB[0] = B[0] * B[1];
+ cA[1] = B[2];
+ cB[1] = 1;
+ rh = A[2].get_d() * std::sqrt(B[3].get_d() *
+ sqr_expr_evaluator<2>::eval<mpt, mpf>(cA, cB).get_d());
+ if (((lh >= 0) && (rh >= 0)) || ((lh <= 0) && (rh <= 0))) {
+ return lh + rh;
+ }
+ cA[0] = A[0] * A[0] * B[0] + A[1] * A[1] * B[1];
+ cA[0] += A[3] * A[3] - A[2] * A[2] * B[2] * B[3];
+ cB[0] = 1;
+ cA[1] = A[3] * A[0] * 2;
+ cB[1] = B[0];
+ cA[2] = A[3] * A[1] * 2;
+ cB[2] = B[1];
+ cA[3] = A[0] * A[1] * 2 - A[2] * A[2] * B[3];
+ cB[3] = B[0] * B[1];
+ numer = sqr_expr_evaluator<4>::eval<mpt, mpf>(cA, cB);
+ return numer / (lh - rh);
+ }
+
+ // Recompute parameters of the circle event using high-precision library.
+ template <typename T>
+ static bool create_circle_event_pss_gmpxx(const site_event<T> &site1,
+ const site_event<T> &site2,
+ const site_event<T> &site3,
+ int point_index,
+ circle_event<T> &c_event) {
+ typedef mpt_wrapper<mpz_class, 8> mpt_type;
+ typedef mpt_wrapper<mpf_class, 2> mpf_type;
+ if (site2.index() == site3.index()) {
+ return false;
+ }
+
+ const point_2d<T> &segm_start1 = site2.point1(true);
+ const point_2d<T> &segm_end1 = site2.point0(true);
+ const point_2d<T> &segm_start2 = site3.point0(true);
+ const point_2d<T> &segm_end2 = site3.point1(true);
+
+ if (point_index == 2) {
+ if (!site2.is_inverse() && site3.is_inverse())
+ return false;
+ if (site2.is_inverse() == site3.is_inverse() &&
+ orientation_test(segm_end1, site1.point0(), segm_end2) != RIGHT_ORIENTATION)
+ return false;
+ }
+
+ static mpt_type a[2], b[2], c[2], cA[4], cB[4], orientation, dx, dy, ix, iy;
+ a[0] = segm_end1.x() - segm_start1.x();
+ b[0] = segm_end1.y() - segm_start1.y();
+ a[1] = segm_end2.x() - segm_start2.x();
+ b[1] = segm_end2.y() - segm_start2.y();
+ orientation = a[1] * b[0] - a[0] * b[1];
+ if (orientation == 0) {
+ double denom = (a[0] * a[0] + b[0] * b[0]).get_d() * 2.0;
+
+ c[0] = b[0] * (segm_start2.x() - segm_start1.x()) -
+ a[0] * (segm_start2.y() - segm_start1.y());
+ dx = a[0] * (site1.y() - segm_start1.y()) -
+ b[0] * (site1.x() - segm_start1.x());
+ dy = b[0] * (site1.x() - segm_start2.x()) -
+ a[0] * (site1.y() - segm_start2.y());
+ cA[0] = b[0] * ((point_index == 2) ? 2 : -2);
+ cB[0] = dx * dy;
+ cA[1] = a[0] * a[0] * (segm_start1.y() + segm_start2.y()) -
+ a[0] * b[0] * (segm_start1.x() + segm_start2.x() - 2.0 * site1.x()) +
+ b[0] * b[0] * (2.0 * site1.y());
+ cB[1] = 1;
+ double c_y = sqr_expr_evaluator<2>::eval<mpt_type, mpf_type>(cA, cB).get_d();
+
+ cA[0] = a[0] * ((point_index == 2) ? 2 : -2);
+ cA[1] = b[0] * b[0] * (segm_start1.x() + segm_start2.x()) -
+ a[0] * b[0] * (segm_start1.y() + segm_start2.y() - 2.0 * site1.y()) +
+ a[0] * a[0] * (2.0 * site1.x());
+ double c_x = sqr_expr_evaluator<2>::eval<mpt_type, mpf_type>(cA, cB).get_d();
+
+ cA[2] = (c[0] < 0) ? -c[0] : c[0];
+ cB[2] = a[0] * a[0] + b[0] * b[0];
+ double lower_x = sqr_expr_evaluator<3>::eval<mpt_type, mpf_type>(cA, cB).get_d();
+ c_event = circle_event<T>(c_x / denom, c_y / denom, lower_x / denom);
+ return true;
+ }
+ c[0] = b[0] * segm_end1.x() - a[0] * segm_end1.y();
+ c[1] = a[1] * segm_end2.y() - b[1] * segm_end2.x();
+ ix = a[0] * c[1] + a[1] * c[0];
+ iy = b[0] * c[1] + b[1] * c[0];
+ dx = ix - orientation * site1.x();
+ dy = iy - orientation * site1.y();
+ if ((dx == 0) && (dy == 0)) {
+ double denom = orientation.get_d();
+ double c_x = ix.get_d() / denom;
+ double c_y = iy.get_d() / denom;
+ c_event = circle_event<T>(c_x, c_y, c_x);
+ return true;
+ }
+
+ double sign = ((point_index == 2) ? 1 : -1) * ((orientation < 0) ? 1 : -1);
+ cA[0] = a[1] * -dx + b[1] * -dy;
+ cA[1] = a[0] * -dx + b[0] * -dy;
+ cA[2] = sign;
+ cA[3] = 0;
+ cB[0] = a[0] * a[0] + b[0] * b[0];
+ cB[1] = a[1] * a[1] + b[1] * b[1];
+ cB[2] = a[0] * a[1] + b[0] * b[1];
+ cB[3] = (a[0] * dy - b[0] * dx) * (a[1] * dy - b[1] * dx) * -2;
+ double denom = sqr_expr_evaluator_pss<mpt_type, mpf_type>(cA, cB).get_d();
+
+ cA[0] = b[1] * (dx * dx + dy * dy) - iy * (dx * a[1] + dy * b[1]);
+ cA[1] = b[0] * (dx * dx + dy * dy) - iy * (dx * a[0] + dy * b[0]);
+ cA[2] = iy * sign;
+ double cy = sqr_expr_evaluator_pss<mpt_type, mpf_type>(cA, cB).get_d();
+
+ cA[0] = a[1] * (dx * dx + dy * dy) - ix * (dx * a[1] + dy * b[1]);
+ cA[1] = a[0] * (dx * dx + dy * dy) - ix * (dx * a[0] + dy * b[0]);
+ cA[2] = ix * sign;
+ double cx = sqr_expr_evaluator_pss<mpt_type, mpf_type>(cA, cB).get_d();
+
+ cA[3] = orientation * (dx * dx + dy * dy) * (denom < 0 ? -1 : 1);
+ double lower_x = sqr_expr_evaluator_pss<mpt_type, mpf_type>(cA, cB).get_d();
+ denom *= orientation.get_d();
+ c_event = circle_event<T>(cx / denom, cy / denom, lower_x / denom);
+ return true;
+ }
+
+ template <typename T>
+ static mpt_wrapper<mpz_class, 8> &mpt_cross(T a0, T b0, T a1, T b1) {
+ static mpt_wrapper<mpz_class, 8> temp[2];
+ temp[0] = a0;
+ temp[1] = b0;
+ temp[0] = temp[0] * b1;
+ temp[1] = temp[1] * a1;
+ temp[0] -= temp[1];
+ return temp[0];
+ }
+
+ // Recompute parameters of the circle event using high-precision library.
+ template <typename T>
+ static bool create_circle_event_sss_gmpxx(const site_event<T> &site1,
+ const site_event<T> &site2,
+ const site_event<T> &site3,
+ circle_event<T> &c_event,
+ bool recompute_c_x = true,
+ bool recompute_c_y = true,
+ bool recompute_lower_x = true,
+ bool recompute_denom = true) {
+ typedef mpt_wrapper<mpz_class, 8> mpt_type;
+ typedef mpt_wrapper<mpf_class, 2> mpf_type;
+ if (site1.index() == site2.index() ||
+ site2.index() == site3.index())
+ return false;
+ static mpt_type a[3], b[3], c[3], sqr_len[4], cross[4];
+ a[0] = site1.x1(true) - site1.x0(true);
+ a[1] = site2.x1(true) - site2.x0(true);
+ a[2] = site3.x1(true) - site3.x0(true);
+
+ b[0] = site1.y1(true) - site1.y0(true);
+ b[1] = site2.y1(true) - site2.y0(true);
+ b[2] = site3.y1(true) - site3.y0(true);
+
+ c[0] = mpt_cross(site1.x0(true), site1.y0(true), site1.x1(true), site1.y1(true));
+ c[1] = mpt_cross(site2.x0(true), site2.y0(true), site2.x1(true), site2.y1(true));
+ c[2] = mpt_cross(site3.x0(true), site3.y0(true), site3.x1(true), site3.y1(true));
+
+ for (int i = 0; i < 3; ++i) {
+ sqr_len[i] = a[i] * a[i] + b[i] * b[i];
+ }
+
+ for (int i = 0; i < 3; ++i) {
+ int j = (i+1) % 3;
+ int k = (i+2) % 3;
+ cross[i] = a[j] * b[k] - a[k] * b[j];
+ }
+ double denom = sqr_expr_evaluator<3>::eval<mpt_type, mpf_type>(cross, sqr_len).get_d();
+
+ if (recompute_c_y) {
+ for (int i = 0; i < 3; ++i) {
+ int j = (i+1) % 3;
+ int k = (i+2) % 3;
+ cross[i] = b[j] * c[k] - b[k] * c[j];
+ }
+ double c_y = sqr_expr_evaluator<3>::eval<mpt_type, mpf_type>(cross, sqr_len).get_d();
+ c_event.y(c_y / denom);
+ }
+
+ if (recompute_c_x || recompute_lower_x) {
+ cross[3] = 0;
+ for (int i = 0; i < 3; ++i) {
+ int j = (i+1) % 3;
+ int k = (i+2) % 3;
+ cross[i] = a[j] * c[k] - a[k] * c[j];
+ if (recompute_lower_x) {
+ cross[3] += cross[i] * b[i];
+ }
+ }
+
+ if (recompute_c_x) {
+ double c_x = sqr_expr_evaluator<3>::eval<mpt_type, mpf_type>(cross, sqr_len).get_d();
+ c_event.x(c_x / denom);
+ }
+
+ if (recompute_lower_x) {
+ sqr_len[3] = 1;
+ double lower_x = sqr_expr_evaluator<4>::eval<mpt_type, mpf_type>(cross, sqr_len).get_d();
+ c_event.lower_x(lower_x / denom);
+ }
+ }
+
+ return true;
+ }
// Find parameters of the inscribed circle that is tangent to three
// point sites.
@@ -1683,7 +1476,6 @@
const site_event<T> &site2,
const site_event<T> &site3,
circle_event<T> &c_event) {
- //return create_circle_event_ppp_gmpxx(site1, site2, site3, c_event);
double dif_x1 = site1.x() - site2.x();
double dif_x2 = site2.x() - site3.x();
double dif_y1 = site1.y() - site2.y();
@@ -1691,28 +1483,36 @@
double orientation = robust_cross_product(dif_x1, dif_y1, dif_x2, dif_y2);
if (orientation_test(orientation) != RIGHT_ORIENTATION)
return false;
- double inv_orientation = 0.5 / orientation;
+ robust_fpt<T> inv_orientation(0.5 / orientation, 3.0);
double sum_x1 = site1.x() + site2.x();
double sum_x2 = site2.x() + site3.x();
double sum_y1 = site1.y() + site2.y();
double sum_y2 = site2.y() + site3.y();
double dif_x3 = site1.x() - site3.x();
double dif_y3 = site1.y() - site3.y();
- epsilon_robust_comparator<T> c_x, c_y;
- c_x += dif_x1 * sum_x1 * dif_y2;
- c_x += dif_y1 * sum_y1 * dif_y2;
- c_x -= dif_x2 * sum_x2 * dif_y1;
- c_x -= dif_y2 * sum_y2 * dif_y1;
- c_y += dif_x2 * sum_x2 * dif_x1;
- c_y += dif_y2 * sum_y2 * dif_x1;
- c_y -= dif_x1 * sum_x1 * dif_x2;
- c_y -= dif_y1 * sum_y1 * dif_x2;
- c_x *= inv_orientation;
- c_y *= inv_orientation;
- epsilon_robust_comparator<T> lower_x(c_x);
- lower_x += sqrt(sqr_distance(dif_x1, dif_y1) * sqr_distance(dif_x2, dif_y2) *
- sqr_distance(dif_x3, dif_y3)) * fabs(inv_orientation);
- c_event = circle_event<double>(c_x, c_y, lower_x);
+ epsilon_robust_comparator< robust_fpt<T> > c_x, c_y;
+ c_x += robust_fpt<T>(dif_x1 * sum_x1 * dif_y2, 2.0);
+ c_x += robust_fpt<T>(dif_y1 * sum_y1 * dif_y2, 2.0);
+ c_x -= robust_fpt<T>(dif_x2 * sum_x2 * dif_y1, 2.0);
+ c_x -= robust_fpt<T>(dif_y2 * sum_y2 * dif_y1, 2.0);
+ c_y += robust_fpt<T>(dif_x2 * sum_x2 * dif_x1, 2.0);
+ c_y += robust_fpt<T>(dif_y2 * sum_y2 * dif_x1, 2.0);
+ c_y -= robust_fpt<T>(dif_x1 * sum_x1 * dif_x2, 2.0);
+ c_y -= robust_fpt<T>(dif_y1 * sum_y1 * dif_x2, 2.0);
+ epsilon_robust_comparator< robust_fpt<T> > lower_x(c_x);
+ lower_x -= robust_fpt<T>(std::sqrt(sqr_distance(dif_x1, dif_y1) *
+ sqr_distance(dif_x2, dif_y2) *
+ sqr_distance(dif_x3, dif_y3)), 5.0);
+ c_event = circle_event<double>(c_x.dif().fpv() * inv_orientation.fpv(),
+ c_y.dif().fpv() * inv_orientation.fpv(),
+ lower_x.dif().fpv() * inv_orientation.fpv());
+ bool recompute_c_x = c_x.dif().ulp() >= 128;
+ bool recompute_c_y = c_y.dif().ulp() >= 128;
+ bool recompute_lower_x = lower_x.dif().ulp() >= 128;
+ if (recompute_c_x || recompute_c_y || recompute_lower_x) {
+ return create_circle_event_ppp_gmpxx(
+ site1, site2, site3, c_event, recompute_c_x, recompute_c_y, recompute_lower_x);
+ }
return true;
}
@@ -1724,7 +1524,7 @@
const site_event<T> &site3,
int segment_index,
circle_event<T> &c_event) {
- //return create_circle_event_pps_gmpxx(site1, site2, site3, segment_index, c_event);
+ return create_circle_event_pps_gmpxx(site1, site2, site3, segment_index, c_event);
if (segment_index != 2) {
kOrientation orient1 = orientation_test(site1.point0(),
site2.point0(), site3.point0(true));
@@ -1757,13 +1557,13 @@
site3.point1().y() - site2.y(),
site2.x() - site3.point1().x());
double denom = robust_cross_product(vec_x, vec_y, line_a, line_b);
- double inv_segm_len = 1.0 / sqrt(sqr_distance(line_a, line_b));
+ double inv_segm_len = 1.0 / std::sqrt(sqr_distance(line_a, line_b));
epsilon_robust_comparator<double> t;
if (orientation_test(denom) == COLLINEAR) {
t += teta / (4.0 * (A + B));
t -= A * B / (teta * (A + B));
} else {
- double det = sqrt((teta * teta + denom * denom) * A * B);
+ double det = std::sqrt((teta * teta + denom * denom) * A * B);
if (segment_index == 2) {
t -= det / (denom * denom);
} else {
@@ -1795,7 +1595,7 @@
const site_event<T> &site3,
int point_index,
circle_event<T> &c_event) {
- //return create_circle_event_pss_gmpxx(site1, site2, site3, point_index, c_event);
+ return create_circle_event_pss_gmpxx(site1, site2, site3, point_index, c_event);
if (site2.index() == site3.index()) {
return false;
}
@@ -1830,9 +1630,9 @@
t -= a1 * ((segm_start1.x() + segm_start2.x()) * 0.5 - site1.x());
t -= b1 * ((segm_start1.y() + segm_start2.y()) * 0.5 - site1.y());
if (point_index == 2) {
- t += sqrt(det);
+ t += std::sqrt(det);
} else {
- t -= sqrt(det);
+ t -= std::sqrt(det);
}
t /= a;
epsilon_robust_comparator<double> c_x, c_y;
@@ -1841,12 +1641,12 @@
c_y += 0.5 * (segm_start1.y() + segm_start2.y());
c_y += b1 * t;
epsilon_robust_comparator<double> lower_x(c_x);
- lower_x += 0.5 * fabs(c) / sqrt(a);
+ lower_x += 0.5 * fabs(c) / std::sqrt(a);
c_event = circle_event<double>(c_x, c_y, lower_x);
return true;
} else {
- double sqr_sum1 = sqrt(a1 * a1 + b1 * b1);
- double sqr_sum2 = sqrt(a2 * a2 + b2 * b2);
+ double sqr_sum1 = std::sqrt(a1 * a1 + b1 * b1);
+ double sqr_sum2 = std::sqrt(a2 * a2 + b2 * b2);
double a = robust_cross_product(a1, b1, -b2, a2);
if (a >= 0) {
a += sqr_sum1 * sqr_sum2;
@@ -1875,9 +1675,9 @@
b -= sqr_sum2 * robust_cross_product(a1, b1, -site1.y(), site1.x());
t -= b;
if (point_index == 2) {
- t += sqrt(det);
+ t += std::sqrt(det);
} else {
- t -= sqrt(det);
+ t -= std::sqrt(det);
}
t /= (a * a);
epsilon_robust_comparator<double> c_x(ix), c_y(iy);
@@ -1888,6 +1688,7 @@
t.abs();
epsilon_robust_comparator<double> lower_x(c_x);
lower_x += t * fabs(orientation);
+
c_event = circle_event<double>(c_x, c_y, lower_x);
}
return true;
@@ -1900,29 +1701,29 @@
const site_event<T> &site2,
const site_event<T> &site3,
circle_event<T> &c_event) {
- //return create_circle_event_sss_gmpxx(site1, site2, site3, c_event);
+ return create_circle_event_sss_gmpxx(site1, site2, site3, c_event);
if (site1.index() == site2.index() ||
site2.index() == site3.index())
return false;
- double a1 = site1.x1(true) - site1.x0(true);
- double b1 = site1.y1(true) - site1.y0(true);
- double c1 = robust_cross_product(site1.x0(true), site1.y0(true),
- site1.x1(true), site1.y1(true));
- double a2 = site2.x1(true) - site2.x0(true);
- double b2 = site2.y1(true) - site2.y0(true);
- double c2 = robust_cross_product(site2.x0(true), site2.y0(true),
- site2.x1(true), site2.y1(true));
- double a3 = site3.x1(true) - site3.x0(true);
- double b3 = site3.y1(true) - site3.y0(true);
- double c3 = robust_cross_product(site3.x0(true), site3.y0(true),
- site3.x1(true), site3.y1(true));
- double len1 = sqrt(a1 * a1 + b1 * b1);
- double len2 = sqrt(a2 * a2 + b2 * b2);
- double len3 = sqrt(a3 * a3 + b3 * b3);
- double cross_12 = robust_cross_product(a1, b1, a2, b2);
- double cross_23 = robust_cross_product(a2, b2, a3, b3);
- double cross_31 = robust_cross_product(a3, b3, a1, b1);
- epsilon_robust_comparator<double> denom, c_x, c_y, r;
+ robust_fpt<T> a1(site1.x1(true) - site1.x0(true), 0.0);
+ robust_fpt<T> b1(site1.y1(true) - site1.y0(true), 0.0);
+ robust_fpt<T> c1(robust_cross_product(site1.x0(true), site1.y0(true), site1.x1(true), site1.y1(true)), 2.0);
+
+ robust_fpt<T> a2(site2.x1(true) - site2.x0(true), 0.0);
+ robust_fpt<T> b2(site2.y1(true) - site2.y0(true), 0.0);
+ robust_fpt<T> c2(robust_cross_product(site2.x0(true), site2.y0(true), site2.x1(true), site2.y1(true)), 2.0);
+
+ robust_fpt<T> a3(site3.x1(true) - site3.x0(true), 0.0);
+ robust_fpt<T> b3(site3.y1(true) - site3.y0(true), 0.0);
+ robust_fpt<T> c3(robust_cross_product(site3.x0(true), site3.y0(true), site3.x1(true), site3.y1(true)), 2.0);
+
+ robust_fpt<T> len1 = (a1 * a1 + b1 * b1).get_sqrt();
+ robust_fpt<T> len2 = (a2 * a2 + b2 * b2).get_sqrt();
+ robust_fpt<T> len3 = (a3 * a3 + b3 * b3).get_sqrt();
+ robust_fpt<T> cross_12(robust_cross_product(a1.fpv(), b1.fpv(), a2.fpv(), b2.fpv()), 2.0);
+ robust_fpt<T> cross_23(robust_cross_product(a2.fpv(), b2.fpv(), a3.fpv(), b3.fpv()), 2.0);
+ robust_fpt<T> cross_31(robust_cross_product(a3.fpv(), b3.fpv(), a1.fpv(), b1.fpv()), 2.0);
+ epsilon_robust_comparator< robust_fpt<T> > denom, c_x, c_y, r;
// denom = cross_12 * len3 + cross_23 * len1 + cross_31 * len2.
denom += cross_12 * len3;
@@ -1946,7 +1747,20 @@
c_y -= b3 * c2 * len1;
c_y += b3 * c1 * len2;
c_y -= b1 * c3 * len2;
- c_event = circle_event<double>(c_x, c_y, c_x + r, denom);
+ epsilon_robust_comparator< robust_fpt<T> > lower_x(c_x + r);
+ bool recompute_c_x = c_x.dif().re() > 128;
+ bool recompute_c_y = c_y.dif().re() > 128;
+ bool recompute_lower_x = lower_x.dif() > 128;
+ bool recompute_denom = denom.dif().re() > 128;
+ c_event = circle_event<double>(c_x.dif().fpv() / denom.dif().fpv(),
+ c_y.dif().fpv() / denom.dif().fpv(),
+ lower_x.dif().fpv() / denom.dif().fpv());
+ if (recompute_c_x || recompute_c_y || recompute_lower_x || recompute_denom) {
+ return create_circle_event_sss_gmpxx(
+ site1, site2, site3, c_event,
+ recompute_c_x, recompute_c_y,
+ recompute_lower_x, recompute_denom);
+ }
return true;
}
Modified: sandbox/SOC/2010/sweepline/boost/sweepline/voronoi_sweepline.hpp
==============================================================================
--- sandbox/SOC/2010/sweepline/boost/sweepline/voronoi_sweepline.hpp (original)
+++ sandbox/SOC/2010/sweepline/boost/sweepline/voronoi_sweepline.hpp 2011-05-02 17:59:42 EDT (Mon, 02 May 2011)
@@ -28,11 +28,13 @@
#endif
#endif
-//#pragma warning(disable:4800)
-//#include <gmpxx.h>
+#pragma warning(disable:4800)
+#include <gmpxx.h>
#include "voronoi_output.hpp"
-#include "detail/mpz_arithmetic.hpp"
+#include "detail/mpt_wrapper.hpp"
+#include "detail/robust_fpt.hpp"
+#include "detail/sqrt_expr_evaluator.hpp"
#include "detail/voronoi_formation.hpp"
namespace boost {
Modified: sandbox/SOC/2010/sweepline/libs/sweepline/test/Jamfile.v2
==============================================================================
--- sandbox/SOC/2010/sweepline/libs/sweepline/test/Jamfile.v2 (original)
+++ sandbox/SOC/2010/sweepline/libs/sweepline/test/Jamfile.v2 2011-05-02 17:59:42 EDT (Mon, 02 May 2011)
@@ -11,7 +11,7 @@
;
-alias "sweepline"
+alias "main"
:
[ run circle_event_creation_test.cpp ]
[ run clipping_test.cpp ]
@@ -19,41 +19,12 @@
[ run event_types_test.cpp ]
[ run node_comparer_test.cpp ]
[ run predicates_test.cpp ]
+ [ run robust_fpt_test.cpp ]
+ [ run sqrt_expr_evaluator_test.cpp ]
[ run sweepline_test.cpp ]
;
-alias "sweepline_benchmark"
+alias "benchmark"
:
[ run benchmark_test.cpp ]
;
-
-path-constant GMP_PATH : "/usr/local" ;
-
-lib gmp
- : # sources
- : # requirements
- <name>gmp <search>$(GMP_PATH)/lib
- : #default-build
- : # usage-requirements
- ;
-
-lib gmpxx
- : # sources
- : # requirements
- <name>gmpxx <search>$(GMP_PATH)/lib
- : # default-build
- : # usage-requirements
- ;
-
-alias gmp_libs : gmp gmpxx : <link>static ;
-
-alias "mpz_arithmetic"
- :
- [ run mpz_arithmetic_test.cpp gmp_libs
- : # args
- : # input files
- : # requirements
- <include>$(GMP_PATH)/include
- : # target-name
- ]
- ;
Modified: sandbox/SOC/2010/sweepline/libs/sweepline/test/benchmark_test.cpp
==============================================================================
--- sandbox/SOC/2010/sweepline/libs/sweepline/test/benchmark_test.cpp (original)
+++ sandbox/SOC/2010/sweepline/libs/sweepline/test/benchmark_test.cpp 2011-05-02 17:59:42 EDT (Mon, 02 May 2011)
@@ -9,55 +9,55 @@
#include <iostream>
#include <stdio.h>
-#include <stdlib.h>
-#include <time.h>
+
+#define BOOST_TEST_MODULE benchmark_test
+#include <boost/random/mersenne_twister.hpp>
+#include <boost/test/test_case_template.hpp>
+#include <boost/timer.hpp>
#include "test_type_list.hpp"
#include "boost/sweepline/voronoi_sweepline.hpp"
using namespace boost::sweepline;
-#define BOOST_TEST_MODULE voronoi_benchmark_test
-#include <boost/test/test_case_template.hpp>
-
#ifdef WIN32
#pragma warning( disable : 4996 )
#endif
BOOST_AUTO_TEST_CASE_TEMPLATE(benchmark_test1, T, test_types) {
typedef T coordinate_type;
- srand(static_cast<unsigned int>(time(NULL)));
+ boost::mt19937 gen(static_cast<unsigned int>(time(NULL)));
+ boost::timer timer;
voronoi_output<coordinate_type> test_output;
FILE *bench_file = fopen("benchmark.txt", "a");
fprintf(bench_file, "Voronoi Segment Sweepline Benchmark Test (time in seconds):\n");
+ fprintf(bench_file, "| Number of points | Number of tests | Time per one test |\n");
#ifdef NDEBUG
int max_points = 1000000;
#else
- int max_points = 1000;
+ int max_points = 100;
#endif
for (int num_points = 10; num_points <= max_points; num_points *= 10) {
std::vector< point_2d<T> > points;
- points.reserve(num_points);
+ points.resize(num_points);
- clock_t start_time = clock();
- int num_times = max_points / num_points;
- for (int cur = 0; cur < num_times; cur++) {
+ timer.restart();
+ int num_tests = max_points / num_points;
+ for (int cur = 0; cur < num_tests; cur++) {
for (int cur_point = 0; cur_point < num_points; cur_point++) {
- points.push_back(make_point_2d<coordinate_type>(
- static_cast<coordinate_type>(rand() % 5000 - 10000),
- static_cast<coordinate_type>(rand() % 5000 - 10000)));
+ points[cur_point] = point_2d<coordinate_type>(
+ static_cast<int>(gen()), static_cast<int>(gen()));
}
build_voronoi(points, test_output);
- points.clear();
}
- clock_t end_time = clock();
- double running_time = static_cast<double>(end_time - start_time) / CLOCKS_PER_SEC / num_times;
+ double elapsed_time = timer.elapsed();
+ double time_per_test = elapsed_time / num_tests;
- fprintf(bench_file,
- "Number of points = %8d; Overall time = %2d; Time per one input = %9.6f.\n",
- num_points, static_cast<int>(end_time - start_time), running_time);
+ fprintf(bench_file, "| %16d ", num_points);
+ fprintf(bench_file, "| %15d ", num_tests);
+ fprintf(bench_file, "| %17.6f |\n", time_per_test);
}
fclose(bench_file);
}
Deleted: sandbox/SOC/2010/sweepline/libs/sweepline/test/mpz_arithmetic_test.cpp
==============================================================================
--- sandbox/SOC/2010/sweepline/libs/sweepline/test/mpz_arithmetic_test.cpp 2011-05-02 17:59:42 EDT (Mon, 02 May 2011)
+++ (empty file)
@@ -1,71 +0,0 @@
-// Boost sweepline library mpz_arithmetic_test.cpp file
-
-// Copyright Andrii Sydorchuk 2010.
-// 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.
-
-#include <stdlib.h>
-#include <time.h>
-
-#pragma warning(disable:4800)
-#include <gmpxx.h>
-#include "test_type_list.hpp"
-
-#include "boost/sweepline/detail/mpz_arithmetic.hpp"
-using namespace boost::sweepline::detail;
-
-#define BOOST_TEST_MODULE mpz_arithmetic_test
-#include <boost/test/test_case_template.hpp>
-
-template <int N>
-struct test_sqr_evaluator {
- template <typename mpt>
- static bool run() {
- bool ret_val = true;
- static mpt A[N], B[N];
- double a[N], b[N];
- for (int i = 0; i < N; ++i) {
- a[i] = rand() % 1000 - 500;
- int rand_number = rand() % 1000 - 500;
- b[i] = rand_number * rand_number;
- }
- int mask = (1 << N);
- for (int i = 0; i < mask; i++) {
- double expected_val = 0.0;
- for (int j = 0; j < N; j++) {
- if (i & (1 << j)) {
- A[j] = a[j];
- B[j] = b[j];
- expected_val += a[j] * sqrt(b[j]);
- } else {
- A[j] = -a[j];
- B[j] = b[j];
- expected_val -= a[j] * sqrt(b[j]);
- }
- }
- double received_val = sqr_expr_evaluator<N>::eval(A, B);
- // TODO(asydorchuk): Change to almost equal check.
- ret_val &= (fabs(expected_val - received_val) <= 1E-9);
- }
- return ret_val;
- }
-};
-
-BOOST_AUTO_TEST_CASE(mpz_sqr_evaluator_test) {
- srand(static_cast<unsigned int>(time(0)));
- typedef mpt_wrapper<mpz_class, 4> mpt_wrapper_type;
- for (int i = 0; i < 1000; i++) {
- BOOST_CHECK(test_sqr_evaluator<2>::run<mpz_class>());
- BOOST_CHECK(test_sqr_evaluator<3>::run<mpz_class>());
- BOOST_CHECK(test_sqr_evaluator<4>::run<mpz_class>());
- }
- for (int i = 0; i < 1000; i++) {
- BOOST_CHECK(test_sqr_evaluator<2>::run<mpt_wrapper_type>());
- BOOST_CHECK(test_sqr_evaluator<3>::run<mpt_wrapper_type>());
- BOOST_CHECK(test_sqr_evaluator<4>::run<mpt_wrapper_type>());
- }
-}
-
Added: sandbox/SOC/2010/sweepline/libs/sweepline/test/robust_fpt_test.cpp
==============================================================================
--- (empty file)
+++ sandbox/SOC/2010/sweepline/libs/sweepline/test/robust_fpt_test.cpp 2011-05-02 17:59:42 EDT (Mon, 02 May 2011)
@@ -0,0 +1,144 @@
+// Boost sweepline library robust_fpt_test.cpp file
+
+// Copyright Andrii Sydorchuk 2010.
+// 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.
+
+#define BOOST_TEST_MODULE robust_fpt_test
+#include <boost/mpl/list.hpp>
+#include <boost/test/test_case_template.hpp>
+
+#include "boost/sweepline/voronoi_sweepline.hpp"
+using namespace boost::sweepline::detail;
+
+typedef boost::mpl::list<double, mpf_class, mpt_wrapper<mpf_class, 8> > test_types;
+
+BOOST_AUTO_TEST_CASE_TEMPLATE(empty_constructor_test, T, test_types) {
+ typedef robust_fpt<T> rfpt_type;
+ rfpt_type a = rfpt_type();
+ BOOST_CHECK_EQUAL(a.fpv() == static_cast<T>(0), true);
+ BOOST_CHECK_EQUAL(a.re() == 0.0, true);
+ BOOST_CHECK_EQUAL(a.ulp() == 0, true);
+}
+
+BOOST_AUTO_TEST_CASE_TEMPLATE(autorounded_constructor_test, T, test_types) {
+ typedef robust_fpt<T> rfpt_type;
+ rfpt_type a(static_cast<T>(10));
+ BOOST_CHECK_EQUAL(a.fpv() == static_cast<T>(10), true);
+ BOOST_CHECK_EQUAL(a.re() == 1.0, true);
+ BOOST_CHECK_EQUAL(a.ulp() == 1.0, true);
+}
+
+BOOST_AUTO_TEST_CASE_TEMPLATE(notrounded_constructor_test, T, test_types) {
+ typedef robust_fpt<T> rfpt_type;
+ rfpt_type a(static_cast<T>(10), false);
+ BOOST_CHECK_EQUAL(a.fpv() == static_cast<T>(10), true);
+ BOOST_CHECK_EQUAL(a.re() == 0.0, true);
+ BOOST_CHECK_EQUAL(a.ulp() == 0.0, true);
+}
+
+BOOST_AUTO_TEST_CASE_TEMPLATE(constructor_test, T, test_types) {
+ typedef robust_fpt<T> rfpt_type;
+ rfpt_type a(static_cast<T>(10), 3.0);
+ BOOST_CHECK_EQUAL(a.fpv() == static_cast<T>(10), true);
+ BOOST_CHECK_EQUAL(a.re() == 3.0, true);
+ BOOST_CHECK_EQUAL(a.ulp() == 3.0, true);
+
+ rfpt_type b(static_cast<T>(10), 2.75);
+ BOOST_CHECK_EQUAL(b.fpv() == static_cast<T>(10), true);
+ BOOST_CHECK_EQUAL(b.re() == 2.75, true);
+ BOOST_CHECK_EQUAL(b.ulp() == 2.75, true);
+}
+
+BOOST_AUTO_TEST_CASE_TEMPLATE(sum_test1, T, test_types) {
+ typedef robust_fpt<T> rfpt_type;
+ rfpt_type a(static_cast<T>(2), 5.0);
+ rfpt_type b(static_cast<T>(3), 4.0);
+ rfpt_type c = a + b;
+ BOOST_CHECK_EQUAL(c.fpv() == static_cast<T>(5), true);
+ BOOST_CHECK_EQUAL(c.re() == 6.0, true);
+ BOOST_CHECK_EQUAL(c.ulp() == 6.0, true);
+
+ c += b;
+ BOOST_CHECK_EQUAL(c.fpv() == static_cast<T>(8), true);
+ BOOST_CHECK_EQUAL(c.re() == 7.0, true);
+ BOOST_CHECK_EQUAL(c.ulp() == 7.0, true);
+}
+
+BOOST_AUTO_TEST_CASE_TEMPLATE(sum_test2, T, test_types) {
+ typedef robust_fpt<T> rfpt_type;
+ rfpt_type a(static_cast<T>(3), 2.0);
+ rfpt_type b(static_cast<T>(-2), 3.0);
+ rfpt_type c = a + b;
+ BOOST_CHECK_EQUAL(c.fpv() == static_cast<T>(1), true);
+ BOOST_CHECK_EQUAL(c.re() == 13.0, true);
+ BOOST_CHECK_EQUAL(c.ulp() == 13.0, true);
+
+ c += b;
+ BOOST_CHECK_EQUAL(c.fpv() == static_cast<T>(-1), true);
+ BOOST_CHECK_EQUAL(c.re() == 20.0, true);
+ BOOST_CHECK_EQUAL(c.ulp() == 20.0, true);
+}
+
+BOOST_AUTO_TEST_CASE_TEMPLATE(minus_test1, T, test_types) {
+ typedef robust_fpt<T> rfpt_type;
+ rfpt_type a(static_cast<T>(2), 5.0);
+ rfpt_type b(static_cast<T>(-3), 4.0);
+ rfpt_type c = a - b;
+ BOOST_CHECK_EQUAL(c.fpv() == static_cast<T>(5), true);
+ BOOST_CHECK_EQUAL(c.re() == 6.0, true);
+ BOOST_CHECK_EQUAL(c.ulp() == 6.0, true);
+
+ c -= b;
+ BOOST_CHECK_EQUAL(c.fpv() == static_cast<T>(8), true);
+ BOOST_CHECK_EQUAL(c.re() == 7.0, true);
+ BOOST_CHECK_EQUAL(c.ulp() == 7.0, true);
+}
+
+BOOST_AUTO_TEST_CASE_TEMPLATE(minus_test2, T, test_types) {
+ typedef robust_fpt<T> rfpt_type;
+ rfpt_type a(static_cast<T>(3), 2.0);
+ rfpt_type b(static_cast<T>(2), 3.0);
+ rfpt_type c = a - b;
+ BOOST_CHECK_EQUAL(c.fpv() == static_cast<T>(1), true);
+ BOOST_CHECK_EQUAL(c.re() == 13.0, true);
+ BOOST_CHECK_EQUAL(c.ulp() == 13.0, true);
+
+ c -= b;
+ BOOST_CHECK_EQUAL(c.fpv() == static_cast<T>(-1), true);
+ BOOST_CHECK_EQUAL(c.re() == 20.0, true);
+ BOOST_CHECK_EQUAL(c.ulp() == 20.0, true);
+}
+
+BOOST_AUTO_TEST_CASE_TEMPLATE(mult_test, T, test_types) {
+ typedef robust_fpt<T> rfpt_type;
+ rfpt_type a(static_cast<T>(2), 3.0);
+ rfpt_type b(static_cast<T>(4));
+ rfpt_type c = a * b;
+ BOOST_CHECK_EQUAL(c.fpv() == static_cast<T>(8), true);
+ BOOST_CHECK_EQUAL(c.re() == 5.0, true);
+ BOOST_CHECK_EQUAL(c.ulp() == 5.0, true);
+
+ c *= b;
+ BOOST_CHECK_EQUAL(c.fpv() == static_cast<T>(32), true);
+ BOOST_CHECK_EQUAL(c.re() == 7.0, true);
+ BOOST_CHECK_EQUAL(c.ulp() == 7.0, true);
+}
+
+BOOST_AUTO_TEST_CASE_TEMPLATE(division_test, T, test_types) {
+ typedef robust_fpt<T> rfpt_type;
+ rfpt_type a(static_cast<T>(2), 3.0);
+ rfpt_type b(static_cast<T>(4));
+ rfpt_type c = a / b;
+ BOOST_CHECK_EQUAL(c.fpv() == static_cast<T>(0.5), true);
+ BOOST_CHECK_EQUAL(c.re() == 5.0, true);
+ BOOST_CHECK_EQUAL(c.ulp() == 5.0, true);
+
+ c /= b;
+ BOOST_CHECK_EQUAL(c.fpv() == static_cast<T>(0.125), true);
+ BOOST_CHECK_EQUAL(c.re() == 7.0, true);
+ BOOST_CHECK_EQUAL(c.ulp() == 7.0, true);
+}
Copied: sandbox/SOC/2010/sweepline/libs/sweepline/test/sqrt_expr_evaluator_test.cpp (from r70966, /sandbox/SOC/2010/sweepline/libs/sweepline/test/mpz_arithmetic_test.cpp)
==============================================================================
--- /sandbox/SOC/2010/sweepline/libs/sweepline/test/mpz_arithmetic_test.cpp (original)
+++ sandbox/SOC/2010/sweepline/libs/sweepline/test/sqrt_expr_evaluator_test.cpp 2011-05-02 17:59:42 EDT (Mon, 02 May 2011)
@@ -1,4 +1,4 @@
-// Boost sweepline library mpz_arithmetic_test.cpp file
+// Boost sweepline library sqrt_expr_evaluator_test.cpp file
// Copyright Andrii Sydorchuk 2010.
// Distributed under the Boost Software License, Version 1.0.
@@ -10,22 +10,18 @@
#include <stdlib.h>
#include <time.h>
-#pragma warning(disable:4800)
-#include <gmpxx.h>
-#include "test_type_list.hpp"
+#define BOOST_TEST_MODULE sqrt_expr_evaluator_test
+#include <boost/test/test_case_template.hpp>
-#include "boost/sweepline/detail/mpz_arithmetic.hpp"
+#include "boost/sweepline/voronoi_sweepline.hpp"
using namespace boost::sweepline::detail;
-#define BOOST_TEST_MODULE mpz_arithmetic_test
-#include <boost/test/test_case_template.hpp>
-
template <int N>
struct test_sqr_evaluator {
- template <typename mpt>
+ template <typename mpz, typename mpf>
static bool run() {
bool ret_val = true;
- static mpt A[N], B[N];
+ static mpz A[N], B[N];
double a[N], b[N];
for (int i = 0; i < N; ++i) {
a[i] = rand() % 1000 - 500;
@@ -46,9 +42,8 @@
expected_val -= a[j] * sqrt(b[j]);
}
}
- double received_val = sqr_expr_evaluator<N>::eval(A, B);
- // TODO(asydorchuk): Change to almost equal check.
- ret_val &= (fabs(expected_val - received_val) <= 1E-9);
+ mpf received_val = (sqr_expr_evaluator<N>::template eval<mpz, mpf>(A, B));
+ ret_val &= almost_equal(expected_val, received_val.get_d(), 1);
}
return ret_val;
}
@@ -56,16 +51,19 @@
BOOST_AUTO_TEST_CASE(mpz_sqr_evaluator_test) {
srand(static_cast<unsigned int>(time(0)));
- typedef mpt_wrapper<mpz_class, 4> mpt_wrapper_type;
+ typedef mpt_wrapper<mpz_class, 4> mpz_wrapper_type;
+ typedef mpt_wrapper<mpf_class, 1> mpf_wrapper_type;
for (int i = 0; i < 1000; i++) {
- BOOST_CHECK(test_sqr_evaluator<2>::run<mpz_class>());
- BOOST_CHECK(test_sqr_evaluator<3>::run<mpz_class>());
- BOOST_CHECK(test_sqr_evaluator<4>::run<mpz_class>());
+ BOOST_CHECK((test_sqr_evaluator<1>::run<mpz_class, mpf_class>()));
+ BOOST_CHECK((test_sqr_evaluator<2>::run<mpz_class, mpf_class>()));
+ BOOST_CHECK((test_sqr_evaluator<3>::run<mpz_class, mpf_class>()));
+ BOOST_CHECK((test_sqr_evaluator<4>::run<mpz_class, mpf_class>()));
}
for (int i = 0; i < 1000; i++) {
- BOOST_CHECK(test_sqr_evaluator<2>::run<mpt_wrapper_type>());
- BOOST_CHECK(test_sqr_evaluator<3>::run<mpt_wrapper_type>());
- BOOST_CHECK(test_sqr_evaluator<4>::run<mpt_wrapper_type>());
- }
+ BOOST_CHECK((test_sqr_evaluator<1>::run<mpz_wrapper_type, mpf_wrapper_type>()));
+ BOOST_CHECK((test_sqr_evaluator<2>::run<mpz_wrapper_type, mpf_wrapper_type>()));
+ BOOST_CHECK((test_sqr_evaluator<3>::run<mpz_wrapper_type, mpf_wrapper_type>()));
+ BOOST_CHECK((test_sqr_evaluator<4>::run<mpz_wrapper_type, mpf_wrapper_type>()));
+ }
}
Modified: sandbox/SOC/2010/sweepline/libs/sweepline/test/sweepline_test.cpp
==============================================================================
--- sandbox/SOC/2010/sweepline/libs/sweepline/test/sweepline_test.cpp (original)
+++ sandbox/SOC/2010/sweepline/libs/sweepline/test/sweepline_test.cpp 2011-05-02 17:59:42 EDT (Mon, 02 May 2011)
@@ -1,4 +1,4 @@
-// Boost sweepline library builder_test.cpp file
+// Boost sweepline library sweepline_test.cpp file
// Copyright Andrii Sydorchuk 2010.
// Distributed under the Boost Software License, Version 1.0.
@@ -15,7 +15,7 @@
using namespace boost::sweepline;
#include "output_verification.hpp"
-#define BOOST_TEST_MODULE voronoi_builder_test
+#define BOOST_TEST_MODULE sweepline_test
#include <boost/test/test_case_template.hpp>
#define CHECK_EQUAL_POINTS(p1, p2) \
@@ -24,10 +24,6 @@
#define VERIFY_VORONOI_OUTPUT(output, mask) BOOST_CHECK_EQUAL(verify_output(output, mask), true)
-#define ALMOST_EQUAL_TEST(a, b, ULP_ERR, ABS_ERR) \
- BOOST_CHECK_EQUAL(detail::almost_equal(a, b, ULP_ERR) || \
- detail::abs_equal(a, b, ABS_ERR), true);
-
// Sites: (0, 0).
BOOST_AUTO_TEST_CASE_TEMPLATE(single_site_test, T, test_types) {
typedef T coordinate_type;
@@ -670,52 +666,52 @@
}
#endif
-//#ifdef NDEBUG
-//BOOST_AUTO_TEST_CASE_TEMPLATE(segment_random_test2, T, test_types) {
-// srand(static_cast<unsigned int>(time(NULL)));
-// voronoi_output<T> test_output_small, test_output_large;
-// std::vector< std::pair< point_2d<T>, point_2d<T> > > segm_vec;
-// int num_segments[] = {10, 100, 1000, 10000};
-// int num_runs[] = {1000, 100, 10, 1};
-// int mod_koef1[] = {100, 1000, 10000, 100000};
-// int mod_koef2[] = {100, 200, 300, 400};
-// int max_value[] = {100, 600, 5150, 50200};
-// int array_length = sizeof(num_segments) / sizeof(int);
-// for (int k = 0; k < 4; k++) {
-// int koef = std::numeric_limits<int>::max() / max_value[k];
-// for (int i = 0; i < num_runs[k]; i++) {
-// for (int j = 0; j < num_segments[k]; j++) {
-// T x1 = (rand() % (mod_koef1[k] / 100)) - mod_koef1[k] / 2;
-// T y1 = (rand() % (mod_koef1[k] / 100)) - mod_koef1[k] / 2;
-// T dx = 0, dy = 0;
-// while (dx == 0 && dy == 0) {
-// dx = (rand() % mod_koef2[k]) - mod_koef2[k] / 2;
-// dy = (rand() % mod_koef2[k]) - mod_koef2[k] / 2;
-// }
-// T x2 = x1 + dx;
-// T y2 = y1 + dy;
-// point_2d<T> point1_small(x1, y1);
-// point_2d<T> point2_small(x2, y2);
-// segm_vec.push_back(std::make_pair(point1_small, point2_small));
-// }
-// remove_intersections(segm_vec);
-// build_voronoi(segm_vec, test_output_small);
-// for (size_t j = 0; j < segm_vec.size(); j++) {
-// segm_vec[j].first.x(segm_vec[j].first.x() * koef);
-// segm_vec[j].first.y(segm_vec[j].first.y() * koef);
-// segm_vec[j].second.x(segm_vec[j].second.x() * koef);
-// segm_vec[j].second.y(segm_vec[j].second.y() * koef);
-// }
-// build_voronoi(segm_vec, test_output_large);
-// VERIFY_VORONOI_OUTPUT(test_output_large, NO_HALF_EDGE_INTERSECTIONS);
-// BOOST_CHECK_EQUAL(test_output_small.num_cell_records(),
-// test_output_large.num_cell_records());
-// BOOST_CHECK_EQUAL(test_output_small.num_vertex_records(),
-// test_output_large.num_vertex_records());
-// BOOST_CHECK_EQUAL(test_output_small.num_edge_records(),
-// test_output_large.num_edge_records());
-// segm_vec.clear();
-// }
-// }
-//}
-//#endif
+#ifdef NDEBUG
+BOOST_AUTO_TEST_CASE_TEMPLATE(segment_random_test2, T, test_types) {
+ srand(static_cast<unsigned int>(time(NULL)));
+ voronoi_output<T> test_output_small, test_output_large;
+ std::vector< std::pair< point_2d<T>, point_2d<T> > > segm_vec;
+ int num_segments[] = {10, 100, 1000, 10000};
+ int num_runs[] = {1000, 100, 10, 1};
+ int mod_koef1[] = {100, 1000, 10000, 100000};
+ int mod_koef2[] = {100, 200, 300, 400};
+ int max_value[] = {100, 600, 5150, 50200};
+ int array_length = sizeof(num_segments) / sizeof(int);
+ for (int k = 0; k < array_length; k++) {
+ int koef = std::numeric_limits<int>::max() / max_value[k];
+ for (int i = 0; i < num_runs[k]; i++) {
+ for (int j = 0; j < num_segments[k]; j++) {
+ T x1 = (rand() % mod_koef1[k]) - mod_koef1[k] / 2;
+ T y1 = (rand() % mod_koef1[k]) - mod_koef1[k] / 2;
+ T dx = 0, dy = 0;
+ while (dx == 0 && dy == 0) {
+ dx = (rand() % mod_koef2[k]) - mod_koef2[k] / 2;
+ dy = (rand() % mod_koef2[k]) - mod_koef2[k] / 2;
+ }
+ T x2 = x1 + dx;
+ T y2 = y1 + dy;
+ point_2d<T> point1_small(x1, y1);
+ point_2d<T> point2_small(x2, y2);
+ segm_vec.push_back(std::make_pair(point1_small, point2_small));
+ }
+ remove_intersections(segm_vec);
+ build_voronoi(segm_vec, test_output_small);
+ for (size_t j = 0; j < segm_vec.size(); j++) {
+ segm_vec[j].first.x(segm_vec[j].first.x() * koef);
+ segm_vec[j].first.y(segm_vec[j].first.y() * koef);
+ segm_vec[j].second.x(segm_vec[j].second.x() * koef);
+ segm_vec[j].second.y(segm_vec[j].second.y() * koef);
+ }
+ build_voronoi(segm_vec, test_output_large);
+ VERIFY_VORONOI_OUTPUT(test_output_large, NO_HALF_EDGE_INTERSECTIONS);
+ BOOST_CHECK_EQUAL(test_output_small.num_cell_records(),
+ test_output_large.num_cell_records());
+ BOOST_CHECK_EQUAL(test_output_small.num_vertex_records(),
+ test_output_large.num_vertex_records());
+ BOOST_CHECK_EQUAL(test_output_small.num_edge_records(),
+ test_output_large.num_edge_records());
+ segm_vec.clear();
+ }
+ }
+}
+#endif