1 // Boost.Polygon library detail/voronoi_ctypes.hpp header file
3 // Copyright Andrii Sydorchuk 2010-2012.
4 // Distributed under the Boost Software License, Version 1.0.
5 // (See accompanying file LICENSE_1_0.txt or copy at
6 // http://www.boost.org/LICENSE_1_0.txt)
8 // See http://www.boost.org for updates, documentation, and revision history.
10 #ifndef BOOST_POLYGON_DETAIL_VORONOI_CTYPES
11 #define BOOST_POLYGON_DETAIL_VORONOI_CTYPES
13 #include <boost/cstdint.hpp>
24 typedef boost::int32_t int32;
25 typedef boost::int64_t int64;
26 typedef boost::uint32_t uint32;
27 typedef boost::uint64_t uint64;
30 // If two floating-point numbers in the same format are ordered (x < y),
31 // then they are ordered the same way when their bits are reinterpreted as
32 // sign-magnitude integers. Values are considered to be almost equal if
33 // their integer bits reinterpretations differ in not more than maxUlps units.
34 template <typename _fpt>
35 struct ulp_comparison;
38 struct ulp_comparison<fpt64> {
45 Result operator()(fpt64 a, fpt64 b, unsigned int maxUlps) const {
48 // Reinterpret double bits as 64-bit signed integer.
49 std::memcpy(&ll_a, &a, sizeof(fpt64));
50 std::memcpy(&ll_b, &b, sizeof(fpt64));
52 // Positive 0.0 is integer zero. Negative 0.0 is 0x8000000000000000.
53 // Map negative zero to an integer zero representation - making it
54 // identical to positive zero - the smallest negative number is
55 // represented by negative one, and downwards from there.
56 if (ll_a < 0x8000000000000000ULL)
57 ll_a = 0x8000000000000000ULL - ll_a;
58 if (ll_b < 0x8000000000000000ULL)
59 ll_b = 0x8000000000000000ULL - ll_b;
61 // Compare 64-bit signed integer representations of input values.
62 // Difference in 1 Ulp is equivalent to a relative error of between
63 // 1/4,000,000,000,000,000 and 1/8,000,000,000,000,000.
65 return (ll_a - ll_b <= maxUlps) ? EQUAL : LESS;
66 return (ll_b - ll_a <= maxUlps) ? EQUAL : MORE;
70 template <typename _fpt>
71 struct extened_exponent_fpt_traits;
74 struct extened_exponent_fpt_traits<fpt64> {
78 MAX_SIGNIFICANT_EXP_DIF = 54
82 // Floating point type wrapper. Allows to extend exponent boundaries to the
83 // integer type range. This class does not handle division by zero, subnormal
85 template <typename _fpt, typename _traits = extened_exponent_fpt_traits<_fpt> >
86 class extended_exponent_fpt {
88 typedef _fpt fpt_type;
89 typedef typename _traits::exp_type exp_type;
91 explicit extended_exponent_fpt(fpt_type val) {
92 val_ = std::frexp(val, &exp_);
95 extended_exponent_fpt(fpt_type val, exp_type exp) {
96 val_ = std::frexp(val, &exp_);
100 bool is_pos() const {
104 bool is_neg() const {
108 bool is_zero() const {
112 extended_exponent_fpt operator-() const {
113 return extended_exponent_fpt(-val_, exp_);
116 extended_exponent_fpt operator+(const extended_exponent_fpt& that) const {
117 if (this->val_ == 0.0 ||
118 that.exp_ > this->exp_ + _traits::MAX_SIGNIFICANT_EXP_DIF) {
121 if (that.val_ == 0.0 ||
122 this->exp_ > that.exp_ + _traits::MAX_SIGNIFICANT_EXP_DIF) {
125 if (this->exp_ >= that.exp_) {
126 exp_type exp_dif = this->exp_ - that.exp_;
127 fpt_type val = std::ldexp(this->val_, exp_dif) + that.val_;
128 return extended_exponent_fpt(val, that.exp_);
130 exp_type exp_dif = that.exp_ - this->exp_;
131 fpt_type val = std::ldexp(that.val_, exp_dif) + this->val_;
132 return extended_exponent_fpt(val, this->exp_);
136 extended_exponent_fpt operator-(const extended_exponent_fpt& that) const {
137 if (this->val_ == 0.0 ||
138 that.exp_ > this->exp_ + _traits::MAX_SIGNIFICANT_EXP_DIF) {
139 return extended_exponent_fpt(-that.val_, that.exp_);
141 if (that.val_ == 0.0 ||
142 this->exp_ > that.exp_ + _traits::MAX_SIGNIFICANT_EXP_DIF) {
145 if (this->exp_ >= that.exp_) {
146 exp_type exp_dif = this->exp_ - that.exp_;
147 fpt_type val = std::ldexp(this->val_, exp_dif) - that.val_;
148 return extended_exponent_fpt(val, that.exp_);
150 exp_type exp_dif = that.exp_ - this->exp_;
151 fpt_type val = std::ldexp(-that.val_, exp_dif) + this->val_;
152 return extended_exponent_fpt(val, this->exp_);
156 extended_exponent_fpt operator*(const extended_exponent_fpt& that) const {
157 fpt_type val = this->val_ * that.val_;
158 exp_type exp = this->exp_ + that.exp_;
159 return extended_exponent_fpt(val, exp);
162 extended_exponent_fpt operator/(const extended_exponent_fpt& that) const {
163 fpt_type val = this->val_ / that.val_;
164 exp_type exp = this->exp_ - that.exp_;
165 return extended_exponent_fpt(val, exp);
168 extended_exponent_fpt& operator+=(const extended_exponent_fpt& that) {
169 return *this = *this + that;
172 extended_exponent_fpt& operator-=(const extended_exponent_fpt& that) {
173 return *this = *this - that;
176 extended_exponent_fpt& operator*=(const extended_exponent_fpt& that) {
177 return *this = *this * that;
180 extended_exponent_fpt& operator/=(const extended_exponent_fpt& that) {
181 return *this = *this / that;
184 extended_exponent_fpt sqrt() const {
191 return extended_exponent_fpt(std::sqrt(val), exp >> 1);
195 return std::ldexp(val_, exp_);
202 typedef extended_exponent_fpt<double> efpt64;
204 template <typename _fpt>
205 extended_exponent_fpt<_fpt> get_sqrt(const extended_exponent_fpt<_fpt>& that) {
209 template <typename _fpt>
210 bool is_pos(const extended_exponent_fpt<_fpt>& that) {
211 return that.is_pos();
214 template <typename _fpt>
215 bool is_neg(const extended_exponent_fpt<_fpt>& that) {
216 return that.is_neg();
219 template <typename _fpt>
220 bool is_zero(const extended_exponent_fpt<_fpt>& that) {
221 return that.is_zero();
224 // Very efficient stack allocated big integer class.
225 // Supports next set of arithmetic operations: +, -, *.
226 template<std::size_t N>
231 extended_int(int32 that) {
233 this->chunks_[0] = that;
235 } else if (that < 0) {
236 this->chunks_[0] = -that;
243 extended_int(int64 that) {
245 this->chunks_[0] = static_cast<uint32>(that);
246 this->chunks_[1] = that >> 32;
247 this->count_ = this->chunks_[1] ? 2 : 1;
248 } else if (that < 0) {
250 this->chunks_[0] = static_cast<uint32>(that);
251 this->chunks_[1] = that >> 32;
252 this->count_ = this->chunks_[1] ? -2 : -1;
258 extended_int(const std::vector<uint32>& chunks, bool plus = true) {
259 this->count_ = static_cast<int32>((std::min)(N, chunks.size()));
260 for (int i = 0; i < this->count_; ++i)
261 this->chunks_[i] = chunks[chunks.size() - i - 1];
263 this->count_ = -this->count_;
266 template<std::size_t M>
267 extended_int(const extended_int<M>& that) {
268 this->count_ = that.count();
269 std::memcpy(this->chunks_, that.chunks(), that.size() * sizeof(uint32));
272 extended_int& operator=(int32 that) {
274 this->chunks_[0] = that;
276 } else if (that < 0) {
277 this->chunks_[0] = -that;
285 extended_int& operator=(int64 that) {
287 this->chunks_[0] = static_cast<uint32>(that);
288 this->chunks_[1] = that >> 32;
289 this->count_ = this->chunks_[1] ? 2 : 1;
290 } else if (that < 0) {
292 this->chunks_[0] = static_cast<uint32>(that);
293 this->chunks_[1] = that >> 32;
294 this->count_ = this->chunks_[1] ? -2 : -1;
301 template<std::size_t M>
302 extended_int& operator=(const extended_int<M>& that) {
303 this->count_ = that.count();
304 std::memcpy(this->chunks_, that.chunks(), that.size() * sizeof(uint32));
308 bool is_pos() const {
309 return this->count_ > 0;
312 bool is_neg() const {
313 return this->count_ < 0;
316 bool is_zero() const {
317 return this->count_ == 0;
320 bool operator==(const extended_int& that) const {
321 if (this->count_ != that.count())
323 for (std::size_t i = 0; i < this->size(); ++i)
324 if (this->chunks_[i] != that.chunks()[i])
329 bool operator!=(const extended_int& that) const {
330 return !(*this == that);
333 bool operator<(const extended_int& that) const {
334 if (this->count_ != that.count())
335 return this->count_ < that.count();
336 std::size_t i = this->size();
341 if (this->chunks_[i] != that.chunks()[i])
342 return (this->chunks_[i] < that.chunks()[i]) ^ (this->count_ < 0);
347 bool operator>(const extended_int& that) const {
351 bool operator<=(const extended_int& that) const {
352 return !(that < *this);
355 bool operator>=(const extended_int& that) const {
356 return !(*this < that);
359 extended_int operator-() const {
360 extended_int ret_val = *this;
366 this->count_ = -this->count_;
369 extended_int operator+(const extended_int& that) const {
370 extended_int ret_val;
371 ret_val.add(*this, that);
375 void add(const extended_int& e1, const extended_int& e2) {
384 if ((e1.count() > 0) ^ (e2.count() > 0)) {
385 dif(e1.chunks(), e1.size(), e2.chunks(), e2.size());
387 add(e1.chunks(), e1.size(), e2.chunks(), e2.size());
390 this->count_ = -this->count_;
393 extended_int operator-(const extended_int& that) const {
394 extended_int ret_val;
395 ret_val.dif(*this, that);
399 void dif(const extended_int& e1, const extended_int& e2) {
402 this->count_ = -this->count_;
409 if ((e1.count() > 0) ^ (e2.count() > 0)) {
410 add(e1.chunks(), e1.size(), e2.chunks(), e2.size());
412 dif(e1.chunks(), e1.size(), e2.chunks(), e2.size());
415 this->count_ = -this->count_;
418 extended_int operator*(int32 that) const {
419 extended_int temp(that);
420 return (*this) * temp;
423 extended_int operator*(int64 that) const {
424 extended_int temp(that);
425 return (*this) * temp;
428 extended_int operator*(const extended_int& that) const {
429 extended_int ret_val;
430 ret_val.mul(*this, that);
434 void mul(const extended_int& e1, const extended_int& e2) {
435 if (!e1.count() || !e2.count()) {
439 mul(e1.chunks(), e1.size(), e2.chunks(), e2.size());
440 if ((e1.count() > 0) ^ (e2.count() > 0))
441 this->count_ = -this->count_;
444 const uint32* chunks() const {
448 int32 count() const {
452 std::size_t size() const {
453 return (std::abs)(count_);
456 std::pair<fpt64, int> p() const {
457 std::pair<fpt64, int> ret_val(0, 0);
458 std::size_t sz = this->size();
463 ret_val.first = static_cast<fpt64>(this->chunks_[0]);
464 } else if (sz == 2) {
465 ret_val.first = static_cast<fpt64>(this->chunks_[1]) *
466 static_cast<fpt64>(0x100000000LL) +
467 static_cast<fpt64>(this->chunks_[0]);
469 for (std::size_t i = 1; i <= 3; ++i) {
470 ret_val.first *= static_cast<fpt64>(0x100000000LL);
471 ret_val.first += static_cast<fpt64>(this->chunks_[sz - i]);
473 ret_val.second = static_cast<int>((sz - 3) << 5);
476 if (this->count_ < 0)
477 ret_val.first = -ret_val.first;
482 std::pair<fpt64, int> p = this->p();
483 return std::ldexp(p.first, p.second);
487 void add(const uint32* c1, std::size_t sz1,
488 const uint32* c2, std::size_t sz2) {
490 add(c2, sz2, c1, sz1);
493 this->count_ = static_cast<int32>(sz1);
495 for (std::size_t i = 0; i < sz2; ++i) {
496 temp += static_cast<uint64>(c1[i]) + static_cast<uint64>(c2[i]);
497 this->chunks_[i] = static_cast<uint32>(temp);
500 for (std::size_t i = sz2; i < sz1; ++i) {
501 temp += static_cast<uint64>(c1[i]);
502 this->chunks_[i] = static_cast<uint32>(temp);
505 if (temp && (this->count_ != N)) {
506 this->chunks_[this->count_] = static_cast<uint32>(temp);
511 void dif(const uint32* c1, std::size_t sz1,
512 const uint32* c2, std::size_t sz2,
515 dif(c2, sz2, c1, sz1, true);
516 this->count_ = -this->count_;
518 } else if ((sz1 == sz2) && !rec) {
521 if (c1[sz1] < c2[sz1]) {
523 dif(c2, sz1, c1, sz1, true);
524 this->count_ = -this->count_;
526 } else if (c1[sz1] > c2[sz1]) {
537 this->count_ = static_cast<int32>(sz1-1);
539 for (std::size_t i = 0; i < sz2; ++i) {
540 this->chunks_[i] = c1[i] - c2[i] - (flag?1:0);
541 flag = (c1[i] < c2[i]) || ((c1[i] == c2[i]) && flag);
543 for (std::size_t i = sz2; i < sz1; ++i) {
544 this->chunks_[i] = c1[i] - (flag?1:0);
545 flag = !c1[i] && flag;
547 if (this->chunks_[this->count_])
551 void mul(const uint32* c1, std::size_t sz1,
552 const uint32* c2, std::size_t sz2) {
553 uint64 cur = 0, nxt, tmp;
554 this->count_ = static_cast<int32>((std::min)(N, sz1 + sz2 - 1));
555 for (std::size_t shift = 0; shift < static_cast<std::size_t>(this->count_);
558 for (std::size_t first = 0; first <= shift; ++first) {
561 std::size_t second = shift - first;
564 tmp = static_cast<uint64>(c1[first]) * static_cast<uint64>(c2[second]);
565 cur += static_cast<uint32>(tmp);
568 this->chunks_[shift] = static_cast<uint32>(cur);
569 cur = nxt + (cur >> 32);
571 if (cur && (this->count_ != N)) {
572 this->chunks_[this->count_] = static_cast<uint32>(cur);
581 template <std::size_t N>
582 bool is_pos(const extended_int<N>& that) {
583 return that.count() > 0;
586 template <std::size_t N>
587 bool is_neg(const extended_int<N>& that) {
588 return that.count() < 0;
591 template <std::size_t N>
592 bool is_zero(const extended_int<N>& that) {
593 return !that.count();
596 struct type_converter_fpt {
597 template <typename T>
598 fpt64 operator()(const T& that) const {
599 return static_cast<fpt64>(that);
602 template <std::size_t N>
603 fpt64 operator()(const extended_int<N>& that) const {
607 fpt64 operator()(const extended_exponent_fpt<fpt64>& that) const {
612 struct type_converter_efpt {
613 template <std::size_t N>
614 extended_exponent_fpt<fpt64> operator()(const extended_int<N>& that) const {
615 std::pair<fpt64, int> p = that.p();
616 return extended_exponent_fpt<fpt64>(p.first, p.second);
620 // Voronoi coordinate type traits make it possible to extend algorithm
621 // input coordinate range to any user provided integer type and algorithm
622 // output coordinate range to any ieee-754 like floating point type.
623 template <typename T>
624 struct voronoi_ctype_traits;
627 struct voronoi_ctype_traits<int32> {
628 typedef int32 int_type;
629 typedef int64 int_x2_type;
630 typedef uint64 uint_x2_type;
631 typedef extended_int<64> big_int_type;
632 typedef fpt64 fpt_type;
633 typedef extended_exponent_fpt<fpt_type> efpt_type;
634 typedef ulp_comparison<fpt_type> ulp_cmp_type;
635 typedef type_converter_fpt to_fpt_converter_type;
636 typedef type_converter_efpt to_efpt_converter_type;
642 #endif // BOOST_POLYGON_DETAIL_VORONOI_CTYPES