Subject: [Boost-commit] svn:boost r76656 - in sandbox/gtl: boost/polygon boost/polygon/detail libs/polygon/test libs/polygon/voronoi_example
From: sydorchuk.andriy_at_[hidden]
Date: 2012-01-23 17:22:46
Author: asydorchuk
Date: 2012-01-23 17:22:43 EST (Mon, 23 Jan 2012)
New Revision: 76656
URL: http://svn.boost.org/trac/boost/changeset/76656
Log:
Added voronoi_ctypes (coordinate types).
Added voronoi_ctype_traits.
Removed output type specification from template arguments of the voronoi_builder class.
Added:
sandbox/gtl/boost/polygon/detail/voronoi_ctypes.hpp (contents, props changed)
Text files modified:
sandbox/gtl/boost/polygon/detail/voronoi_predicates.hpp | 129 +-
sandbox/gtl/boost/polygon/detail/voronoi_robust_fpt.hpp | 1510 ++++++++++-----------------------------
sandbox/gtl/boost/polygon/detail/voronoi_structures.hpp | 6
sandbox/gtl/boost/polygon/voronoi.hpp | 21
sandbox/gtl/boost/polygon/voronoi_builder.hpp | 87 +-
sandbox/gtl/boost/polygon/voronoi_diagram.hpp | 21
sandbox/gtl/libs/polygon/test/voronoi_benchmark_test.cpp | 9
sandbox/gtl/libs/polygon/voronoi_example/voronoi_visualizer.cpp | 5
8 files changed, 542 insertions(+), 1246 deletions(-)
Added: sandbox/gtl/boost/polygon/detail/voronoi_ctypes.hpp
==============================================================================
--- (empty file)
+++ sandbox/gtl/boost/polygon/detail/voronoi_ctypes.hpp 2012-01-23 17:22:43 EST (Mon, 23 Jan 2012)
@@ -0,0 +1,744 @@
+// Boost.Polygon library detail/voronoi_ctypes.hpp header file
+
+// Copyright Andrii Sydorchuk 2010-2012.
+// Distributed under the Boost Software License, Version 1.0.
+// (See accompanying file LICENSE_1_0.txt or copy at
+// http://www.boost.org/LICENSE_1_0.txt)
+
+// See http://www.boost.org for updates, documentation, and revision history.
+
+#ifndef BOOST_POLYGON_VORONOI_DETAIL_CTYPES
+#define BOOST_POLYGON_VORONOI_DETAIL_CTYPES
+
+#include <cmath>
+#include <boost/cstdint.hpp>
+
+namespace boost {
+namespace polygon {
+namespace detail {
+
+typedef boost::int32_t int32;
+typedef boost::int64_t int64;
+typedef boost::uint32_t uint32;
+typedef boost::uint64_t uint64;
+typedef double fpt64;
+
+// 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 reinterpretations differ in not more than maxUlps units.
+template <typename _fpt>
+struct ulp_comparison;
+
+template <>
+struct ulp_comparison<fpt64> {
+ enum kResult {
+ LESS = -1,
+ EQUAL = 0,
+ MORE = 1
+ };
+
+ kResult operator()(fpt64 a, fpt64 b, unsigned int maxUlps) const {
+ uint64 ll_a, ll_b;
+
+ // Reinterpret double bits as 64-bit signed integer.
+ memcpy(&ll_a, &a, sizeof(fpt64));
+ memcpy(&ll_b, &b, sizeof(fpt64));
+
+ // 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) ? EQUAL : LESS;
+ return (ll_b - ll_a <= maxUlps) ? EQUAL : MORE;
+ }
+};
+
+// Manages exponent of the floating-point value.
+template <typename _fpt>
+struct fpt_exponent_accessor;
+
+template <>
+class fpt_exponent_accessor<fpt64> {
+public:
+ static const int64 kExponentMask;
+ static const int64 kSignedMantissaMask;
+ static const int64 kMinExponent;
+ static const int64 kMaxExponent;
+ static const int64 kMaxSignificantExpDif;
+
+ static int64 set_exponent(fpt64& value, int64 exponent) {
+ int64 bits;
+ memcpy(&bits, &value, sizeof(fpt64));
+ int64 exp = ((bits & kExponentMask) >> 52) - 1023;
+ if (exp == exponent) {
+ return exp;
+ }
+ bits = (bits & kSignedMantissaMask) | ((exponent + 1023) << 52);
+ memcpy(&value, &bits, sizeof(fpt64));
+ return exp;
+ }
+};
+
+const int64 fpt_exponent_accessor<fpt64>::kExponentMask = 0x7ff0000000000000LL;
+const int64 fpt_exponent_accessor<fpt64>::kSignedMantissaMask = 0x800fffffffffffffLL;
+const int64 fpt_exponent_accessor<fpt64>::kMinExponent = -1023LL;
+const int64 fpt_exponent_accessor<fpt64>::kMaxExponent = 1024LL;
+const int64 fpt_exponent_accessor<fpt64>::kMaxSignificantExpDif = 54;
+
+// Allows to extend floating-point type exponent boundaries to the 64 bit
+// integer range. This class does not handle division by zero, subnormal
+// numbers or NaNs.
+template <typename _fpt>
+class extended_exponent_fpt {
+public:
+ typedef _fpt fpt_type;
+ typedef int64 exp_type;
+ typedef fpt_exponent_accessor<fpt_type> fea;
+
+ explicit extended_exponent_fpt(fpt_type value) {
+ if (value == 0.0) {
+ exponent_ = 0;
+ value_ = 0.0;
+ } else {
+ exponent_ = fea::set_exponent(value, 0);
+ value_ = value;
+ }
+ }
+
+ extended_exponent_fpt(fpt_type value, exp_type exponent) {
+ if (value == 0.0) {
+ exponent_ = 0;
+ value_ = 0.0;
+ } else {
+ exponent_ = fea::set_exponent(value, 0) + exponent;
+ value_ = value;
+ }
+ }
+
+ bool is_pos() const {
+ return value_ > 0;
+ }
+
+ bool is_neg() const {
+ return value_ < 0;
+ }
+
+ bool is_zero() const {
+ return value_ == 0;
+ }
+
+ extended_exponent_fpt operator-() const {
+ return extended_exponent_fpt(-value_, exponent_);
+ }
+
+ extended_exponent_fpt operator+(const extended_exponent_fpt& that) const {
+ if (this->value_ == 0.0 ||
+ that.exponent_ > this->exponent_ + fea::kMaxSignificantExpDif) {
+ return that;
+ }
+ if (that.value_ == 0.0 ||
+ this->exponent_ > that.exponent_ + fea::kMaxSignificantExpDif) {
+ return *this;
+ }
+ if (this->exponent_ >= that.exponent_) {
+ exp_type exp_dif = this->exponent_ - that.exponent_;
+ fpt_type value = this->value_;
+ fea::set_exponent(value, exp_dif);
+ value += that.value_;
+ return extended_exponent_fpt(value, that.exponent_);
+ } else {
+ exp_type exp_dif = that.exponent_ - this->exponent_;
+ fpt_type value = that.value_;
+ fea::set_exponent(value, exp_dif);
+ value += this->value_;
+ return extended_exponent_fpt(value, this->exponent_);
+ }
+ }
+
+ extended_exponent_fpt operator-(const extended_exponent_fpt& that) const {
+ if (this->value_ == 0.0 ||
+ that.exponent_ > this->exponent_ + fea::kMaxSignificantExpDif) {
+ return extended_exponent_fpt(-that.value_, that.exponent_);
+ }
+ if (that.value_ == 0.0 ||
+ this->exponent_ > that.exponent_ + fea::kMaxSignificantExpDif) {
+ return *this;
+ }
+ if (this->exponent_ >= that.exponent_) {
+ exp_type exp_dif = this->exponent_ - that.exponent_;
+ fpt_type value = this->value_;
+ fea::set_exponent(value, exp_dif);
+ value -= that.value_;
+ return extended_exponent_fpt(value, that.exponent_);
+ } else {
+ exp_type exp_dif = that.exponent_ - this->exponent_;
+ fpt_type value = -that.value_;
+ fea::set_exponent(value, exp_dif);
+ value += this->value_;
+ return extended_exponent_fpt(value, this->exponent_);
+ }
+ }
+
+ extended_exponent_fpt operator*(const extended_exponent_fpt& that) const {
+ fpt_type value = this->value_ * that.value_;
+ exp_type exponent = this->exponent_ + that.exponent_;
+ return extended_exponent_fpt(value, exponent);
+ }
+
+ extended_exponent_fpt operator/(const extended_exponent_fpt& that) const {
+ fpt_type value = this->value_ / that.value_;
+ exp_type exponent = this->exponent_ - that.exponent_;
+ return extended_exponent_fpt(value, exponent);
+ }
+
+ extended_exponent_fpt& operator+=(const extended_exponent_fpt& that) {
+ return *this = *this + that;
+ }
+
+ extended_exponent_fpt& operator-=(const extended_exponent_fpt& that) {
+ return *this = *this - that;
+ }
+
+ extended_exponent_fpt& operator*=(const extended_exponent_fpt& that) {
+ return *this = *this * that;
+ }
+
+ extended_exponent_fpt& operator/=(const extended_exponent_fpt& that) {
+ return *this = *this / that;
+ }
+
+ extended_exponent_fpt sqrt() const {
+ fpt_type val = value_;
+ exp_type exp = exponent_;
+ if (exp & 1) {
+ val *= 2.0;
+ --exp;
+ }
+ return extended_exponent_fpt(get_sqrt(val), exp >> 1);
+ }
+
+ fpt_type d() const {
+ fpt_type ret_val = value_;
+ exp_type exp = exponent_;
+ if (ret_val == 0.0) {
+ return ret_val;
+ }
+ if (exp >= fea::kMaxExponent) {
+ ret_val = 1.0;
+ exp = fea::kMaxExponent;
+ } else if (exp <= fea::kMinExponent) {
+ ret_val = 1.0;
+ exp = fea::kMinExponent;
+ }
+ fea::set_exponent(ret_val, exp);
+ return ret_val;
+ }
+
+private:
+ fpt_type value_;
+ exp_type exponent_;
+};
+typedef extended_exponent_fpt<double> efpt64;
+
+template <typename _fpt>
+extended_exponent_fpt<_fpt> get_sqrt(const extended_exponent_fpt<_fpt>& that) {
+ return that.sqrt();
+}
+
+template <typename _fpt>
+bool is_pos(const extended_exponent_fpt<_fpt>& that) {
+ return that.is_pos();
+}
+
+template <typename _fpt>
+bool is_neg(const extended_exponent_fpt<_fpt>& that) {
+ return that.is_neg();
+}
+
+template <typename _fpt>
+bool is_zero(const extended_exponent_fpt<_fpt>& that) {
+ return that.is_zero();
+}
+
+template<size_t N>
+class extended_int {
+public:
+ static const uint64 kUInt64LowMask;
+ static const uint64 kUInt64HighMask;
+
+ extended_int() {}
+
+ extended_int(int32 that) {
+ if (that > 0) {
+ this->chunks_[0] = that;
+ this->count_ = 1;
+ } else if (that < 0) {
+ this->chunks_[0] = -that;
+ this->count_ = -1;
+ } else {
+ this->chunks_[0] = this->count_ = 0;
+ }
+ }
+
+ extended_int(int64 that) {
+ if (that > 0) {
+ this->chunks_[0] = static_cast<uint32>(that & kUInt64LowMask);
+ uint32 high = (that & kUInt64HighMask) >> 32;
+ if (high) {
+ this->chunks_[1] = high;
+ this->count_ = 2;
+ } else {
+ this->count_ = 1;
+ }
+ } else if (that < 0) {
+ that = -that;
+ this->chunks_[0] = static_cast<uint32>(that & kUInt64LowMask);
+ uint32 high = (that & kUInt64HighMask) >> 32;
+ if (high) {
+ this->chunks_[1] = high;
+ this->count_ = -2;
+ } else {
+ this->count_ = -1;
+ }
+ } else {
+ this->chunks_[0] = this->count_ = 0;
+ }
+ }
+
+ extended_int(const std::vector<uint32>& chunks, bool plus = true) {
+ this->count_ = static_cast<int32>(chunks.size());
+ for (size_t i = 0; i < chunks.size(); ++i) {
+ this->chunks_[i] = chunks[chunks.size() - i - 1];
+ }
+ if (!plus) {
+ this->count_ = -this->count_;
+ }
+ }
+
+ template <size_t M>
+ extended_int(const extended_int<M>& that) {
+ this->count_ = that.count();
+ memcpy(this->chunks_, that.chunks(), that.size() * sizeof(uint32));
+ }
+
+ extended_int<N>& operator=(int32 that) {
+ if (that > 0) {
+ this->chunks_[0] = that;
+ this->count_ = 1;
+ } else if (that < 0) {
+ this->chunks_[0] = -that;
+ this->count_ = -1;
+ } else {
+ this->chunks_[0] = this->count_ = 0;
+ }
+ return *this;
+ }
+
+ extended_int<N>& operator=(int64 that) {
+ if (that > 0) {
+ this->chunks_[0] = static_cast<uint32>(that & kUInt64LowMask);
+ uint32 high = (that & kUInt64HighMask) >> 32;
+ if (high) {
+ this->chunks_[1] = high;
+ this->count_ = 2;
+ } else {
+ this->count_ = 1;
+ }
+ } else if (that < 0) {
+ that = -that;
+ this->chunks_[0] = static_cast<uint32>(that & kUInt64LowMask);
+ uint32 high = (that & kUInt64HighMask) >> 32;
+ if (high) {
+ this->chunks_[1] = high;
+ this->count_ = -2;
+ } else {
+ this->count_ = -1;
+ }
+ } else {
+ this->chunks_[0] = this->count_ = 0;
+ }
+ return *this;
+ }
+
+ template <size_t M>
+ extended_int<N>& operator=(const extended_int<M>& that) {
+ this->count_ = that.count();
+ size_t sz = (std::min)(N, that.size());
+ memcpy(this->chunks_, that.chunks(), sz * sizeof(uint32));
+ return *this;
+ }
+
+ bool is_pos() const {
+ return this->count_ > 0;
+ }
+
+ bool is_neg() const {
+ return this->count_ < 0;
+ }
+
+ bool is_zero() const {
+ return this->count_ == 0;
+ }
+
+ template <size_t M>
+ bool operator==(const extended_int<M>& that) const {
+ if (this->count_ != that.count()) {
+ return false;
+ }
+ for (size_t i = 0; i < this->size(); ++i) {
+ if (this->chunks_[i] != that.chunks()[i]) {
+ return false;
+ }
+ }
+ return true;
+ }
+
+ template <size_t M>
+ bool operator!=(const extended_int<M>& that) const {
+ return !(*this == that);
+ }
+
+ template <size_t M>
+ bool operator<(const extended_int<M>& that) const {
+ if (this->count_ != that.count()) {
+ return this->count_ < that.count();
+ }
+ size_t i = this->size();
+ if (!i) {
+ return false;
+ }
+ do {
+ --i;
+ if (this->chunks_[i] != that.chunks()[i]) {
+ return (this->chunks_[i] < that.chunks()[i]) ^ (this->count_ < 0);
+ }
+ } while (i);
+ return false;
+ }
+
+ template <size_t M>
+ bool operator>(const extended_int<M>& that) const {
+ return that < *this;
+ }
+
+ template <size_t M>
+ bool operator<=(const extended_int<M>& that) const {
+ return !(that < *this);
+ }
+
+ template <size_t M>
+ bool operator>=(const extended_int<M>& that) const {
+ return !(*this < that);
+ }
+
+ extended_int<N> operator-() const {
+ extended_int<N> ret_val = *this;
+ ret_val.neg();
+ return ret_val;
+ }
+
+ void neg() {
+ this->count_ = -this->count_;
+ }
+
+ template <size_t M>
+ extended_int<(N>M?N:M)> operator+(const extended_int<M>& that) const {
+ extended_int<(N>M?N:M)> ret_val;
+ ret_val.add(*this, that);
+ return ret_val;
+ }
+
+ template <size_t N1, size_t N2>
+ void add(const extended_int<N1>& e1, const extended_int<N2>& e2) {
+ if (!e1.count()) {
+ *this = e2;
+ return;
+ }
+ if (!e2.count()) {
+ *this = e1;
+ return;
+ }
+ if ((e1.count() > 0) ^ (e2.count() > 0)) {
+ dif(e1.chunks(), e1.size(), e2.chunks(), e2.size());
+ } else {
+ add(e1.chunks(), e1.size(), e2.chunks(), e2.size());
+ }
+ if (e1.count() < 0) {
+ this->count_ = -this->count_;
+ }
+ }
+
+ template <size_t M>
+ extended_int<(N>M?N:M)> operator-(const extended_int<M>& that) const {
+ extended_int<(N>M?N:M)> ret_val;
+ ret_val.dif(*this, that);
+ return ret_val;
+ }
+
+ template <size_t N1, size_t N2>
+ void dif(const extended_int<N1>& e1, const extended_int<N2> &e2) {
+ if (!e1.count()) {
+ *this = e2;
+ this->count_ = -this->count_;
+ return;
+ }
+ if (!e2.count()) {
+ *this = e1;
+ return;
+ }
+ if ((e1.count() > 0) ^ (e2.count() > 0)) {
+ add(e1.chunks(), e1.size(), e2.chunks(), e2.size());
+ } else {
+ dif(e1.chunks(), e1.size(), e2.chunks(), e2.size());
+ }
+ if (e1.count() < 0) {
+ this->count_ = -this->count_;
+ }
+ }
+
+ extended_int<N> operator*(int32 that) const {
+ extended_int<N> temp(that);
+ return (*this) * temp;
+ }
+
+ extended_int<N> operator*(int64 that) const {
+ extended_int<N> temp(that);
+ return (*this) * temp;
+ }
+
+ template <size_t M>
+ extended_int<(N>M?N:M)> operator*(const extended_int<M>& that) const {
+ extended_int<(N>M?N:M)> ret_val;
+ ret_val.mul(*this, that);
+ return ret_val;
+ }
+
+ template <size_t N1, size_t N2>
+ void mul(const extended_int<N1>& e1, const extended_int<N2>& e2) {
+ if (!e1.count() || !e2.count()) {
+ this->chunks_[0] = this->count_ = 0;
+ return;
+ }
+ mul(e1.chunks(), e1.size(), e2.chunks(), e2.size());
+ if ((e1.count() > 0) ^ (e2.count() > 0)) {
+ this->count_ = -this->count_;
+ }
+ }
+
+ const uint32* chunks() const {
+ return chunks_;
+ }
+
+ int32 count() const {
+ return count_;
+ }
+
+ size_t size() const {
+ return (std::abs)(count_);
+ }
+
+ std::pair<fpt64, int64> p() const {
+ std::pair<fpt64, int64> ret_val(0, 0);
+ size_t sz = this->size();
+ if (!sz) {
+ return ret_val;
+ } else {
+ if (sz == 1) {
+ ret_val.first = static_cast<fpt64>(this->chunks_[0]);
+ ret_val.second = 0;
+ } else if (sz == 2) {
+ ret_val.first = static_cast<fpt64>(this->chunks_[1]) *
+ static_cast<fpt64>(0x100000000LL) +
+ static_cast<fpt64>(this->chunks_[0]);
+ ret_val.second = 0;
+ } else {
+ for (size_t i = 1; i <= 3; ++i) {
+ ret_val.first *= static_cast<fpt64>(0x100000000LL);
+ ret_val.first += static_cast<fpt64>(this->chunks_[sz - i]);
+ }
+ ret_val.second = (sz - 3) << 5;
+ }
+ }
+ if (this->count_ < 0) {
+ ret_val.first = -ret_val.first;
+ }
+ return ret_val;
+ }
+
+ fpt64 d() const {
+ std::pair<fpt64, int64> p = this->p();
+ extended_exponent_fpt<fpt64> efpt(p.first, p.second);
+ return efpt.d();
+ }
+
+private:
+ void add(const uint32* c1, size_t sz1, const uint32* c2, size_t sz2) {
+ if (sz1 < sz2) {
+ add(c2, sz2, c1, sz1);
+ return;
+ }
+ this->count_ = sz1;
+ uint64 temp = 0;
+ for (size_t i = 0; i < sz2; ++i) {
+ temp += static_cast<uint64>(c1[i]) +
+ static_cast<uint64>(c2[i]);
+ this->chunks_[i] = static_cast<uint32>(temp & kUInt64LowMask);
+ temp = (temp & kUInt64HighMask) >> 32;
+ }
+ for (size_t i = sz2; i < sz1; ++i) {
+ temp += static_cast<uint64>(c1[i]);
+ this->chunks_[i] = static_cast<uint32>(temp & kUInt64LowMask);
+ temp = (temp & kUInt64HighMask) >> 32;
+ }
+ if (temp && (this->count_ != N)) {
+ this->chunks_[this->count_] = static_cast<uint32>(temp & kUInt64LowMask);
+ ++this->count_;
+ }
+ }
+
+ void dif(const uint32* c1, size_t sz1, const uint32* c2, size_t sz2, bool rec = false) {
+ if (sz1 < sz2) {
+ dif(c2, sz2, c1, sz1, true);
+ this->count_ = -this->count_;
+ return;
+ } else if ((sz1 == sz2) && !rec) {
+ do {
+ --sz1;
+ if (c1[sz1] < c2[sz1]) {
+ ++sz1;
+ dif(c2, sz1, c1, sz1, true);
+ this->count_ = -this->count_;
+ return;
+ } else if (c1[sz1] > c2[sz1]) {
+ ++sz1;
+ break;
+ } else if (!sz1) {
+ this->chunks_[0] = this->count_ = 0;
+ return;
+ }
+ } while (sz1);
+ sz2 = sz1;
+ }
+ this->count_ = sz1-1;
+ bool flag = false;
+ for (size_t i = 0; i < sz2; ++i) {
+ this->chunks_[i] = c1[i] - c2[i] - (flag?1:0);
+ flag = (c1[i] < c2[i]) || ((c1[i] == c2[i]) && flag);
+ }
+ for (size_t i = sz2; i < sz1; ++i) {
+ this->chunks_[i] = c1[i] - (flag?1:0);
+ flag = !c1[i] && flag;
+ }
+ if (this->chunks_[this->count_]) {
+ ++this->count_;
+ }
+ }
+
+ void mul(const uint32* c1, size_t sz1, const uint32* c2, size_t sz2) {
+ uint64 cur = 0, nxt, tmp;
+ this->count_ = static_cast<int32>((std::min)(N, sz1 + sz2 - 1));
+ for (size_t shift = 0; shift < static_cast<size_t>(this->count_); ++shift) {
+ nxt = 0;
+ for (size_t first = 0; first <= shift; ++first) {
+ if (first >= sz1) {
+ break;
+ }
+ size_t second = shift - first;
+ if (second >= sz2) {
+ continue;
+ }
+ tmp = static_cast<uint64>(c1[first]) *
+ static_cast<uint64>(c2[second]);
+ cur += tmp & kUInt64LowMask;
+ nxt += (tmp & kUInt64HighMask) >> 32;
+ }
+ this->chunks_[shift] = static_cast<uint32>(cur & kUInt64LowMask);
+ cur = nxt + ((cur & kUInt64HighMask) >> 32);
+ }
+ if (cur && (this->count_ != N)) {
+ this->chunks_[this->count_] = static_cast<uint32>(cur & kUInt64LowMask);
+ ++this->count_;
+ }
+ }
+
+ uint32 chunks_[N];
+ int32 count_;
+};
+
+template <size_t N>
+const uint64 extended_int<N>::kUInt64LowMask = 0x00000000ffffffffULL;
+template <size_t N>
+const uint64 extended_int<N>::kUInt64HighMask = 0xffffffff00000000ULL;
+
+template <size_t N>
+bool is_pos(const extended_int<N>& that) {
+ return that.count() > 0;
+}
+
+template <size_t N>
+bool is_neg(const extended_int<N>& that) {
+ return that.count() < 0;
+}
+
+template <size_t N>
+bool is_zero(const extended_int<N>& that) {
+ return !that.count();
+}
+
+struct type_converter_fpt {
+ template <typename T>
+ fpt64 operator()(const T& that) const {
+ return static_cast<fpt64>(that);
+ }
+
+ template <size_t N>
+ fpt64 operator()(const extended_int<N>& that) const {
+ return that.d();
+ }
+
+ template <>
+ fpt64 operator()< extended_exponent_fpt<fpt64> >(const extended_exponent_fpt<fpt64>& that) const {
+ return that.d();
+ }
+};
+
+struct type_converter_efpt {
+ template <size_t N>
+ extended_exponent_fpt<fpt64> operator()(const extended_int<N>& that) const {
+ std::pair<fpt64, int64> p = that.p();
+ return extended_exponent_fpt<fpt64>(p.first, p.second);
+ }
+};
+
+template <typename T>
+struct voronoi_ctype_traits;
+
+template <>
+struct voronoi_ctype_traits<int32> {
+ typedef int32 int_type;
+ typedef uint32 uint_type;
+ typedef int64 int_x2_type;
+ typedef uint64 uint_x2_type;
+ typedef extended_int<128> big_int_type;
+ typedef fpt64 fpt_type;
+ typedef extended_exponent_fpt<fpt_type> efpt_type;
+ typedef ulp_comparison<fpt_type> ulp_cmp_type;
+ typedef type_converter_fpt to_fpt_converter_type;
+ typedef type_converter_efpt to_efpt_converter_type;
+};
+
+} // detail
+} // polygon
+} // boost
+
+#endif
\ No newline at end of file
Modified: sandbox/gtl/boost/polygon/detail/voronoi_predicates.hpp
==============================================================================
--- sandbox/gtl/boost/polygon/detail/voronoi_predicates.hpp (original)
+++ sandbox/gtl/boost/polygon/detail/voronoi_predicates.hpp 2012-01-23 17:22:43 EST (Mon, 23 Jan 2012)
@@ -1,14 +1,14 @@
// Boost.Polygon library detail/voronoi_predicates.hpp header file
-// Copyright Andrii Sydorchuk 2010-2011.
+// Copyright Andrii Sydorchuk 2010-2012.
// Distributed under the Boost Software License, Version 1.0.
// (See accompanying file LICENSE_1_0.txt or copy at
// http://www.boost.org/LICENSE_1_0.txt)
// See http://www.boost.org for updates, documentation, and revision history.
-#ifndef BOOST_POLYGON_VORONOI_CALC_UTILS
-#define BOOST_POLYGON_VORONOI_CALC_UTILS
+#ifndef BOOST_POLYGON_VORONOI_DETAIL_PREDICATES
+#define BOOST_POLYGON_VORONOI_DETAIL_PREDICATES
#include "voronoi_robust_fpt.hpp"
@@ -16,24 +16,21 @@
namespace polygon {
namespace detail {
-template <typename T>
-class voronoi_predicates;
-
// Predicate utilities. Operates with the coordinate types that could
// be converted to the 32-bit signed integer without precision loss.
-template <>
-class voronoi_predicates<int32> {
+template <typename CTYPE_TRAITS>
+class voronoi_predicates {
public:
- typedef int32 int_type;
- typedef uint32 uint_type;
- typedef int64 int_x2_type;
- typedef uint64 uint_x2_type;
- typedef eint4096 int_x128_type;
- typedef fpt64 fpt_type;
- typedef efpt64 efpt_type;
- typedef ulp_comparison<fpt_type> ulp_cmp_type;
- typedef type_converter_fpt to_fpt_converter;
- typedef type_converter_efpt to_efpt_converter;
+ typedef typename CTYPE_TRAITS::int_type int_type;
+ typedef typename CTYPE_TRAITS::uint_type uint_type;
+ typedef typename CTYPE_TRAITS::int_x2_type int_x2_type;
+ typedef typename CTYPE_TRAITS::uint_x2_type uint_x2_type;
+ typedef typename CTYPE_TRAITS::big_int_type big_int_type;
+ typedef typename CTYPE_TRAITS::fpt_type fpt_type;
+ typedef typename CTYPE_TRAITS::efpt_type efpt_type;
+ typedef typename CTYPE_TRAITS::ulp_cmp_type ulp_cmp_type;
+ typedef typename CTYPE_TRAITS::to_fpt_converter_type to_fpt_converter;
+ typedef typename CTYPE_TRAITS::to_efpt_converter_type to_efpt_converter;
static const unsigned int ULPS;
static const unsigned int ULPSx2;
@@ -529,7 +526,7 @@
typedef typename Site::point_type point_type;
typedef Site site_type;
typedef Circle circle_type;
- typedef robust_sqrt_expr<int_x128_type, efpt_type, to_efpt_converter> robust_sqrt_expr_type;
+ typedef robust_sqrt_expr<big_int_type, efpt_type, to_efpt_converter> robust_sqrt_expr_type;
void ppp(const site_type &site1,
const site_type &site2,
@@ -538,7 +535,7 @@
bool recompute_c_x = true,
bool recompute_c_y = true,
bool recompute_lower_x = true) {
- int_x128_type dif_x[3], dif_y[3], sum_x[2], sum_y[2];
+ big_int_type dif_x[3], dif_y[3], sum_x[2], sum_y[2];
dif_x[0] = static_cast<int_x2_type>(site1.x()) -
static_cast<int_x2_type>(site2.x());
dif_x[1] = static_cast<int_x2_type>(site2.x()) -
@@ -560,16 +557,16 @@
sum_y[1] = static_cast<int_x2_type>(site2.y()) +
static_cast<int_x2_type>(site3.y());
fpt_type inv_denom = to_fpt(0.5) / to_fpt(dif_x[0] * dif_y[1] - dif_x[1] * dif_y[0]);
- int_x128_type numer1 = dif_x[0] * sum_x[0] + dif_y[0] * sum_y[0];
- int_x128_type numer2 = dif_x[1] * sum_x[1] + dif_y[1] * sum_y[1];
+ big_int_type numer1 = dif_x[0] * sum_x[0] + dif_y[0] * sum_y[0];
+ big_int_type numer2 = dif_x[1] * sum_x[1] + dif_y[1] * sum_y[1];
if (recompute_c_x || recompute_lower_x) {
- int_x128_type c_x = numer1 * dif_y[1] - numer2 * dif_y[0];
+ big_int_type c_x = numer1 * dif_y[1] - numer2 * dif_y[0];
circle.x(to_fpt(c_x) * inv_denom);
if (recompute_lower_x) {
// Evaluate radius of the circle.
- int_x128_type sqr_r = (dif_x[0] * dif_x[0] + dif_y[0] * dif_y[0]) *
+ big_int_type sqr_r = (dif_x[0] * dif_x[0] + dif_y[0] * dif_y[0]) *
(dif_x[1] * dif_x[1] + dif_y[1] * dif_y[1]) *
(dif_x[2] * dif_x[2] + dif_y[2] * dif_y[2]);
fpt_type r = get_sqrt(to_fpt(sqr_r));
@@ -584,7 +581,7 @@
circle.lower_x(circle.x() - r * inv_denom);
}
} else {
- int_x128_type numer = c_x * c_x - sqr_r;
+ big_int_type numer = c_x * c_x - sqr_r;
fpt_type lower_x = to_fpt(numer) * inv_denom /
(to_fpt(c_x) + r);
circle.lower_x(lower_x);
@@ -593,7 +590,7 @@
}
if (recompute_c_y) {
- int_x128_type c_y = numer2 * dif_x[0] - numer1 * dif_x[1];
+ big_int_type c_y = numer2 * dif_x[0] - numer1 * dif_x[1];
circle.y(to_fpt(c_y) * inv_denom);
}
}
@@ -607,38 +604,38 @@
bool recompute_c_x = true,
bool recompute_c_y = true,
bool recompute_lower_x = true) {
- int_x128_type cA[4], cB[4];
- int_x128_type line_a = static_cast<int_x2_type>(site3.point1(true).y()) -
+ big_int_type cA[4], cB[4];
+ big_int_type line_a = static_cast<int_x2_type>(site3.point1(true).y()) -
static_cast<int_x2_type>(site3.point0(true).y());
- int_x128_type line_b = static_cast<int_x2_type>(site3.point0(true).x()) -
+ big_int_type line_b = static_cast<int_x2_type>(site3.point0(true).x()) -
static_cast<int_x2_type>(site3.point1(true).x());
- int_x128_type segm_len = line_a * line_a + line_b * line_b;
- int_x128_type vec_x = static_cast<int_x2_type>(site2.y()) -
+ big_int_type segm_len = line_a * line_a + line_b * line_b;
+ big_int_type vec_x = static_cast<int_x2_type>(site2.y()) -
static_cast<int_x2_type>(site1.y());
- int_x128_type vec_y = static_cast<int_x2_type>(site1.x()) -
+ big_int_type vec_y = static_cast<int_x2_type>(site1.x()) -
static_cast<int_x2_type>(site2.x());
- int_x128_type sum_x = static_cast<int_x2_type>(site1.x()) +
+ big_int_type sum_x = static_cast<int_x2_type>(site1.x()) +
static_cast<int_x2_type>(site2.x());
- int_x128_type sum_y = static_cast<int_x2_type>(site1.y()) +
+ big_int_type sum_y = static_cast<int_x2_type>(site1.y()) +
static_cast<int_x2_type>(site2.y());
- int_x128_type teta = line_a * vec_x + line_b * vec_y;
- int_x128_type denom = vec_x * line_b - vec_y * line_a;
+ big_int_type teta = line_a * vec_x + line_b * vec_y;
+ big_int_type denom = vec_x * line_b - vec_y * line_a;
- int_x128_type dif0 = static_cast<int_x2_type>(site3.point1().y()) -
+ big_int_type dif0 = static_cast<int_x2_type>(site3.point1().y()) -
static_cast<int_x2_type>(site1.y());
- int_x128_type dif1 = static_cast<int_x2_type>(site1.x()) -
+ big_int_type dif1 = static_cast<int_x2_type>(site1.x()) -
static_cast<int_x2_type>(site3.point1().x());
- int_x128_type A = line_a * dif1 - line_b * dif0;
+ big_int_type A = line_a * dif1 - line_b * dif0;
dif0 = static_cast<int_x2_type>(site3.point1().y()) -
static_cast<int_x2_type>(site2.y());
dif1 = static_cast<int_x2_type>(site2.x()) -
static_cast<int_x2_type>(site3.point1().x());
- int_x128_type B = line_a * dif1 - line_b * dif0;
- int_x128_type sum_AB = A + B;
+ big_int_type B = line_a * dif1 - line_b * dif0;
+ big_int_type sum_AB = A + B;
if (is_zero(denom)) {
- int_x128_type numer = teta * teta - sum_AB * sum_AB;
- int_x128_type denom = teta * sum_AB;
+ big_int_type numer = teta * teta - sum_AB * sum_AB;
+ big_int_type 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;
@@ -658,7 +655,7 @@
return;
}
- int_x128_type det = (teta * teta + denom * denom) * A * B * 4;
+ big_int_type det = (teta * teta + denom * denom) * A * B * 4;
fpt_type inv_denom_sqr = to_fpt(1.0) / to_fpt(denom);
inv_denom_sqr *= inv_denom_sqr;
@@ -704,7 +701,7 @@
bool recompute_c_x = true,
bool recompute_c_y = true,
bool recompute_lower_x = true) {
- int_x128_type a[2], b[2], c[2], cA[4], cB[4];
+ big_int_type a[2], b[2], c[2], cA[4], cB[4];
const point_type &segm_start1 = site2.point1(true);
const point_type &segm_end1 = site2.point0(true);
const point_type &segm_start2 = site3.point0(true);
@@ -717,18 +714,18 @@
static_cast<int_x2_type>(segm_start2.x());
b[1] = static_cast<int_x2_type>(segm_end2.y()) -
static_cast<int_x2_type>(segm_start2.y());
- int_x128_type orientation = a[1] * b[0] - a[0] * b[1];
+ big_int_type orientation = a[1] * b[0] - a[0] * b[1];
if (is_zero(orientation)) {
fpt_type denom = to_fpt(2.0) * to_fpt(a[0] * a[0] + b[0] * b[0]);
c[0] = b[0] * (static_cast<int_x2_type>(segm_start2.x()) -
static_cast<int_x2_type>(segm_start1.x())) -
a[0] * (static_cast<int_x2_type>(segm_start2.y()) -
static_cast<int_x2_type>(segm_start1.y()));
- int_x128_type dx = a[0] * (static_cast<int_x2_type>(site1.y()) -
+ big_int_type dx = a[0] * (static_cast<int_x2_type>(site1.y()) -
static_cast<int_x2_type>(segm_start1.y())) -
b[0] * (static_cast<int_x2_type>(site1.x()) -
static_cast<int_x2_type>(segm_start1.x()));
- int_x128_type dy = b[0] * (static_cast<int_x2_type>(site1.x()) -
+ big_int_type dy = b[0] * (static_cast<int_x2_type>(site1.x()) -
static_cast<int_x2_type>(segm_start2.x())) -
a[0] * (static_cast<int_x2_type>(site1.y()) -
static_cast<int_x2_type>(segm_start2.y()));
@@ -774,10 +771,10 @@
a[0] * segm_end1.y();
c[1] = a[1] * segm_end2.y() -
b[1] * segm_end2.x();
- int_x128_type ix = a[0] * c[1] + a[1] * c[0];
- int_x128_type iy = b[0] * c[1] + b[1] * c[0];
- int_x128_type dx = ix - orientation * site1.x();
- int_x128_type dy = iy - orientation * site1.y();
+ big_int_type ix = a[0] * c[1] + a[1] * c[0];
+ big_int_type iy = b[0] * c[1] + b[1] * c[0];
+ big_int_type dx = ix - orientation * site1.x();
+ big_int_type dy = iy - orientation * site1.y();
if (is_zero(dx) && is_zero(dy)) {
fpt_type denom = to_fpt(orientation);
fpt_type c_x = to_fpt(ix) / denom;
@@ -786,7 +783,7 @@
return;
}
- int_x128_type sign = ((point_index == 2) ? 1 : -1) * (is_neg(orientation) ? 1 : -1);
+ big_int_type sign = ((point_index == 2) ? 1 : -1) * (is_neg(orientation) ? 1 : -1);
cA[0] = a[1] * -dx + b[1] * -dy;
cA[1] = a[0] * -dx + b[0] * -dy;
cA[2] = sign;
@@ -795,14 +792,14 @@
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;
- fpt_type temp = to_fpt(sqrt_expr_evaluator_pss4<int_x128_type, efpt_type>(cA, cB));
+ fpt_type temp = to_fpt(sqrt_expr_evaluator_pss4<big_int_type, efpt_type>(cA, cB));
fpt_type denom = temp * to_fpt(orientation);
if (recompute_c_y) {
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;
- fpt_type cy = to_fpt(sqrt_expr_evaluator_pss4<int_x128_type, efpt_type>(cA, cB));
+ fpt_type cy = to_fpt(sqrt_expr_evaluator_pss4<big_int_type, efpt_type>(cA, cB));
c_event.y(cy / denom);
}
@@ -812,13 +809,13 @@
cA[2] = ix * sign;
if (recompute_c_x) {
- fpt_type cx = to_fpt(sqrt_expr_evaluator_pss4<int_x128_type, efpt_type>(cA, cB));
+ fpt_type cx = to_fpt(sqrt_expr_evaluator_pss4<big_int_type, efpt_type>(cA, cB));
c_event.x(cx / denom);
}
if (recompute_lower_x) {
cA[3] = orientation * (dx * dx + dy * dy) * (is_neg(temp) ? -1 : 1);
- fpt_type lower_x = to_fpt(sqrt_expr_evaluator_pss4<int_x128_type, efpt_type>(cA, cB));
+ fpt_type lower_x = to_fpt(sqrt_expr_evaluator_pss4<big_int_type, efpt_type>(cA, cB));
c_event.lower_x(lower_x / denom);
}
}
@@ -832,7 +829,7 @@
bool recompute_c_x = true,
bool recompute_c_y = true,
bool recompute_lower_x = true) {
- int_x128_type a[3], b[3], c[3], cA[4], cB[4];
+ big_int_type a[3], b[3], c[3], cA[4], cB[4];
// cA - corresponds to the cross product.
// cB - corresponds to the squared length.
a[0] = static_cast<int_x2_type>(site1.x1(true)) -
@@ -1381,12 +1378,16 @@
}
};
-const unsigned int voronoi_predicates<int>::ULPS = 64;
-const unsigned int voronoi_predicates<int>::ULPSx2 = 128;
-const voronoi_predicates<int>::fpt_type voronoi_predicates<int>::fULPS =
- voronoi_predicates<int>::ULPS;
-const voronoi_predicates<int>::fpt_type voronoi_predicates<int>::fULPSx2 =
- voronoi_predicates<int>::ULPSx2;
+template <typename CTYPE_TRAITS>
+const unsigned int voronoi_predicates<CTYPE_TRAITS>::ULPS = 64;
+template <typename CTYPE_TRAITS>
+const unsigned int voronoi_predicates<CTYPE_TRAITS>::ULPSx2 = 128;
+template <typename CTYPE_TRAITS>
+const typename voronoi_predicates<CTYPE_TRAITS>::fpt_type voronoi_predicates<CTYPE_TRAITS>::fULPS =
+ voronoi_predicates<CTYPE_TRAITS>::ULPS;
+template <typename CTYPE_TRAITS>
+const typename voronoi_predicates<CTYPE_TRAITS>::fpt_type voronoi_predicates<CTYPE_TRAITS>::fULPSx2 =
+ voronoi_predicates<CTYPE_TRAITS>::ULPSx2;
} // detail
} // polygon
Modified: sandbox/gtl/boost/polygon/detail/voronoi_robust_fpt.hpp
==============================================================================
--- sandbox/gtl/boost/polygon/detail/voronoi_robust_fpt.hpp (original)
+++ sandbox/gtl/boost/polygon/detail/voronoi_robust_fpt.hpp 2012-01-23 17:22:43 EST (Mon, 23 Jan 2012)
@@ -1,19 +1,17 @@
// Boost.Polygon library detail/voronoi_robust_fpt.hpp header file
-// Copyright Andrii Sydorchuk 2010-2011.
+// Copyright Andrii Sydorchuk 2010-2012.
// Distributed under the Boost Software License, Version 1.0.
// (See accompanying file LICENSE_1_0.txt or copy at
// http://www.boost.org/LICENSE_1_0.txt)
// See http://www.boost.org for updates, documentation, and revision history.
-#ifndef BOOST_POLYGON_VORONOI_ROBUST_FPT
-#define BOOST_POLYGON_VORONOI_ROBUST_FPT
+#ifndef BOOST_POLYGON_VORONOI_DETAIL_ROBUST_FPT
+#define BOOST_POLYGON_VORONOI_DETAIL_ROBUST_FPT
#include <cmath>
-#include <boost/cstdint.hpp>
-
// Geometry predicates with floating-point variables usually require
// high-precision predicates to retrieve the correct result.
// Epsilon robust predicates give the result within some epsilon relative
@@ -56,1181 +54,467 @@
namespace boost {
namespace polygon {
namespace detail {
- typedef boost::int32_t int32;
- typedef boost::int64_t int64;
- typedef boost::uint32_t uint32;
- typedef boost::uint64_t uint64;
- typedef double fpt64;
-
- // 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 reinterpretations differ in not more than maxUlps units.
- template <typename _fpt>
- struct ulp_comparison;
-
- template <>
- struct ulp_comparison<fpt64> {
- enum kResult {
- LESS = -1,
- EQUAL = 0,
- MORE = 1
- };
-
- kResult operator()(fpt64 a, fpt64 b, unsigned int maxUlps) const {
- uint64 ll_a, ll_b;
-
- // Reinterpret double bits as 64-bit signed integer.
- memcpy(&ll_a, &a, sizeof(fpt64));
- memcpy(&ll_b, &b, sizeof(fpt64));
-
- // 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) ? EQUAL : LESS;
- return (ll_b - ll_a <= maxUlps) ? EQUAL : MORE;
- }
- };
- template <typename T>
- T get_sqrt(const T& that) {
- return (std::sqrt)(that);
- }
+template <typename T>
+T get_sqrt(const T& that) {
+ return (std::sqrt)(that);
+}
- template <typename T>
- bool is_pos(const T& that) {
- return that > 0;
- }
+template <typename T>
+bool is_pos(const T& that) {
+ return that > 0;
+}
- template <typename T>
- bool is_neg(const T& that) {
- return that < 0;
- }
+template <typename T>
+bool is_neg(const T& that) {
+ return that < 0;
+}
- template <typename T>
- bool is_zero(const T& that) {
- return that == 0;
- }
+template <typename T>
+bool is_zero(const T& that) {
+ return that == 0;
+}
- template <typename _fpt>
- class robust_fpt {
- public:
- typedef _fpt floating_point_type;
- typedef _fpt relative_error_type;
+template <typename _fpt>
+class robust_fpt {
+public:
+ typedef _fpt floating_point_type;
+ typedef _fpt relative_error_type;
- // Rounding error is at most 1 EPS.
- static const relative_error_type ROUNDING_ERROR;
+ // Rounding error is at most 1 EPS.
+ static const relative_error_type ROUNDING_ERROR;
- robust_fpt() : fpv_(0.0), re_(0.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) {}
+ robust_fpt() : fpv_(0.0), re_(0.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_; }
-
- robust_fpt& operator=(const robust_fpt &that) {
- this->fpv_ = that.fpv_;
- this->re_ = that.re_;
- return *this;
- }
+ floating_point_type fpv() const { return fpv_; }
+ relative_error_type re() const { return re_; }
+ relative_error_type ulp() const { return re_; }
- bool has_pos_value() const {
- return is_pos(fpv_);
- }
+ robust_fpt& operator=(const robust_fpt &that) {
+ this->fpv_ = that.fpv_;
+ this->re_ = that.re_;
+ return *this;
+ }
- bool has_neg_value() const {
- return is_neg(fpv_);
- }
+ bool has_pos_value() const {
+ return is_pos(fpv_);
+ }
- bool has_zero_value() const {
- return is_zero(fpv_);
- }
+ bool has_neg_value() const {
+ return is_neg(fpv_);
+ }
- robust_fpt operator-() const {
- return robust_fpt(-fpv_, re_);
- }
+ bool has_zero_value() const {
+ return is_zero(fpv_);
+ }
+
+ robust_fpt operator-() const {
+ return robust_fpt(-fpv_, re_);
+ }
- robust_fpt& operator+=(const robust_fpt &that) {
- floating_point_type fpv = this->fpv_ + that.fpv_;
- if ((!is_neg(this->fpv_) && !is_neg(that.fpv_)) ||
- (!is_pos(this->fpv_) && !is_pos(that.fpv_)))
- this->re_ = (std::max)(this->re_, that.re_) + ROUNDING_ERROR;
- else {
- floating_point_type temp =
- (this->fpv_ * this->re_ - that.fpv_ * that.re_) / fpv;
- if (is_neg(temp)) {
- temp = -temp;
- }
- this->re_ = temp + ROUNDING_ERROR;
+ robust_fpt& operator+=(const robust_fpt &that) {
+ floating_point_type fpv = this->fpv_ + that.fpv_;
+ if ((!is_neg(this->fpv_) && !is_neg(that.fpv_)) ||
+ (!is_pos(this->fpv_) && !is_pos(that.fpv_)))
+ this->re_ = (std::max)(this->re_, that.re_) + ROUNDING_ERROR;
+ else {
+ floating_point_type temp =
+ (this->fpv_ * this->re_ - that.fpv_ * that.re_) / fpv;
+ if (is_neg(temp)) {
+ temp = -temp;
}
- this->fpv_ = fpv;
- return *this;
+ this->re_ = temp + ROUNDING_ERROR;
}
+ this->fpv_ = fpv;
+ return *this;
+ }
- robust_fpt& operator-=(const robust_fpt &that) {
- floating_point_type fpv = this->fpv_ - that.fpv_;
- if ((!is_neg(this->fpv_) && !is_pos(that.fpv_)) ||
- (!is_pos(this->fpv_) && !is_neg(that.fpv_)))
- this->re_ = (std::max)(this->re_, that.re_) + ROUNDING_ERROR;
- else {
- floating_point_type temp =
- (this->fpv_ * this->re_ + that.fpv_ * that.re_) / fpv;
- if (is_neg(temp)) {
- temp = -temp;
- }
- this->re_ = temp + ROUNDING_ERROR;
+ robust_fpt& operator-=(const robust_fpt &that) {
+ floating_point_type fpv = this->fpv_ - that.fpv_;
+ if ((!is_neg(this->fpv_) && !is_pos(that.fpv_)) ||
+ (!is_pos(this->fpv_) && !is_neg(that.fpv_)))
+ this->re_ = (std::max)(this->re_, that.re_) + ROUNDING_ERROR;
+ else {
+ floating_point_type temp =
+ (this->fpv_ * this->re_ + that.fpv_ * that.re_) / fpv;
+ if (is_neg(temp)) {
+ temp = -temp;
}
- this->fpv_ = fpv;
- return *this;
+ this->re_ = 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) {
- 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 ((!is_neg(this->fpv_) && !is_neg(that.fpv_)) ||
- (!is_pos(this->fpv_) && !is_pos(that.fpv_)))
- re = (std::max)(this->re_, that.re_) + ROUNDING_ERROR;
- else {
- floating_point_type temp =
- (this->fpv_ * this->re_ - that.fpv_ * that.re_) / fpv;
- if (is_neg(temp)) {
- temp = -temp;
- }
- re = temp + ROUNDING_ERROR;
+ robust_fpt operator+(const robust_fpt &that) const {
+ floating_point_type fpv = this->fpv_ + that.fpv_;
+ relative_error_type re;
+ if ((!is_neg(this->fpv_) && !is_neg(that.fpv_)) ||
+ (!is_pos(this->fpv_) && !is_pos(that.fpv_)))
+ re = (std::max)(this->re_, that.re_) + ROUNDING_ERROR;
+ else {
+ floating_point_type temp =
+ (this->fpv_ * this->re_ - that.fpv_ * that.re_) / fpv;
+ if (is_neg(temp)) {
+ temp = -temp;
}
- return robust_fpt(fpv, re);
+ re = 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 ((!is_neg(this->fpv_) && !is_pos(that.fpv_)) ||
- (!is_pos(this->fpv_) && !is_neg(that.fpv_)))
- re = (std::max)(this->re_, that.re_) + ROUNDING_ERROR;
- else {
- floating_point_type temp =
- (this->fpv_ * this->re_ + that.fpv_ * that.re_) / fpv;
- if (is_neg(temp)) {
- temp = -temp;
- }
- re = temp + ROUNDING_ERROR;
+ robust_fpt operator-(const robust_fpt &that) const {
+ floating_point_type fpv = this->fpv_ - that.fpv_;
+ relative_error_type re;
+ if ((!is_neg(this->fpv_) && !is_pos(that.fpv_)) ||
+ (!is_pos(this->fpv_) && !is_neg(that.fpv_)))
+ re = (std::max)(this->re_, that.re_) + ROUNDING_ERROR;
+ else {
+ floating_point_type temp =
+ (this->fpv_ * this->re_ + that.fpv_ * that.re_) / fpv;
+ if (is_neg(temp)) {
+ temp = -temp;
}
- 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);
+ re = temp + ROUNDING_ERROR;
}
-
- 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 sqrt() const {
- return robust_fpt(get_sqrt(fpv_),
- re_ * static_cast<relative_error_type>(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;
-
- template <typename T>
- robust_fpt<T> get_sqrt(const robust_fpt<T>& that) {
- return that.sqrt();
+ return robust_fpt(fpv, re);
}
- template <typename T>
- bool is_pos(const robust_fpt<T>& that) {
- return that.has_pos_value();
+ 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);
}
- template <typename T>
- bool is_neg(const robust_fpt<T>& that) {
- return that.has_neg_value();
+ 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);
}
- template <typename T>
- bool is_zero(const robust_fpt<T>& that) {
- return that.has_zero_value();
+ robust_fpt sqrt() const {
+ return robust_fpt(get_sqrt(fpv_),
+ re_ * static_cast<relative_error_type>(0.5) +
+ ROUNDING_ERROR);
}
- // robust_dif consists of two not negative values: value1 and value2.
- // The resulting expression is equal to the value1 - value2.
- // 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. The structure implicitly avoids difference computation.
- template <typename T>
- class robust_dif {
- public:
- robust_dif() :
- positive_sum_(0),
- negative_sum_(0) {}
+private:
+ floating_point_type fpv_;
+ relative_error_type re_;
+};
- robust_dif(const T &value) :
- positive_sum_((value>0)?value:0),
- negative_sum_((value<0)?-value:0) {}
+template <typename T>
+const typename robust_fpt<T>::relative_error_type
+ robust_fpt<T>::ROUNDING_ERROR = 1;
- robust_dif(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_;
- }
+template <typename T>
+robust_fpt<T> get_sqrt(const robust_fpt<T>& that) {
+ return that.sqrt();
+}
- T neg() const {
- return negative_sum_;
- }
-
- robust_dif<T> operator-() const {
- return robust_dif(negative_sum_, positive_sum_);
- }
-
- robust_dif<T> &operator+=(const T &val) {
- if (!is_neg(val))
- positive_sum_ += val;
- else
- negative_sum_ -= val;
- return *this;
- }
-
- robust_dif<T> &operator+=(const robust_dif<T> &that) {
- positive_sum_ += that.positive_sum_;
- negative_sum_ += that.negative_sum_;
- return *this;
- }
+template <typename T>
+bool is_pos(const robust_fpt<T>& that) {
+ return that.has_pos_value();
+}
- robust_dif<T> &operator-=(const T &val) {
- if (!is_neg(val))
- negative_sum_ += val;
- else
- positive_sum_ -= val;
- return *this;
- }
+template <typename T>
+bool is_neg(const robust_fpt<T>& that) {
+ return that.has_neg_value();
+}
- robust_dif<T> &operator-=(const robust_dif<T> &that) {
- positive_sum_ += that.negative_sum_;
- negative_sum_ += that.positive_sum_;
- return *this;
- }
+template <typename T>
+bool is_zero(const robust_fpt<T>& that) {
+ return that.has_zero_value();
+}
- robust_dif<T> &operator*=(const T &val) {
- if (!is_neg(val)) {
- positive_sum_ *= val;
- negative_sum_ *= val;
- } else {
- positive_sum_ *= -val;
- negative_sum_ *= -val;
- swap();
- }
- return *this;
- }
+// robust_dif consists of two not negative values: value1 and value2.
+// The resulting expression is equal to the value1 - value2.
+// 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. The structure implicitly avoids difference computation.
+template <typename T>
+class robust_dif {
+public:
+ robust_dif() :
+ positive_sum_(0),
+ negative_sum_(0) {}
- robust_dif<T> &operator*=(const robust_dif<T> &that) {
- T positive_sum = this->positive_sum_ * that.positive_sum_ +
- this->negative_sum_ * that.negative_sum_;
- T negative_sum = this->positive_sum_ * that.negative_sum_ +
- this->negative_sum_ * that.positive_sum_;
- positive_sum_ = positive_sum;
- negative_sum_ = negative_sum;
- return *this;
- }
+ robust_dif(const T &value) :
+ positive_sum_((value>0)?value:0),
+ negative_sum_((value<0)?-value:0) {}
- robust_dif<T> &operator/=(const T &val) {
- if (!is_neg(val)) {
- positive_sum_ /= val;
- negative_sum_ /= val;
- } else {
- positive_sum_ /= -val;
- negative_sum_ /= -val;
- swap();
- }
- return *this;
- }
+ robust_dif(const T &pos, const T &neg) :
+ positive_sum_(pos),
+ negative_sum_(neg) {}
- private:
- void swap() {
- (std::swap)(positive_sum_, negative_sum_);
- }
+ T dif() const {
+ return positive_sum_ - negative_sum_;
+ }
- T positive_sum_;
- T negative_sum_;
- };
-
- template<typename T>
- robust_dif<T> operator+(const robust_dif<T>& lhs,
- const robust_dif<T>& rhs) {
- return robust_dif<T>(lhs.pos() + rhs.pos(),
- lhs.neg() + rhs.neg());
+ T pos() const {
+ return positive_sum_;
}
- template<typename T>
- robust_dif<T> operator+(const robust_dif<T>& lhs, const T& rhs) {
- if (!is_neg(rhs)) {
- return robust_dif<T>(lhs.pos() + rhs, lhs.neg());
- } else {
- return robust_dif<T>(lhs.pos(), lhs.neg() - rhs);
- }
+ T neg() const {
+ return negative_sum_;
}
- template<typename T>
- robust_dif<T> operator+(const T& lhs, const robust_dif<T>& rhs) {
- if (!is_neg(lhs)) {
- return robust_dif<T>(lhs + rhs.pos(), rhs.neg());
- } else {
- return robust_dif<T>(rhs.pos(), rhs.neg() - lhs);
- }
+ robust_dif<T> operator-() const {
+ return robust_dif(negative_sum_, positive_sum_);
}
- template<typename T>
- robust_dif<T> operator-(const robust_dif<T>& lhs,
- const robust_dif<T>& rhs) {
- return robust_dif<T>(lhs.pos() + rhs.neg(), lhs.neg() + rhs.pos());
+ robust_dif<T> &operator+=(const T &val) {
+ if (!is_neg(val))
+ positive_sum_ += val;
+ else
+ negative_sum_ -= val;
+ return *this;
}
- template<typename T>
- robust_dif<T> operator-(const robust_dif<T>& lhs, const T& rhs) {
- if (!is_neg(rhs)) {
- return robust_dif<T>(lhs.pos(), lhs.neg() + rhs);
- } else {
- return robust_dif<T>(lhs.pos() - rhs, lhs.neg());
- }
+ robust_dif<T> &operator+=(const robust_dif<T> &that) {
+ positive_sum_ += that.positive_sum_;
+ negative_sum_ += that.negative_sum_;
+ return *this;
}
- template<typename T>
- robust_dif<T> operator-(const T& lhs, const robust_dif<T>& rhs) {
- if (!is_neg(lhs)) {
- return robust_dif<T>(lhs + rhs.neg(), rhs.pos());
- } else {
- return robust_dif<T>(rhs.neg(), rhs.pos() - lhs);
- }
+ robust_dif<T> &operator-=(const T &val) {
+ if (!is_neg(val))
+ negative_sum_ += val;
+ else
+ positive_sum_ -= val;
+ return *this;
}
- template<typename T>
- robust_dif<T> operator*(const robust_dif<T>& lhs,
- const robust_dif<T>& rhs) {
- T res_pos = lhs.pos() * rhs.pos() + lhs.neg() * rhs.neg();
- T res_neg = lhs.pos() * rhs.neg() + lhs.neg() * rhs.pos();
- return robust_dif<T>(res_pos, res_neg);
+ robust_dif<T> &operator-=(const robust_dif<T> &that) {
+ positive_sum_ += that.negative_sum_;
+ negative_sum_ += that.positive_sum_;
+ return *this;
}
- template<typename T>
- robust_dif<T> operator*(const robust_dif<T>& lhs, const T& val) {
+ robust_dif<T> &operator*=(const T &val) {
if (!is_neg(val)) {
- return robust_dif<T>(lhs.pos() * val, lhs.neg() * val);
+ positive_sum_ *= val;
+ negative_sum_ *= val;
} else {
- return robust_dif<T>(-lhs.neg() * val, -lhs.pos() * val);
+ positive_sum_ *= -val;
+ negative_sum_ *= -val;
+ swap();
}
+ return *this;
}
- template<typename T>
- robust_dif<T> operator*(const T& val, const robust_dif<T>& rhs) {
- if (!is_neg(val)) {
- return robust_dif<T>(val * rhs.pos(), val * rhs.neg());
- } else {
- return robust_dif<T>(-val * rhs.neg(), -val * rhs.pos());
- }
+ robust_dif<T> &operator*=(const robust_dif<T> &that) {
+ T positive_sum = this->positive_sum_ * that.positive_sum_ +
+ this->negative_sum_ * that.negative_sum_;
+ T negative_sum = this->positive_sum_ * that.negative_sum_ +
+ this->negative_sum_ * that.positive_sum_;
+ positive_sum_ = positive_sum;
+ negative_sum_ = negative_sum;
+ return *this;
}
- template<typename T>
- robust_dif<T> operator/(const robust_dif<T>& lhs, const T& val) {
+ robust_dif<T> &operator/=(const T &val) {
if (!is_neg(val)) {
- return robust_dif<T>(lhs.pos() / val, lhs.neg() / val);
+ positive_sum_ /= val;
+ negative_sum_ /= val;
} else {
- return robust_dif<T>(-lhs.neg() / val, -lhs.pos() / val);
- }
- }
-
- // Manages exponent of the floating-point value.
- template <typename _fpt>
- struct fpt_exponent_accessor;
-
- template <>
- class fpt_exponent_accessor<fpt64> {
- public:
- static const int64 kExponentMask;
- static const int64 kSignedMantissaMask;
- static const int64 kMinExponent;
- static const int64 kMaxExponent;
- static const int64 kMaxSignificantExpDif;
-
- static int64 set_exponent(fpt64& value, int64 exponent) {
- int64 bits;
- memcpy(&bits, &value, sizeof(fpt64));
- int64 exp = ((bits & kExponentMask) >> 52) - 1023;
- if (exp == exponent) {
- return exp;
- }
- bits = (bits & kSignedMantissaMask) | ((exponent + 1023) << 52);
- memcpy(&value, &bits, sizeof(fpt64));
- return exp;
- }
- };
-
- const int64 fpt_exponent_accessor<fpt64>::kExponentMask = 0x7ff0000000000000LL;
- const int64 fpt_exponent_accessor<fpt64>::kSignedMantissaMask = 0x800fffffffffffffLL;
- const int64 fpt_exponent_accessor<fpt64>::kMinExponent = -1023LL;
- const int64 fpt_exponent_accessor<fpt64>::kMaxExponent = 1024LL;
- const int64 fpt_exponent_accessor<fpt64>::kMaxSignificantExpDif = 54;
-
- // Allows to extend floating-point type exponent boundaries to the 64 bit
- // integer range. This class does not handle division by zero, subnormal
- // numbers or NaNs.
- template <typename _fpt>
- class extended_exponent_fpt {
- public:
- typedef _fpt fpt_type;
- typedef int64 exp_type;
- typedef fpt_exponent_accessor<fpt_type> fea;
-
- explicit extended_exponent_fpt(fpt_type value) {
- if (value == 0.0) {
- exponent_ = 0;
- value_ = 0.0;
- } else {
- exponent_ = fea::set_exponent(value, 0);
- value_ = value;
- }
- }
-
- extended_exponent_fpt(fpt_type value, exp_type exponent) {
- if (value == 0.0) {
- exponent_ = 0;
- value_ = 0.0;
- } else {
- exponent_ = fea::set_exponent(value, 0) + exponent;
- value_ = value;
- }
- }
-
- bool is_pos() const {
- return value_ > 0;
- }
-
- bool is_neg() const {
- return value_ < 0;
- }
-
- bool is_zero() const {
- return value_ == 0;
- }
-
- extended_exponent_fpt operator-() const {
- return extended_exponent_fpt(-value_, exponent_);
- }
-
- extended_exponent_fpt operator+(const extended_exponent_fpt& that) const {
- if (this->value_ == 0.0 ||
- that.exponent_ > this->exponent_ + fea::kMaxSignificantExpDif) {
- return that;
- }
- if (that.value_ == 0.0 ||
- this->exponent_ > that.exponent_ + fea::kMaxSignificantExpDif) {
- return *this;
- }
- if (this->exponent_ >= that.exponent_) {
- exp_type exp_dif = this->exponent_ - that.exponent_;
- fpt_type value = this->value_;
- fea::set_exponent(value, exp_dif);
- value += that.value_;
- return extended_exponent_fpt(value, that.exponent_);
- } else {
- exp_type exp_dif = that.exponent_ - this->exponent_;
- fpt_type value = that.value_;
- fea::set_exponent(value, exp_dif);
- value += this->value_;
- return extended_exponent_fpt(value, this->exponent_);
- }
- }
-
- extended_exponent_fpt operator-(const extended_exponent_fpt& that) const {
- if (this->value_ == 0.0 ||
- that.exponent_ > this->exponent_ + fea::kMaxSignificantExpDif) {
- return extended_exponent_fpt(-that.value_, that.exponent_);
- }
- if (that.value_ == 0.0 ||
- this->exponent_ > that.exponent_ + fea::kMaxSignificantExpDif) {
- return *this;
- }
- if (this->exponent_ >= that.exponent_) {
- exp_type exp_dif = this->exponent_ - that.exponent_;
- fpt_type value = this->value_;
- fea::set_exponent(value, exp_dif);
- value -= that.value_;
- return extended_exponent_fpt(value, that.exponent_);
- } else {
- exp_type exp_dif = that.exponent_ - this->exponent_;
- fpt_type value = -that.value_;
- fea::set_exponent(value, exp_dif);
- value += this->value_;
- return extended_exponent_fpt(value, this->exponent_);
- }
- }
-
- extended_exponent_fpt operator*(const extended_exponent_fpt& that) const {
- fpt_type value = this->value_ * that.value_;
- exp_type exponent = this->exponent_ + that.exponent_;
- return extended_exponent_fpt(value, exponent);
- }
-
- extended_exponent_fpt operator/(const extended_exponent_fpt& that) const {
- fpt_type value = this->value_ / that.value_;
- exp_type exponent = this->exponent_ - that.exponent_;
- return extended_exponent_fpt(value, exponent);
- }
-
- extended_exponent_fpt& operator+=(const extended_exponent_fpt& that) {
- return *this = *this + that;
- }
-
- extended_exponent_fpt& operator-=(const extended_exponent_fpt& that) {
- return *this = *this - that;
- }
-
- extended_exponent_fpt& operator*=(const extended_exponent_fpt& that) {
- return *this = *this * that;
- }
-
- extended_exponent_fpt& operator/=(const extended_exponent_fpt& that) {
- return *this = *this / that;
- }
-
- extended_exponent_fpt sqrt() const {
- fpt_type val = value_;
- exp_type exp = exponent_;
- if (exp & 1) {
- val *= 2.0;
- --exp;
- }
- return extended_exponent_fpt(get_sqrt(val), exp >> 1);
- }
-
- fpt_type d() const {
- fpt_type ret_val = value_;
- exp_type exp = exponent_;
- if (ret_val == 0.0) {
- return ret_val;
- }
- if (exp >= fea::kMaxExponent) {
- ret_val = 1.0;
- exp = fea::kMaxExponent;
- } else if (exp <= fea::kMinExponent) {
- ret_val = 1.0;
- exp = fea::kMinExponent;
- }
- fea::set_exponent(ret_val, exp);
- return ret_val;
- }
-
- private:
- fpt_type value_;
- exp_type exponent_;
- };
- typedef extended_exponent_fpt<double> efpt64;
-
- template <typename _fpt>
- extended_exponent_fpt<_fpt> get_sqrt(const extended_exponent_fpt<_fpt>& that) {
- return that.sqrt();
- }
-
- template <typename _fpt>
- bool is_pos(const extended_exponent_fpt<_fpt>& that) {
- return that.is_pos();
- }
-
- template <typename _fpt>
- bool is_neg(const extended_exponent_fpt<_fpt>& that) {
- return that.is_neg();
- }
-
- template <typename _fpt>
- bool is_zero(const extended_exponent_fpt<_fpt>& that) {
- return that.is_zero();
- }
-
- template<size_t N>
- class extended_int {
- public:
- static const uint64 kUInt64LowMask;
- static const uint64 kUInt64HighMask;
-
- extended_int() {}
-
- extended_int(int32 that) {
- if (that > 0) {
- this->chunks_[0] = that;
- this->count_ = 1;
- } else if (that < 0) {
- this->chunks_[0] = -that;
- this->count_ = -1;
- } else {
- this->chunks_[0] = this->count_ = 0;
- }
- }
-
- extended_int(int64 that) {
- if (that > 0) {
- this->chunks_[0] = static_cast<uint32>(that & kUInt64LowMask);
- uint32 high = (that & kUInt64HighMask) >> 32;
- if (high) {
- this->chunks_[1] = high;
- this->count_ = 2;
- } else {
- this->count_ = 1;
- }
- } else if (that < 0) {
- that = -that;
- this->chunks_[0] = static_cast<uint32>(that & kUInt64LowMask);
- uint32 high = (that & kUInt64HighMask) >> 32;
- if (high) {
- this->chunks_[1] = high;
- this->count_ = -2;
- } else {
- this->count_ = -1;
- }
- } else {
- this->chunks_[0] = this->count_ = 0;
- }
- }
-
- extended_int(const std::vector<uint32>& chunks, bool plus = true) {
- this->count_ = static_cast<int32>(chunks.size());
- for (size_t i = 0; i < chunks.size(); ++i) {
- this->chunks_[i] = chunks[chunks.size() - i - 1];
- }
- if (!plus) {
- this->count_ = -this->count_;
- }
- }
-
- template <size_t M>
- extended_int(const extended_int<M>& that) {
- this->count_ = that.count();
- memcpy(this->chunks_, that.chunks(), that.size() * sizeof(uint32));
- }
-
- extended_int<N>& operator=(int32 that) {
- if (that > 0) {
- this->chunks_[0] = that;
- this->count_ = 1;
- } else if (that < 0) {
- this->chunks_[0] = -that;
- this->count_ = -1;
- } else {
- this->chunks_[0] = this->count_ = 0;
- }
- return *this;
- }
-
- extended_int<N>& operator=(int64 that) {
- if (that > 0) {
- this->chunks_[0] = static_cast<uint32>(that & kUInt64LowMask);
- uint32 high = (that & kUInt64HighMask) >> 32;
- if (high) {
- this->chunks_[1] = high;
- this->count_ = 2;
- } else {
- this->count_ = 1;
- }
- } else if (that < 0) {
- that = -that;
- this->chunks_[0] = static_cast<uint32>(that & kUInt64LowMask);
- uint32 high = (that & kUInt64HighMask) >> 32;
- if (high) {
- this->chunks_[1] = high;
- this->count_ = -2;
- } else {
- this->count_ = -1;
- }
- } else {
- this->chunks_[0] = this->count_ = 0;
- }
- return *this;
- }
-
- template <size_t M>
- extended_int<N>& operator=(const extended_int<M>& that) {
- this->count_ = that.count();
- size_t sz = (std::min)(N, that.size());
- memcpy(this->chunks_, that.chunks(), sz * sizeof(uint32));
- return *this;
- }
-
- bool is_pos() const {
- return this->count_ > 0;
- }
-
- bool is_neg() const {
- return this->count_ < 0;
- }
-
- bool is_zero() const {
- return this->count_ == 0;
- }
-
- template <size_t M>
- bool operator==(const extended_int<M>& that) const {
- if (this->count_ != that.count()) {
- return false;
- }
- for (size_t i = 0; i < this->size(); ++i) {
- if (this->chunks_[i] != that.chunks()[i]) {
- return false;
- }
- }
- return true;
- }
-
- template <size_t M>
- bool operator!=(const extended_int<M>& that) const {
- return !(*this == that);
- }
-
- template <size_t M>
- bool operator<(const extended_int<M>& that) const {
- if (this->count_ != that.count()) {
- return this->count_ < that.count();
- }
- size_t i = this->size();
- if (!i) {
- return false;
- }
- do {
- --i;
- if (this->chunks_[i] != that.chunks()[i]) {
- return (this->chunks_[i] < that.chunks()[i]) ^ (this->count_ < 0);
- }
- } while (i);
- return false;
- }
-
- template <size_t M>
- bool operator>(const extended_int<M>& that) const {
- return that < *this;
- }
-
- template <size_t M>
- bool operator<=(const extended_int<M>& that) const {
- return !(that < *this);
- }
-
- template <size_t M>
- bool operator>=(const extended_int<M>& that) const {
- return !(*this < that);
- }
-
- extended_int<N> operator-() const {
- extended_int<N> ret_val = *this;
- ret_val.neg();
- return ret_val;
- }
-
- void neg() {
- this->count_ = -this->count_;
- }
-
- template <size_t M>
- extended_int<(N>M?N:M)> operator+(const extended_int<M>& that) const {
- extended_int<(N>M?N:M)> ret_val;
- ret_val.add(*this, that);
- return ret_val;
- }
-
- template <size_t N1, size_t N2>
- void add(const extended_int<N1>& e1, const extended_int<N2>& e2) {
- if (!e1.count()) {
- *this = e2;
- return;
- }
- if (!e2.count()) {
- *this = e1;
- return;
- }
- if ((e1.count() > 0) ^ (e2.count() > 0)) {
- dif(e1.chunks(), e1.size(), e2.chunks(), e2.size());
- } else {
- add(e1.chunks(), e1.size(), e2.chunks(), e2.size());
- }
- if (e1.count() < 0) {
- this->count_ = -this->count_;
- }
- }
-
- template <size_t M>
- extended_int<(N>M?N:M)> operator-(const extended_int<M>& that) const {
- extended_int<(N>M?N:M)> ret_val;
- ret_val.dif(*this, that);
- return ret_val;
- }
-
- template <size_t N1, size_t N2>
- void dif(const extended_int<N1>& e1, const extended_int<N2> &e2) {
- if (!e1.count()) {
- *this = e2;
- this->count_ = -this->count_;
- return;
- }
- if (!e2.count()) {
- *this = e1;
- return;
- }
- if ((e1.count() > 0) ^ (e2.count() > 0)) {
- add(e1.chunks(), e1.size(), e2.chunks(), e2.size());
- } else {
- dif(e1.chunks(), e1.size(), e2.chunks(), e2.size());
- }
- if (e1.count() < 0) {
- this->count_ = -this->count_;
- }
- }
-
- extended_int<N> operator*(int32 that) const {
- extended_int<N> temp(that);
- return (*this) * temp;
- }
-
- extended_int<N> operator*(int64 that) const {
- extended_int<N> temp(that);
- return (*this) * temp;
- }
-
- template <size_t M>
- extended_int<(N>M?N:M)> operator*(const extended_int<M>& that) const {
- extended_int<(N>M?N:M)> ret_val;
- ret_val.mul(*this, that);
- return ret_val;
- }
-
- template <size_t N1, size_t N2>
- void mul(const extended_int<N1>& e1, const extended_int<N2>& e2) {
- if (!e1.count() || !e2.count()) {
- this->chunks_[0] = this->count_ = 0;
- return;
- }
- mul(e1.chunks(), e1.size(), e2.chunks(), e2.size());
- if ((e1.count() > 0) ^ (e2.count() > 0)) {
- this->count_ = -this->count_;
- }
- }
-
- const uint32* chunks() const {
- return chunks_;
- }
-
- int32 count() const {
- return count_;
- }
-
- size_t size() const {
- return (std::abs)(count_);
- }
-
- std::pair<fpt64, int64> p() const {
- std::pair<fpt64, int64> ret_val(0, 0);
- size_t sz = this->size();
- if (!sz) {
- return ret_val;
- } else {
- if (sz == 1) {
- ret_val.first = static_cast<fpt64>(this->chunks_[0]);
- ret_val.second = 0;
- } else if (sz == 2) {
- ret_val.first = static_cast<fpt64>(this->chunks_[1]) *
- static_cast<fpt64>(0x100000000LL) +
- static_cast<fpt64>(this->chunks_[0]);
- ret_val.second = 0;
- } else {
- for (size_t i = 1; i <= 3; ++i) {
- ret_val.first *= static_cast<fpt64>(0x100000000LL);
- ret_val.first += static_cast<fpt64>(this->chunks_[sz - i]);
- }
- ret_val.second = (sz - 3) << 5;
- }
- }
- if (this->count_ < 0) {
- ret_val.first = -ret_val.first;
- }
- return ret_val;
- }
-
- fpt64 d() const {
- std::pair<fpt64, int64> p = this->p();
- extended_exponent_fpt<fpt64> efpt(p.first, p.second);
- return efpt.d();
- }
-
- private:
- void add(const uint32* c1, size_t sz1, const uint32* c2, size_t sz2) {
- if (sz1 < sz2) {
- add(c2, sz2, c1, sz1);
- return;
- }
- this->count_ = sz1;
- uint64 temp = 0;
- for (size_t i = 0; i < sz2; ++i) {
- temp += static_cast<uint64>(c1[i]) +
- static_cast<uint64>(c2[i]);
- this->chunks_[i] = static_cast<uint32>(temp & kUInt64LowMask);
- temp = (temp & kUInt64HighMask) >> 32;
- }
- for (size_t i = sz2; i < sz1; ++i) {
- temp += static_cast<uint64>(c1[i]);
- this->chunks_[i] = static_cast<uint32>(temp & kUInt64LowMask);
- temp = (temp & kUInt64HighMask) >> 32;
- }
- if (temp && (this->count_ != N)) {
- this->chunks_[this->count_] = static_cast<uint32>(temp & kUInt64LowMask);
- ++this->count_;
- }
- }
-
- void dif(const uint32* c1, size_t sz1, const uint32* c2, size_t sz2, bool rec = false) {
- if (sz1 < sz2) {
- dif(c2, sz2, c1, sz1, true);
- this->count_ = -this->count_;
- return;
- } else if ((sz1 == sz2) && !rec) {
- do {
- --sz1;
- if (c1[sz1] < c2[sz1]) {
- ++sz1;
- dif(c2, sz1, c1, sz1, true);
- this->count_ = -this->count_;
- return;
- } else if (c1[sz1] > c2[sz1]) {
- ++sz1;
- break;
- } else if (!sz1) {
- this->chunks_[0] = this->count_ = 0;
- return;
- }
- } while (sz1);
- sz2 = sz1;
- }
- this->count_ = sz1-1;
- bool flag = false;
- for (size_t i = 0; i < sz2; ++i) {
- this->chunks_[i] = c1[i] - c2[i] - (flag?1:0);
- flag = (c1[i] < c2[i]) || ((c1[i] == c2[i]) && flag);
- }
- for (size_t i = sz2; i < sz1; ++i) {
- this->chunks_[i] = c1[i] - (flag?1:0);
- flag = !c1[i] && flag;
- }
- if (this->chunks_[this->count_]) {
- ++this->count_;
- }
- }
-
- void mul(const uint32* c1, size_t sz1, const uint32* c2, size_t sz2) {
- uint64 cur = 0, nxt, tmp;
- this->count_ = static_cast<int32>((std::min)(N, sz1 + sz2 - 1));
- for (size_t shift = 0; shift < static_cast<size_t>(this->count_); ++shift) {
- nxt = 0;
- for (size_t first = 0; first <= shift; ++first) {
- if (first >= sz1) {
- break;
- }
- size_t second = shift - first;
- if (second >= sz2) {
- continue;
- }
- tmp = static_cast<uint64>(c1[first]) *
- static_cast<uint64>(c2[second]);
- cur += tmp & kUInt64LowMask;
- nxt += (tmp & kUInt64HighMask) >> 32;
- }
- this->chunks_[shift] = static_cast<uint32>(cur & kUInt64LowMask);
- cur = nxt + ((cur & kUInt64HighMask) >> 32);
- }
- if (cur && (this->count_ != N)) {
- this->chunks_[this->count_] = static_cast<uint32>(cur & kUInt64LowMask);
- ++this->count_;
- }
- }
-
- uint32 chunks_[N];
- int32 count_;
- };
-
- template <size_t N>
- const uint64 extended_int<N>::kUInt64LowMask = 0x00000000ffffffffULL;
- template <size_t N>
- const uint64 extended_int<N>::kUInt64HighMask = 0xffffffff00000000ULL;
-
- typedef extended_int<1> eint32;
- typedef extended_int<2> eint64;
- typedef extended_int<4> eint128;
- typedef extended_int<8> eint256;
- typedef extended_int<16> eint512;
- typedef extended_int<32> eint1024;
- typedef extended_int<64> eint2048;
- typedef extended_int<128> eint4096;
- typedef extended_int<256> eint8192;
- typedef extended_int<512> eint16384;
- typedef extended_int<1024> eint32768;
-
- template <size_t N>
- bool is_pos(const extended_int<N>& that) {
- return that.count() > 0;
- }
-
- template <size_t N>
- bool is_neg(const extended_int<N>& that) {
- return that.count() < 0;
- }
-
- template <size_t N>
- bool is_zero(const extended_int<N>& that) {
- return !that.count();
- }
-
- struct type_converter_fpt {
- template <typename T>
- fpt64 operator()(const T& that) const {
- return static_cast<fpt64>(that);
- }
-
- template <size_t N>
- fpt64 operator()(const extended_int<N>& that) const {
- return that.d();
- }
-
- template <>
- fpt64 operator()< extended_exponent_fpt<fpt64> >(const extended_exponent_fpt<fpt64>& that) const {
- return that.d();
- }
- };
-
- struct type_converter_efpt {
- template <size_t N>
- extended_exponent_fpt<fpt64> operator()(const extended_int<N>& that) const {
- std::pair<fpt64, int64> p = that.p();
- return extended_exponent_fpt<fpt64>(p.first, p.second);
- }
- };
-
- // Used to compute expressions that operate with sqrts with predefined
- // relative error. Evaluates expressions of the next type:
- // sum(i = 1 .. n)(A[i] * sqrt(B[i])), 1 <= n <= 4.
- template <typename _int, typename _fpt, typename _converter>
- class robust_sqrt_expr {
- public:
- static const unsigned int EVAL1_MAX_RELATIVE_ERROR;
- static const unsigned int EVAL2_MAX_RELATIVE_ERROR;
- static const unsigned int EVAL3_MAX_RELATIVE_ERROR;
- static const unsigned int EVAL4_MAX_RELATIVE_ERROR;
-
- // Evaluates expression (re = 4 EPS):
- // A[0] * sqrt(B[0]).
- _fpt eval1(_int *A, _int *B) {
- _fpt a = convert(A[0]);
- _fpt b = convert(B[0]);
- return a * get_sqrt(b);
- }
-
- // Evaluates expression (re = 7 EPS):
- // A[0] * sqrt(B[0]) + A[1] * sqrt(B[1]).
- _fpt eval2(_int *A, _int *B) {
- _fpt a = eval1(A, B);
- _fpt b = eval1(A + 1, B + 1);
- if ((!is_neg(a) && !is_neg(b)) ||
- (!is_pos(a) && !is_pos(b)))
- return a + b;
- return convert(A[0] * A[0] * B[0] - A[1] * A[1] * B[1]) / (a - b);
- }
-
- // Evaluates expression (re = 16 EPS):
- // A[0] * sqrt(B[0]) + A[1] * sqrt(B[1]) + A[2] * sqrt(B[2]).
- _fpt eval3(_int *A, _int *B) {
- _fpt a = eval2(A, B);
- _fpt b = eval1(A + 2, B + 2);
- if ((!is_neg(a) && !is_neg(b)) ||
- (!is_pos(a) && !is_pos(b)))
- return a + b;
- tA[3] = A[0] * A[0] * B[0] + A[1] * A[1] * B[1] - A[2] * A[2] * B[2];
- tB[3] = 1;
- tA[4] = A[0] * A[1] * 2;
- tB[4] = B[0] * B[1];
- return eval2(tA + 3, tB + 3) / (a - b);
- }
-
-
- // Evaluates expression (re = 25 EPS):
- // A[0] * sqrt(B[0]) + A[1] * sqrt(B[1]) +
- // A[2] * sqrt(B[2]) + A[3] * sqrt(B[3]).
- _fpt eval4(_int *A, _int *B) {
- _fpt a = eval2(A, B);
- _fpt b = eval2(A + 2, B + 2);
- if ((!is_neg(a) && !is_neg(b)) ||
- (!is_pos(a) && !is_pos(b)))
- return a + b;
- tA[0] = A[0] * A[0] * B[0] + A[1] * A[1] * B[1] -
- A[2] * A[2] * B[2] - A[3] * A[3] * B[3];
- tB[0] = 1;
- tA[1] = A[0] * A[1] * 2;
- tB[1] = B[0] * B[1];
- tA[2] = A[2] * A[3] * -2;
- tB[2] = B[2] * B[3];
- return eval3(tA, tB) / (a - b);
- }
-
- private:
- _int tA[5];
- _int tB[5];
- _converter convert;
- };
-
- template <typename _int, typename _fpt, typename _converter>
- const unsigned int robust_sqrt_expr<_int, _fpt, _converter>::EVAL1_MAX_RELATIVE_ERROR = 4;
- template <typename _int, typename _fpt, typename _converter>
- const unsigned int robust_sqrt_expr<_int, _fpt, _converter>::EVAL2_MAX_RELATIVE_ERROR = 7;
- template <typename _int, typename _fpt, typename _converter>
- const unsigned int robust_sqrt_expr<_int, _fpt, _converter>::EVAL3_MAX_RELATIVE_ERROR = 16;
- template <typename _int, typename _fpt, typename _converter>
- const unsigned int robust_sqrt_expr<_int, _fpt, _converter>::EVAL4_MAX_RELATIVE_ERROR = 25;
+ positive_sum_ /= -val;
+ negative_sum_ /= -val;
+ swap();
+ }
+ return *this;
+ }
+
+private:
+ void swap() {
+ (std::swap)(positive_sum_, negative_sum_);
+ }
+
+ T positive_sum_;
+ T negative_sum_;
+};
+
+template<typename T>
+robust_dif<T> operator+(const robust_dif<T>& lhs,
+ const robust_dif<T>& rhs) {
+ return robust_dif<T>(lhs.pos() + rhs.pos(),
+ lhs.neg() + rhs.neg());
+}
+
+template<typename T>
+robust_dif<T> operator+(const robust_dif<T>& lhs, const T& rhs) {
+ if (!is_neg(rhs)) {
+ return robust_dif<T>(lhs.pos() + rhs, lhs.neg());
+ } else {
+ return robust_dif<T>(lhs.pos(), lhs.neg() - rhs);
+ }
+}
+
+template<typename T>
+robust_dif<T> operator+(const T& lhs, const robust_dif<T>& rhs) {
+ if (!is_neg(lhs)) {
+ return robust_dif<T>(lhs + rhs.pos(), rhs.neg());
+ } else {
+ return robust_dif<T>(rhs.pos(), rhs.neg() - lhs);
+ }
+}
+
+template<typename T>
+robust_dif<T> operator-(const robust_dif<T>& lhs,
+ const robust_dif<T>& rhs) {
+ return robust_dif<T>(lhs.pos() + rhs.neg(), lhs.neg() + rhs.pos());
+}
+
+template<typename T>
+robust_dif<T> operator-(const robust_dif<T>& lhs, const T& rhs) {
+ if (!is_neg(rhs)) {
+ return robust_dif<T>(lhs.pos(), lhs.neg() + rhs);
+ } else {
+ return robust_dif<T>(lhs.pos() - rhs, lhs.neg());
+ }
+}
+
+template<typename T>
+robust_dif<T> operator-(const T& lhs, const robust_dif<T>& rhs) {
+ if (!is_neg(lhs)) {
+ return robust_dif<T>(lhs + rhs.neg(), rhs.pos());
+ } else {
+ return robust_dif<T>(rhs.neg(), rhs.pos() - lhs);
+ }
+}
+
+template<typename T>
+robust_dif<T> operator*(const robust_dif<T>& lhs,
+ const robust_dif<T>& rhs) {
+ T res_pos = lhs.pos() * rhs.pos() + lhs.neg() * rhs.neg();
+ T res_neg = lhs.pos() * rhs.neg() + lhs.neg() * rhs.pos();
+ return robust_dif<T>(res_pos, res_neg);
+}
+
+template<typename T>
+robust_dif<T> operator*(const robust_dif<T>& lhs, const T& val) {
+ if (!is_neg(val)) {
+ return robust_dif<T>(lhs.pos() * val, lhs.neg() * val);
+ } else {
+ return robust_dif<T>(-lhs.neg() * val, -lhs.pos() * val);
+ }
+}
+
+template<typename T>
+robust_dif<T> operator*(const T& val, const robust_dif<T>& rhs) {
+ if (!is_neg(val)) {
+ return robust_dif<T>(val * rhs.pos(), val * rhs.neg());
+ } else {
+ return robust_dif<T>(-val * rhs.neg(), -val * rhs.pos());
+ }
+}
+
+template<typename T>
+robust_dif<T> operator/(const robust_dif<T>& lhs, const T& val) {
+ if (!is_neg(val)) {
+ return robust_dif<T>(lhs.pos() / val, lhs.neg() / val);
+ } else {
+ return robust_dif<T>(-lhs.neg() / val, -lhs.pos() / val);
+ }
+}
+
+// Used to compute expressions that operate with sqrts with predefined
+// relative error. Evaluates expressions of the next type:
+// sum(i = 1 .. n)(A[i] * sqrt(B[i])), 1 <= n <= 4.
+template <typename _int, typename _fpt, typename _converter>
+class robust_sqrt_expr {
+public:
+ static const unsigned int EVAL1_MAX_RELATIVE_ERROR;
+ static const unsigned int EVAL2_MAX_RELATIVE_ERROR;
+ static const unsigned int EVAL3_MAX_RELATIVE_ERROR;
+ static const unsigned int EVAL4_MAX_RELATIVE_ERROR;
+
+ // Evaluates expression (re = 4 EPS):
+ // A[0] * sqrt(B[0]).
+ _fpt eval1(_int *A, _int *B) {
+ _fpt a = convert(A[0]);
+ _fpt b = convert(B[0]);
+ return a * get_sqrt(b);
+ }
+
+ // Evaluates expression (re = 7 EPS):
+ // A[0] * sqrt(B[0]) + A[1] * sqrt(B[1]).
+ _fpt eval2(_int *A, _int *B) {
+ _fpt a = eval1(A, B);
+ _fpt b = eval1(A + 1, B + 1);
+ if ((!is_neg(a) && !is_neg(b)) ||
+ (!is_pos(a) && !is_pos(b)))
+ return a + b;
+ return convert(A[0] * A[0] * B[0] - A[1] * A[1] * B[1]) / (a - b);
+ }
+
+ // Evaluates expression (re = 16 EPS):
+ // A[0] * sqrt(B[0]) + A[1] * sqrt(B[1]) + A[2] * sqrt(B[2]).
+ _fpt eval3(_int *A, _int *B) {
+ _fpt a = eval2(A, B);
+ _fpt b = eval1(A + 2, B + 2);
+ if ((!is_neg(a) && !is_neg(b)) ||
+ (!is_pos(a) && !is_pos(b)))
+ return a + b;
+ tA[3] = A[0] * A[0] * B[0] + A[1] * A[1] * B[1] - A[2] * A[2] * B[2];
+ tB[3] = 1;
+ tA[4] = A[0] * A[1] * 2;
+ tB[4] = B[0] * B[1];
+ return eval2(tA + 3, tB + 3) / (a - b);
+ }
+
+
+ // Evaluates expression (re = 25 EPS):
+ // A[0] * sqrt(B[0]) + A[1] * sqrt(B[1]) +
+ // A[2] * sqrt(B[2]) + A[3] * sqrt(B[3]).
+ _fpt eval4(_int *A, _int *B) {
+ _fpt a = eval2(A, B);
+ _fpt b = eval2(A + 2, B + 2);
+ if ((!is_neg(a) && !is_neg(b)) ||
+ (!is_pos(a) && !is_pos(b)))
+ return a + b;
+ tA[0] = A[0] * A[0] * B[0] + A[1] * A[1] * B[1] -
+ A[2] * A[2] * B[2] - A[3] * A[3] * B[3];
+ tB[0] = 1;
+ tA[1] = A[0] * A[1] * 2;
+ tB[1] = B[0] * B[1];
+ tA[2] = A[2] * A[3] * -2;
+ tB[2] = B[2] * B[3];
+ return eval3(tA, tB) / (a - b);
+ }
+
+private:
+ _int tA[5];
+ _int tB[5];
+ _converter convert;
+};
+
+template <typename _int, typename _fpt, typename _converter>
+const unsigned int robust_sqrt_expr<_int, _fpt, _converter>::EVAL1_MAX_RELATIVE_ERROR = 4;
+template <typename _int, typename _fpt, typename _converter>
+const unsigned int robust_sqrt_expr<_int, _fpt, _converter>::EVAL2_MAX_RELATIVE_ERROR = 7;
+template <typename _int, typename _fpt, typename _converter>
+const unsigned int robust_sqrt_expr<_int, _fpt, _converter>::EVAL3_MAX_RELATIVE_ERROR = 16;
+template <typename _int, typename _fpt, typename _converter>
+const unsigned int robust_sqrt_expr<_int, _fpt, _converter>::EVAL4_MAX_RELATIVE_ERROR = 25;
} // detail
} // polygon
} // boost
Modified: sandbox/gtl/boost/polygon/detail/voronoi_structures.hpp
==============================================================================
--- sandbox/gtl/boost/polygon/detail/voronoi_structures.hpp (original)
+++ sandbox/gtl/boost/polygon/detail/voronoi_structures.hpp 2012-01-23 17:22:43 EST (Mon, 23 Jan 2012)
@@ -1,14 +1,14 @@
// Boost.Polygon library detail/voronoi_structures.hpp header file
-// Copyright Andrii Sydorchuk 2010-2011.
+// Copyright Andrii Sydorchuk 2010-2012.
// Distributed under the Boost Software License, Version 1.0.
// (See accompanying file LICENSE_1_0.txt or copy at
// http://www.boost.org/LICENSE_1_0.txt)
// See http://www.boost.org for updates, documentation, and revision history.
-#ifndef BOOST_POLYGON_VORONOI_STRUCTURES
-#define BOOST_POLYGON_VORONOI_STRUCTURES
+#ifndef BOOST_POLYGON_VORONOI_DETAIL_STRUCTURES
+#define BOOST_POLYGON_VORONOI_DETAIL_STRUCTURES
#include <list>
#include <queue>
Modified: sandbox/gtl/boost/polygon/voronoi.hpp
==============================================================================
--- sandbox/gtl/boost/polygon/voronoi.hpp (original)
+++ sandbox/gtl/boost/polygon/voronoi.hpp 2012-01-23 17:22:43 EST (Mon, 23 Jan 2012)
@@ -1,6 +1,6 @@
// Boost.Polygon library voronoi.hpp header file
-// Copyright Andrii Sydorchuk 2010-2011.
+// Copyright Andrii Sydorchuk 2010-2012.
// Distributed under the Boost Software License, Version 1.0.
// (See accompanying file LICENSE_1_0.txt or copy at
// http://www.boost.org/LICENSE_1_0.txt)
@@ -10,7 +10,6 @@
#ifndef BOOST_POLYGON_VORONOI
#define BOOST_POLYGON_VORONOI
-#include "polygon.hpp"
#include "voronoi_builder.hpp"
#include "voronoi_diagram.hpp"
@@ -26,28 +25,28 @@
// the signed integer range [-2^31, 2^31-1].
// Complexity - O(N*logN), memory usage - O(N), where N is the total
// number of points and segments.
- template <typename T, typename PC, typename VD>
+ template <typename PC, typename VD>
static inline void construct_voronoi_points(
- const PC &points, VD &output) {
- voronoi_builder<T, VD> builder;
+ const PC &points, VD *output) {
+ default_voronoi_builder builder;
builder.insert_points(points.begin(), points.end());
builder.construct(output);
builder.clear();
}
- template <typename T, typename SC, typename VD>
+ template <typename SC, typename VD>
static inline void construct_voronoi_segments(
- const SC &segments, VD &output) {
- voronoi_builder<T, VD> builder;
+ const SC &segments, VD *output) {
+ default_voronoi_builder builder;
builder.insert_segments(segments.begin(), segments.end());
builder.construct(output);
builder.clear();
}
- template <typename T, typename PC, typename SC, typename VD>
+ template <typename PC, typename SC, typename VD>
static inline void construct_voronoi(
- const PC &points, const SC &segments, VD &output) {
- voronoi_builder<T, VD> builder;
+ const PC &points, const SC &segments, VD *output) {
+ default_voronoi_builder builder;
builder.insert_sites(points.begin(), points.end(),
segments.begin(), segments.end());
builder.construct(output);
Modified: sandbox/gtl/boost/polygon/voronoi_builder.hpp
==============================================================================
--- sandbox/gtl/boost/polygon/voronoi_builder.hpp (original)
+++ sandbox/gtl/boost/polygon/voronoi_builder.hpp 2012-01-23 17:22:43 EST (Mon, 23 Jan 2012)
@@ -1,6 +1,6 @@
// Boost.Polygon library voronoi_builder.hpp header file
-// Copyright Andrii Sydorchuk 2010-2011.
+// Copyright Andrii Sydorchuk 2010-2012.
// Distributed under the Boost Software License, Version 1.0.
// (See accompanying file LICENSE_1_0.txt or copy at
// http://www.boost.org/LICENSE_1_0.txt)
@@ -14,8 +14,7 @@
#include <map>
#include <vector>
-#include "polygon.hpp"
-
+#include "detail/voronoi_ctypes.hpp"
#include "detail/voronoi_predicates.hpp"
#include "detail/voronoi_structures.hpp"
@@ -61,11 +60,11 @@
// defined comparison functor is used to make the beach line correctly
// ordered. Epsilon-based and high-precision approaches are used to guarantee
// efficiency and robustness of the algorithm implementation.
- template <typename T, typename VD>
+ template <typename T,
+ typename CTT = detail::voronoi_ctype_traits<T>,
+ typename VP = detail::voronoi_predicates<CTT> >
class voronoi_builder {
public:
- typedef VD output_type;
-
voronoi_builder() {}
template <typename PointIterator>
@@ -108,13 +107,12 @@
}
// Run the sweepline algorithm.
- void construct(output_type &output) {
+ template <typename OUTPUT>
+ void construct(OUTPUT *output) {
// Init structures.
- output_ = &output;
- output_->clear();
- output_->reserve(site_events_.size());
+ output->reserve(site_events_.size());
init_sites_queue();
- init_beach_line();
+ init_beach_line(output);
// The algorithm stops when there are no events to process.
// The circle events with the same rightmost point as the next
@@ -123,14 +121,14 @@
while (!circle_events_.empty() ||
!(site_event_iterator_ == site_events_.end())) {
if (circle_events_.empty()) {
- process_site_event();
+ process_site_event(output);
} else if (site_event_iterator_ == site_events_.end()) {
- process_circle_event();
+ process_circle_event(output);
} else {
if (event_comparison(*site_event_iterator_, circle_events_.top().first)) {
- process_site_event();
+ process_site_event(output);
} else {
- process_circle_event();
+ process_circle_event(output);
}
}
while (!circle_events_.empty() && !circle_events_.top().first.is_active()) {
@@ -141,7 +139,7 @@
beach_line_.clear();
// Clean the output (remove zero-length edges).
- output_->clean();
+ output->clean();
}
void clear() {
@@ -149,9 +147,8 @@
}
private:
- typedef detail::voronoi_predicates<T> VP;
- typedef typename VP::int_type int_type;
- typedef typename VP::fpt_type fpt_type;
+ typedef typename CTT::int_type int_type;
+ typedef typename CTT::fpt_type fpt_type;
typedef detail::point_2d<int_type> point_type;
typedef detail::site_event<int_type> site_event_type;
@@ -166,7 +163,7 @@
typedef typename VP::
template circle_formation_predicate<site_event_type, circle_event_type>
circle_formation_predicate_type;
- typedef typename output_type::voronoi_edge_type edge_type;
+ typedef void edge_type;
typedef detail::beach_line_node_key<site_event_type> key_type;
typedef detail::beach_line_node_data<edge_type, circle_event_type> value_type;
typedef typename VP::
@@ -204,11 +201,12 @@
site_event_iterator_ = site_events_.begin();
}
- void init_beach_line() {
+ template <typename OUTPUT>
+ void init_beach_line(OUTPUT *output) {
if (site_events_.empty()) return;
if (site_events_.size() == 1) {
// Handle one input site case.
- output_->process_single_site(site_events_[0]);
+ output->process_single_site(site_events_[0]);
++site_event_iterator_;
} else {
int skip = 0;
@@ -224,26 +222,28 @@
if (skip == 1) {
// Init the beach line with the two first sites.
- init_beach_line_default();
+ init_beach_line_default(output);
} else {
// Init the beach line with the sites situated on the same
// vertical line. This could be a set of point and vertical
// segment sites.
- init_beach_line_collinear_sites();
+ init_beach_line_collinear_sites(output);
}
}
}
// Init the beach line with the two first sites.
// The first site is always a point.
- void init_beach_line_default() {
+ template <typename OUTPUT>
+ void init_beach_line_default(OUTPUT *output) {
// Get the first and the second site events.
site_event_iterator_type it_first = site_events_.begin();
site_event_iterator_type it_second = site_events_.begin();
++it_second;
// Update the beach line.
- insert_new_arc(*it_first, *it_first, *it_second, beach_line_.begin());
+ insert_new_arc(*it_first, *it_first, *it_second,
+ beach_line_.begin(), output);
// The second site has been already processed.
// Move the site's iterator.
@@ -251,7 +251,8 @@
}
// Insert bisectors for all collinear initial sites.
- void init_beach_line_collinear_sites() {
+ template <typename OUTPUT>
+ void init_beach_line_collinear_sites(OUTPUT *output) {
site_event_iterator_type it_first = site_events_.begin();
site_event_iterator_type it_second = site_events_.begin();
++it_second;
@@ -260,7 +261,7 @@
key_type new_node(*it_first, *it_second);
// Update the output.
- edge_type *edge = output_->insert_new_edge(*it_first, *it_second);
+ edge_type *edge = output->insert_new_edge(*it_first, *it_second).first;
// Insert a new bisector into the beach line.
beach_line_.insert(
@@ -279,7 +280,8 @@
}
}
- void process_site_event() {
+ template <typename OUTPUT>
+ void process_site_event(OUTPUT *output) {
// Get the site event to process.
site_event_type site_event = *site_event_iterator_;
@@ -324,7 +326,7 @@
// Insert new nodes into the beach line. Update the output.
beach_line_iterator new_node_it =
- insert_new_arc(site_arc, site_arc, site_event, it);
+ insert_new_arc(site_arc, site_arc, site_event, it, output);
// Add a candidate circle to the circle event queue.
// There could be only one new circle event formed by
@@ -338,7 +340,7 @@
const site_event_type &site_arc = it->first.left_site();
// Insert new nodes into the beach line. Update the output.
- insert_new_arc(site_arc, site_arc, site_event, it);
+ insert_new_arc(site_arc, site_arc, site_event, it, output);
// If the site event is a segment, update its direction.
if (site_event.is_segment()) {
@@ -364,7 +366,7 @@
// Insert new nodes into the beach line. Update the output.
beach_line_iterator new_node_it =
- insert_new_arc(site_arc1, site_arc2, site_event, it);
+ insert_new_arc(site_arc1, site_arc2, site_event, it, output);
// Add candidate circles to the circle event queue.
// There could be up to two circle events formed by
@@ -393,7 +395,8 @@
// (B, C) bisector and change (A, B) bisector to the (A, C). That's
// why we use const_cast there and take all the responsibility that
// map data structure keeps correct ordering.
- void process_circle_event() {
+ template <typename OUTPUT>
+ void process_circle_event(OUTPUT *output) {
// Get the topmost circle event.
const event_type &e = circle_events_.top();
const circle_event_type &circle_event = e.first;
@@ -422,8 +425,8 @@
const_cast<key_type &>(it_first->first).right_site(site3);
// Insert the new bisector into the beach line.
- it_first->second.edge(output_->insert_new_edge(site1, site3, circle_event,
- bisector1, bisector2));
+ it_first->second.edge(output->insert_new_edge(site1, site3, circle_event,
+ bisector1, bisector2));
// Remove the (B, C) bisector node from the beach line.
beach_line_.erase(it_last);
@@ -452,10 +455,12 @@
}
// Insert new nodes into the beach line. Update the output.
+ template <typename OUTPUT>
beach_line_iterator insert_new_arc(const site_event_type &site_arc1,
const site_event_type &site_arc2,
const site_event_type &site_event,
- const beach_line_iterator &position) {
+ const beach_line_iterator &position,
+ OUTPUT *output) {
// Create two new bisectors with opposite directions.
key_type new_left_node(site_arc1, site_event);
key_type new_right_node(site_event, site_arc2);
@@ -466,9 +471,10 @@
}
// Update the output.
- edge_type *edge = output_->insert_new_edge(site_arc2, site_event);
+ std::pair<edge_type*, edge_type*> edges =
+ output->insert_new_edge(site_arc2, site_event);
beach_line_iterator it = beach_line_.insert(position,
- typename beach_line_type::value_type(new_right_node, value_type(edge->twin())));
+ beach_line_type::value_type(new_right_node, value_type(edges.second)));
if (site_event.is_segment()) {
// Update the beach line with temporary bisector, that will
@@ -484,7 +490,7 @@
}
beach_line_.insert(position,
- typename beach_line_type::value_type(new_left_node, value_type(edge)));
+ typename beach_line_type::value_type(new_left_node, value_type(edges.first)));
return it;
}
@@ -520,13 +526,14 @@
end_point_comparison > end_points_;
circle_event_queue_type circle_events_;
beach_line_type beach_line_;
- output_type *output_;
circle_formation_predicate_type circle_formation_predicate_;
//Disallow copy constructor and operator=
voronoi_builder(const voronoi_builder&);
void operator=(const voronoi_builder&);
};
+
+ typedef voronoi_builder<detail::int32> default_voronoi_builder;
} // polygon
} // boost
Modified: sandbox/gtl/boost/polygon/voronoi_diagram.hpp
==============================================================================
--- sandbox/gtl/boost/polygon/voronoi_diagram.hpp (original)
+++ sandbox/gtl/boost/polygon/voronoi_diagram.hpp 2012-01-23 17:22:43 EST (Mon, 23 Jan 2012)
@@ -1,6 +1,6 @@
// Boost.Polygon library voronoi_diagram.hpp header file
-// Copyright Andrii Sydorchuk 2010-2011.
+// Copyright Andrii Sydorchuk 2010-2012.
// Distributed under the Boost Software License, Version 1.0.
// (See accompanying file LICENSE_1_0.txt or copy at
// http://www.boost.org/LICENSE_1_0.txt)
@@ -846,8 +846,8 @@
// Takes as input left and right sites that form a new bisector.
// Returns a pointer to a new half-edge.
template <typename SEvent>
- voronoi_edge_type *insert_new_edge(const SEvent &site1,
- const SEvent &site2) {
+ std::pair<void*, void*> insert_new_edge(const SEvent &site1,
+ const SEvent &site2) {
// Get sites' indices.
int site_index1 = site1.index();
int site_index2 = site2.index();
@@ -892,7 +892,7 @@
edge2.twin(&edge1);
// Return a pointer to the new half-edge.
- return &edge1;
+ return std::make_pair(&edge1, &edge2);
}
// Insert a new half-edge into the output data structure with the
@@ -902,11 +902,14 @@
// pointers to those half-edges. Half-edges' direction goes out of the
// new voronoi vertex point. Returns a pointer to the new half-edge.
template <typename SEvent, typename CEvent>
- voronoi_edge_type *insert_new_edge(const SEvent &site1,
- const SEvent &site3,
- const CEvent &circle,
- voronoi_edge_type *edge12,
- voronoi_edge_type *edge23) {
+ void *insert_new_edge(const SEvent &site1,
+ const SEvent &site3,
+ const CEvent &circle,
+ void *data12,
+ void *data23) {
+ voronoi_edge_type *edge12 = static_cast<voronoi_edge_type*>(data12);
+ voronoi_edge_type *edge23 = static_cast<voronoi_edge_type*>(data23);
+
// Add a new voronoi vertex.
vertex_records_.push_back(voronoi_vertex_type(
point_type(circle.x(), circle.y()), edge12));
Modified: sandbox/gtl/libs/polygon/test/voronoi_benchmark_test.cpp
==============================================================================
--- sandbox/gtl/libs/polygon/test/voronoi_benchmark_test.cpp (original)
+++ sandbox/gtl/libs/polygon/test/voronoi_benchmark_test.cpp 2012-01-23 17:22:43 EST (Mon, 23 Jan 2012)
@@ -1,6 +1,6 @@
// Boost.Polygon library voronoi_benchmark_test.cpp file
-// Copyright Andrii Sydorchuk 2010-2011.
+// Copyright Andrii Sydorchuk 2010-2012.
// Distributed under the Boost Software License, Version 1.0.
// (See accompanying file LICENSE_1_0.txt or copy at
// http://www.boost.org/LICENSE_1_0.txt)
@@ -55,12 +55,13 @@
timer.restart();
int num_tests = max_points / num_points;
for (int cur = 0; cur < num_tests; cur++) {
+ test_output.clear();
for (int cur_point = 0; cur_point < num_points; cur_point++) {
x = static_cast<coordinate_type>(gen());
y = static_cast<coordinate_type>(gen());
points[cur_point] = point_type(x, y);
}
- construct_voronoi_points<coordinate_type>(points, test_output);
+ construct_voronoi_points(points, &test_output);
}
double elapsed_time = timer.elapsed();
double time_per_test = elapsed_time / num_tests;
@@ -97,7 +98,7 @@
for (int i = 0; i < POINT_RUNS; ++i) {
voronoi_diagram<double> test_output;
timer.restart();
- construct_voronoi_points<coordinate_type>(points, test_output);
+ construct_voronoi_points(points, &test_output);
periods.push_back(timer.elapsed());
}
}
@@ -137,7 +138,7 @@
for (int i = 0; i < SEGMENT_RUNS; ++i) {
voronoi_diagram<double> test_output;
timer.restart();
- construct_voronoi_segments<coordinate_type>(segments, test_output);
+ construct_voronoi_segments(segments, &test_output);
periods.push_back(timer.elapsed());
}
}
Modified: sandbox/gtl/libs/polygon/voronoi_example/voronoi_visualizer.cpp
==============================================================================
--- sandbox/gtl/libs/polygon/voronoi_example/voronoi_visualizer.cpp (original)
+++ sandbox/gtl/libs/polygon/voronoi_example/voronoi_visualizer.cpp 2012-01-23 17:22:43 EST (Mon, 23 Jan 2012)
@@ -1,6 +1,6 @@
// Boost.Polygon library voronoi_visualizer.cpp file
-// Copyright Andrii Sydorchuk 2010-2011.
+// Copyright Andrii Sydorchuk 2010-2012.
// Distributed under the Boost Software License, Version 1.0.
// (See accompanying file LICENSE_1_0.txt or copy at
// http://www.boost.org/LICENSE_1_0.txt)
@@ -57,7 +57,8 @@
in_stream.flush();
// Build voronoi diagram.
- construct_voronoi<int>(point_sites, segment_sites, vd_);
+ vd_.clear();
+ construct_voronoi(point_sites, segment_sites, &vd_);
brect_ = voronoi_helper<coordinate_type>::get_view_rectangle(
vd_.bounding_rectangle());