summaryrefslogtreecommitdiff
path: root/boost/polygon/detail/voronoi_ctypes.hpp
diff options
context:
space:
mode:
Diffstat (limited to 'boost/polygon/detail/voronoi_ctypes.hpp')
-rw-r--r--boost/polygon/detail/voronoi_ctypes.hpp642
1 files changed, 642 insertions, 0 deletions
diff --git a/boost/polygon/detail/voronoi_ctypes.hpp b/boost/polygon/detail/voronoi_ctypes.hpp
new file mode 100644
index 0000000000..d8580dda4e
--- /dev/null
+++ b/boost/polygon/detail/voronoi_ctypes.hpp
@@ -0,0 +1,642 @@
+// 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_DETAIL_VORONOI_CTYPES
+#define BOOST_POLYGON_DETAIL_VORONOI_CTYPES
+
+#include <boost/cstdint.hpp>
+
+#include <cmath>
+#include <cstring>
+#include <utility>
+#include <vector>
+
+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 bits reinterpretations differ in not more than maxUlps units.
+template <typename _fpt>
+struct ulp_comparison;
+
+template <>
+struct ulp_comparison<fpt64> {
+ enum Result {
+ LESS = -1,
+ EQUAL = 0,
+ MORE = 1
+ };
+
+ Result operator()(fpt64 a, fpt64 b, unsigned int maxUlps) const {
+ uint64 ll_a, ll_b;
+
+ // Reinterpret double bits as 64-bit signed integer.
+ std::memcpy(&ll_a, &a, sizeof(fpt64));
+ std::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 _fpt>
+struct extened_exponent_fpt_traits;
+
+template <>
+struct extened_exponent_fpt_traits<fpt64> {
+ public:
+ typedef int exp_type;
+ enum {
+ MAX_SIGNIFICANT_EXP_DIF = 54
+ };
+};
+
+// Floating point type wrapper. Allows to extend exponent boundaries to the
+// integer type range. This class does not handle division by zero, subnormal
+// numbers or NaNs.
+template <typename _fpt, typename _traits = extened_exponent_fpt_traits<_fpt> >
+class extended_exponent_fpt {
+ public:
+ typedef _fpt fpt_type;
+ typedef typename _traits::exp_type exp_type;
+
+ explicit extended_exponent_fpt(fpt_type val) {
+ val_ = std::frexp(val, &exp_);
+ }
+
+ extended_exponent_fpt(fpt_type val, exp_type exp) {
+ val_ = std::frexp(val, &exp_);
+ exp_ += exp;
+ }
+
+ bool is_pos() const {
+ return val_ > 0;
+ }
+
+ bool is_neg() const {
+ return val_ < 0;
+ }
+
+ bool is_zero() const {
+ return val_ == 0;
+ }
+
+ extended_exponent_fpt operator-() const {
+ return extended_exponent_fpt(-val_, exp_);
+ }
+
+ extended_exponent_fpt operator+(const extended_exponent_fpt& that) const {
+ if (this->val_ == 0.0 ||
+ that.exp_ > this->exp_ + _traits::MAX_SIGNIFICANT_EXP_DIF) {
+ return that;
+ }
+ if (that.val_ == 0.0 ||
+ this->exp_ > that.exp_ + _traits::MAX_SIGNIFICANT_EXP_DIF) {
+ return *this;
+ }
+ if (this->exp_ >= that.exp_) {
+ exp_type exp_dif = this->exp_ - that.exp_;
+ fpt_type val = std::ldexp(this->val_, exp_dif) + that.val_;
+ return extended_exponent_fpt(val, that.exp_);
+ } else {
+ exp_type exp_dif = that.exp_ - this->exp_;
+ fpt_type val = std::ldexp(that.val_, exp_dif) + this->val_;
+ return extended_exponent_fpt(val, this->exp_);
+ }
+ }
+
+ extended_exponent_fpt operator-(const extended_exponent_fpt& that) const {
+ if (this->val_ == 0.0 ||
+ that.exp_ > this->exp_ + _traits::MAX_SIGNIFICANT_EXP_DIF) {
+ return extended_exponent_fpt(-that.val_, that.exp_);
+ }
+ if (that.val_ == 0.0 ||
+ this->exp_ > that.exp_ + _traits::MAX_SIGNIFICANT_EXP_DIF) {
+ return *this;
+ }
+ if (this->exp_ >= that.exp_) {
+ exp_type exp_dif = this->exp_ - that.exp_;
+ fpt_type val = std::ldexp(this->val_, exp_dif) - that.val_;
+ return extended_exponent_fpt(val, that.exp_);
+ } else {
+ exp_type exp_dif = that.exp_ - this->exp_;
+ fpt_type val = std::ldexp(-that.val_, exp_dif) + this->val_;
+ return extended_exponent_fpt(val, this->exp_);
+ }
+ }
+
+ extended_exponent_fpt operator*(const extended_exponent_fpt& that) const {
+ fpt_type val = this->val_ * that.val_;
+ exp_type exp = this->exp_ + that.exp_;
+ return extended_exponent_fpt(val, exp);
+ }
+
+ extended_exponent_fpt operator/(const extended_exponent_fpt& that) const {
+ fpt_type val = this->val_ / that.val_;
+ exp_type exp = this->exp_ - that.exp_;
+ return extended_exponent_fpt(val, exp);
+ }
+
+ 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 = val_;
+ exp_type exp = exp_;
+ if (exp & 1) {
+ val *= 2.0;
+ --exp;
+ }
+ return extended_exponent_fpt(std::sqrt(val), exp >> 1);
+ }
+
+ fpt_type d() const {
+ return std::ldexp(val_, exp_);
+ }
+
+ private:
+ fpt_type val_;
+ exp_type exp_;
+};
+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();
+}
+
+// Very efficient stack allocated big integer class.
+// Supports next set of arithmetic operations: +, -, *.
+template<std::size_t N>
+class extended_int {
+ public:
+ 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->count_ = 0;
+ }
+ }
+
+ extended_int(int64 that) {
+ if (that > 0) {
+ this->chunks_[0] = static_cast<uint32>(that);
+ this->chunks_[1] = that >> 32;
+ this->count_ = this->chunks_[1] ? 2 : 1;
+ } else if (that < 0) {
+ that = -that;
+ this->chunks_[0] = static_cast<uint32>(that);
+ this->chunks_[1] = that >> 32;
+ this->count_ = this->chunks_[1] ? -2 : -1;
+ } else {
+ this->count_ = 0;
+ }
+ }
+
+ extended_int(const std::vector<uint32>& chunks, bool plus = true) {
+ this->count_ = static_cast<int32>((std::min)(N, chunks.size()));
+ for (int i = 0; i < this->count_; ++i)
+ this->chunks_[i] = chunks[chunks.size() - i - 1];
+ if (!plus)
+ this->count_ = -this->count_;
+ }
+
+ template<std::size_t M>
+ extended_int(const extended_int<M>& that) {
+ this->count_ = that.count();
+ std::memcpy(this->chunks_, that.chunks(), that.size() * sizeof(uint32));
+ }
+
+ extended_int& 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->count_ = 0;
+ }
+ return *this;
+ }
+
+ extended_int& operator=(int64 that) {
+ if (that > 0) {
+ this->chunks_[0] = static_cast<uint32>(that);
+ this->chunks_[1] = that >> 32;
+ this->count_ = this->chunks_[1] ? 2 : 1;
+ } else if (that < 0) {
+ that = -that;
+ this->chunks_[0] = static_cast<uint32>(that);
+ this->chunks_[1] = that >> 32;
+ this->count_ = this->chunks_[1] ? -2 : -1;
+ } else {
+ this->count_ = 0;
+ }
+ return *this;
+ }
+
+ template<std::size_t M>
+ extended_int& operator=(const extended_int<M>& that) {
+ this->count_ = that.count();
+ std::memcpy(this->chunks_, that.chunks(), that.size() * 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;
+ }
+
+ bool operator==(const extended_int& that) const {
+ if (this->count_ != that.count())
+ return false;
+ for (std::size_t i = 0; i < this->size(); ++i)
+ if (this->chunks_[i] != that.chunks()[i])
+ return false;
+ return true;
+ }
+
+ bool operator!=(const extended_int& that) const {
+ return !(*this == that);
+ }
+
+ bool operator<(const extended_int& that) const {
+ if (this->count_ != that.count())
+ return this->count_ < that.count();
+ std::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;
+ }
+
+ bool operator>(const extended_int& that) const {
+ return that < *this;
+ }
+
+ bool operator<=(const extended_int& that) const {
+ return !(that < *this);
+ }
+
+ bool operator>=(const extended_int& that) const {
+ return !(*this < that);
+ }
+
+ extended_int operator-() const {
+ extended_int ret_val = *this;
+ ret_val.neg();
+ return ret_val;
+ }
+
+ void neg() {
+ this->count_ = -this->count_;
+ }
+
+ extended_int operator+(const extended_int& that) const {
+ extended_int ret_val;
+ ret_val.add(*this, that);
+ return ret_val;
+ }
+
+ void add(const extended_int& e1, const extended_int& 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_;
+ }
+
+ extended_int operator-(const extended_int& that) const {
+ extended_int ret_val;
+ ret_val.dif(*this, that);
+ return ret_val;
+ }
+
+ void dif(const extended_int& e1, const extended_int& 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 operator*(int32 that) const {
+ extended_int temp(that);
+ return (*this) * temp;
+ }
+
+ extended_int operator*(int64 that) const {
+ extended_int temp(that);
+ return (*this) * temp;
+ }
+
+ extended_int operator*(const extended_int& that) const {
+ extended_int ret_val;
+ ret_val.mul(*this, that);
+ return ret_val;
+ }
+
+ void mul(const extended_int& e1, const extended_int& e2) {
+ if (!e1.count() || !e2.count()) {
+ 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_;
+ }
+
+ std::size_t size() const {
+ return (std::abs)(count_);
+ }
+
+ std::pair<fpt64, int> p() const {
+ std::pair<fpt64, int> ret_val(0, 0);
+ std::size_t sz = this->size();
+ if (!sz) {
+ return ret_val;
+ } else {
+ if (sz == 1) {
+ ret_val.first = static_cast<fpt64>(this->chunks_[0]);
+ } else if (sz == 2) {
+ ret_val.first = static_cast<fpt64>(this->chunks_[1]) *
+ static_cast<fpt64>(0x100000000LL) +
+ static_cast<fpt64>(this->chunks_[0]);
+ } else {
+ for (std::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, int> p = this->p();
+ return std::ldexp(p.first, p.second);
+ }
+
+ private:
+ void add(const uint32* c1, std::size_t sz1,
+ const uint32* c2, std::size_t sz2) {
+ if (sz1 < sz2) {
+ add(c2, sz2, c1, sz1);
+ return;
+ }
+ this->count_ = sz1;
+ uint64 temp = 0;
+ for (std::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);
+ temp >>= 32;
+ }
+ for (std::size_t i = sz2; i < sz1; ++i) {
+ temp += static_cast<uint64>(c1[i]);
+ this->chunks_[i] = static_cast<uint32>(temp);
+ temp >>= 32;
+ }
+ if (temp && (this->count_ != N)) {
+ this->chunks_[this->count_] = static_cast<uint32>(temp);
+ ++this->count_;
+ }
+ }
+
+ void dif(const uint32* c1, std::size_t sz1,
+ const uint32* c2, std::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;
+ }
+ } while (sz1);
+ if (!sz1) {
+ this->count_ = 0;
+ return;
+ }
+ sz2 = sz1;
+ }
+ this->count_ = sz1-1;
+ bool flag = false;
+ for (std::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 (std::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, std::size_t sz1,
+ const uint32* c2, std::size_t sz2) {
+ uint64 cur = 0, nxt, tmp;
+ this->count_ = static_cast<int32>((std::min)(N, sz1 + sz2 - 1));
+ for (std::size_t shift = 0; shift < static_cast<std::size_t>(this->count_);
+ ++shift) {
+ nxt = 0;
+ for (std::size_t first = 0; first <= shift; ++first) {
+ if (first >= sz1)
+ break;
+ std::size_t second = shift - first;
+ if (second >= sz2)
+ continue;
+ tmp = static_cast<uint64>(c1[first]) * static_cast<uint64>(c2[second]);
+ cur += static_cast<uint32>(tmp);
+ nxt += tmp >> 32;
+ }
+ this->chunks_[shift] = static_cast<uint32>(cur);
+ cur = nxt + (cur >> 32);
+ }
+ if (cur && (this->count_ != N)) {
+ this->chunks_[this->count_] = static_cast<uint32>(cur);
+ ++this->count_;
+ }
+ }
+
+ uint32 chunks_[N];
+ int32 count_;
+};
+
+template <std::size_t N>
+bool is_pos(const extended_int<N>& that) {
+ return that.count() > 0;
+}
+
+template <std::size_t N>
+bool is_neg(const extended_int<N>& that) {
+ return that.count() < 0;
+}
+
+template <std::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 <std::size_t N>
+ fpt64 operator()(const extended_int<N>& that) const {
+ return that.d();
+ }
+
+ fpt64 operator()(const extended_exponent_fpt<fpt64>& that) const {
+ return that.d();
+ }
+};
+
+struct type_converter_efpt {
+ template <std::size_t N>
+ extended_exponent_fpt<fpt64> operator()(const extended_int<N>& that) const {
+ std::pair<fpt64, int> p = that.p();
+ return extended_exponent_fpt<fpt64>(p.first, p.second);
+ }
+};
+
+// Voronoi coordinate type traits make it possible to extend algorithm
+// input coordinate range to any user provided integer type and algorithm
+// output coordinate range to any ieee-754 like floating point type.
+template <typename T>
+struct voronoi_ctype_traits;
+
+template <>
+struct voronoi_ctype_traits<int32> {
+ typedef int32 int_type;
+ typedef int64 int_x2_type;
+ typedef uint64 uint_x2_type;
+ typedef extended_int<64> 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 // BOOST_POLYGON_DETAIL_VORONOI_CTYPES
generated by cgit v1.2.3 (git 2.25.1) at 2024-06-07 22:58:38 +0000