1 // Boost.Polygon library detail/voronoi_ctypes.hpp header file
2
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)
7
8 // See http://www.boost.org for updates, documentation, and revision history.
9
10 #ifndef BOOST_POLYGON_DETAIL_VORONOI_CTYPES
11 #define BOOST_POLYGON_DETAIL_VORONOI_CTYPES
12
13 #include <boost/cstdint.hpp>
14
15 #include <algorithm>
16 #include <cmath>
17 #include <cstring>
18 #include <utility>
19 #include <vector>
20
21 namespace boost {
22 namespace polygon {
23 namespace detail {
24
25 typedef boost::int32_t int32;
26 typedef boost::int64_t int64;
27 typedef boost::uint32_t uint32;
28 typedef boost::uint64_t uint64;
29 typedef double fpt64;
30
31 // If two floating-point numbers in the same format are ordered (x < y),
32 // then they are ordered the same way when their bits are reinterpreted as
33 // sign-magnitude integers. Values are considered to be almost equal if
34 // their integer bits reinterpretations differ in not more than maxUlps units.
35 template <typename _fpt>
36 struct ulp_comparison;
37
38 template <>
39 struct ulp_comparison<fpt64> {
40 enum Result {
41 LESS = -1,
42 EQUAL = 0,
43 MORE = 1
44 };
45
operator ()boost::polygon::detail::ulp_comparison46 Result operator()(fpt64 a, fpt64 b, unsigned int maxUlps) const {
47 uint64 ll_a, ll_b;
48
49 // Reinterpret double bits as 64-bit signed integer.
50 std::memcpy(&ll_a, &a, sizeof(fpt64));
51 std::memcpy(&ll_b, &b, sizeof(fpt64));
52
53 // Positive 0.0 is integer zero. Negative 0.0 is 0x8000000000000000.
54 // Map negative zero to an integer zero representation - making it
55 // identical to positive zero - the smallest negative number is
56 // represented by negative one, and downwards from there.
57 if (ll_a < 0x8000000000000000ULL)
58 ll_a = 0x8000000000000000ULL - ll_a;
59 if (ll_b < 0x8000000000000000ULL)
60 ll_b = 0x8000000000000000ULL - ll_b;
61
62 // Compare 64-bit signed integer representations of input values.
63 // Difference in 1 Ulp is equivalent to a relative error of between
64 // 1/4,000,000,000,000,000 and 1/8,000,000,000,000,000.
65 if (ll_a > ll_b)
66 return (ll_a - ll_b <= maxUlps) ? EQUAL : LESS;
67 return (ll_b - ll_a <= maxUlps) ? EQUAL : MORE;
68 }
69 };
70
71 template <typename _fpt>
72 struct extened_exponent_fpt_traits;
73
74 template <>
75 struct extened_exponent_fpt_traits<fpt64> {
76 public:
77 typedef int exp_type;
78 enum {
79 MAX_SIGNIFICANT_EXP_DIF = 54
80 };
81 };
82
83 // Floating point type wrapper. Allows to extend exponent boundaries to the
84 // integer type range. This class does not handle division by zero, subnormal
85 // numbers or NaNs.
86 template <typename _fpt, typename _traits = extened_exponent_fpt_traits<_fpt> >
87 class extended_exponent_fpt {
88 public:
89 typedef _fpt fpt_type;
90 typedef typename _traits::exp_type exp_type;
91
extended_exponent_fpt(fpt_type val)92 explicit extended_exponent_fpt(fpt_type val) {
93 val_ = std::frexp(val, &exp_);
94 }
95
extended_exponent_fpt(fpt_type val,exp_type exp)96 extended_exponent_fpt(fpt_type val, exp_type exp) {
97 val_ = std::frexp(val, &exp_);
98 exp_ += exp;
99 }
100
is_pos() const101 bool is_pos() const {
102 return val_ > 0;
103 }
104
is_neg() const105 bool is_neg() const {
106 return val_ < 0;
107 }
108
is_zero() const109 bool is_zero() const {
110 return val_ == 0;
111 }
112
operator -() const113 extended_exponent_fpt operator-() const {
114 return extended_exponent_fpt(-val_, exp_);
115 }
116
operator +(const extended_exponent_fpt & that) const117 extended_exponent_fpt operator+(const extended_exponent_fpt& that) const {
118 if (this->val_ == 0.0 ||
119 that.exp_ > this->exp_ + _traits::MAX_SIGNIFICANT_EXP_DIF) {
120 return that;
121 }
122 if (that.val_ == 0.0 ||
123 this->exp_ > that.exp_ + _traits::MAX_SIGNIFICANT_EXP_DIF) {
124 return *this;
125 }
126 if (this->exp_ >= that.exp_) {
127 exp_type exp_dif = this->exp_ - that.exp_;
128 fpt_type val = std::ldexp(this->val_, exp_dif) + that.val_;
129 return extended_exponent_fpt(val, that.exp_);
130 } else {
131 exp_type exp_dif = that.exp_ - this->exp_;
132 fpt_type val = std::ldexp(that.val_, exp_dif) + this->val_;
133 return extended_exponent_fpt(val, this->exp_);
134 }
135 }
136
operator -(const extended_exponent_fpt & that) const137 extended_exponent_fpt operator-(const extended_exponent_fpt& that) const {
138 if (this->val_ == 0.0 ||
139 that.exp_ > this->exp_ + _traits::MAX_SIGNIFICANT_EXP_DIF) {
140 return extended_exponent_fpt(-that.val_, that.exp_);
141 }
142 if (that.val_ == 0.0 ||
143 this->exp_ > that.exp_ + _traits::MAX_SIGNIFICANT_EXP_DIF) {
144 return *this;
145 }
146 if (this->exp_ >= that.exp_) {
147 exp_type exp_dif = this->exp_ - that.exp_;
148 fpt_type val = std::ldexp(this->val_, exp_dif) - that.val_;
149 return extended_exponent_fpt(val, that.exp_);
150 } else {
151 exp_type exp_dif = that.exp_ - this->exp_;
152 fpt_type val = std::ldexp(-that.val_, exp_dif) + this->val_;
153 return extended_exponent_fpt(val, this->exp_);
154 }
155 }
156
operator *(const extended_exponent_fpt & that) const157 extended_exponent_fpt operator*(const extended_exponent_fpt& that) const {
158 fpt_type val = this->val_ * that.val_;
159 exp_type exp = this->exp_ + that.exp_;
160 return extended_exponent_fpt(val, exp);
161 }
162
operator /(const extended_exponent_fpt & that) const163 extended_exponent_fpt operator/(const extended_exponent_fpt& that) const {
164 fpt_type val = this->val_ / that.val_;
165 exp_type exp = this->exp_ - that.exp_;
166 return extended_exponent_fpt(val, exp);
167 }
168
operator +=(const extended_exponent_fpt & that)169 extended_exponent_fpt& operator+=(const extended_exponent_fpt& that) {
170 return *this = *this + that;
171 }
172
operator -=(const extended_exponent_fpt & that)173 extended_exponent_fpt& operator-=(const extended_exponent_fpt& that) {
174 return *this = *this - that;
175 }
176
operator *=(const extended_exponent_fpt & that)177 extended_exponent_fpt& operator*=(const extended_exponent_fpt& that) {
178 return *this = *this * that;
179 }
180
operator /=(const extended_exponent_fpt & that)181 extended_exponent_fpt& operator/=(const extended_exponent_fpt& that) {
182 return *this = *this / that;
183 }
184
sqrt() const185 extended_exponent_fpt sqrt() const {
186 fpt_type val = val_;
187 exp_type exp = exp_;
188 if (exp & 1) {
189 val *= 2.0;
190 --exp;
191 }
192 return extended_exponent_fpt(std::sqrt(val), exp >> 1);
193 }
194
d() const195 fpt_type d() const {
196 return std::ldexp(val_, exp_);
197 }
198
199 private:
200 fpt_type val_;
201 exp_type exp_;
202 };
203 typedef extended_exponent_fpt<double> efpt64;
204
205 template <typename _fpt>
get_sqrt(const extended_exponent_fpt<_fpt> & that)206 extended_exponent_fpt<_fpt> get_sqrt(const extended_exponent_fpt<_fpt>& that) {
207 return that.sqrt();
208 }
209
210 template <typename _fpt>
is_pos(const extended_exponent_fpt<_fpt> & that)211 bool is_pos(const extended_exponent_fpt<_fpt>& that) {
212 return that.is_pos();
213 }
214
215 template <typename _fpt>
is_neg(const extended_exponent_fpt<_fpt> & that)216 bool is_neg(const extended_exponent_fpt<_fpt>& that) {
217 return that.is_neg();
218 }
219
220 template <typename _fpt>
is_zero(const extended_exponent_fpt<_fpt> & that)221 bool is_zero(const extended_exponent_fpt<_fpt>& that) {
222 return that.is_zero();
223 }
224
225 // Very efficient stack allocated big integer class.
226 // Supports next set of arithmetic operations: +, -, *.
227 template<std::size_t N>
228 class extended_int {
229 public:
extended_int()230 extended_int() {}
231
extended_int(int32 that)232 extended_int(int32 that) {
233 if (that > 0) {
234 this->chunks_[0] = that;
235 this->count_ = 1;
236 } else if (that < 0) {
237 this->chunks_[0] = -that;
238 this->count_ = -1;
239 } else {
240 this->count_ = 0;
241 }
242 }
243
extended_int(int64 that)244 extended_int(int64 that) {
245 if (that > 0) {
246 this->chunks_[0] = static_cast<uint32>(that);
247 this->chunks_[1] = that >> 32;
248 this->count_ = this->chunks_[1] ? 2 : 1;
249 } else if (that < 0) {
250 that = -that;
251 this->chunks_[0] = static_cast<uint32>(that);
252 this->chunks_[1] = that >> 32;
253 this->count_ = this->chunks_[1] ? -2 : -1;
254 } else {
255 this->count_ = 0;
256 }
257 }
258
extended_int(const std::vector<uint32> & chunks,bool plus=true)259 extended_int(const std::vector<uint32>& chunks, bool plus = true) {
260 this->count_ = static_cast<int32>((std::min)(N, chunks.size()));
261 for (int i = 0; i < this->count_; ++i)
262 this->chunks_[i] = chunks[chunks.size() - i - 1];
263 if (!plus)
264 this->count_ = -this->count_;
265 }
266
267 template<std::size_t M>
extended_int(const extended_int<M> & that)268 extended_int(const extended_int<M>& that) {
269 this->count_ = that.count();
270 std::memcpy(this->chunks_, that.chunks(), that.size() * sizeof(uint32));
271 }
272
operator =(int32 that)273 extended_int& operator=(int32 that) {
274 if (that > 0) {
275 this->chunks_[0] = that;
276 this->count_ = 1;
277 } else if (that < 0) {
278 this->chunks_[0] = -that;
279 this->count_ = -1;
280 } else {
281 this->count_ = 0;
282 }
283 return *this;
284 }
285
operator =(int64 that)286 extended_int& operator=(int64 that) {
287 if (that > 0) {
288 this->chunks_[0] = static_cast<uint32>(that);
289 this->chunks_[1] = that >> 32;
290 this->count_ = this->chunks_[1] ? 2 : 1;
291 } else if (that < 0) {
292 that = -that;
293 this->chunks_[0] = static_cast<uint32>(that);
294 this->chunks_[1] = that >> 32;
295 this->count_ = this->chunks_[1] ? -2 : -1;
296 } else {
297 this->count_ = 0;
298 }
299 return *this;
300 }
301
302 template<std::size_t M>
operator =(const extended_int<M> & that)303 extended_int& operator=(const extended_int<M>& that) {
304 this->count_ = that.count();
305 std::memcpy(this->chunks_, that.chunks(), that.size() * sizeof(uint32));
306 return *this;
307 }
308
is_pos() const309 bool is_pos() const {
310 return this->count_ > 0;
311 }
312
is_neg() const313 bool is_neg() const {
314 return this->count_ < 0;
315 }
316
is_zero() const317 bool is_zero() const {
318 return this->count_ == 0;
319 }
320
operator ==(const extended_int & that) const321 bool operator==(const extended_int& that) const {
322 if (this->count_ != that.count())
323 return false;
324 for (std::size_t i = 0; i < this->size(); ++i)
325 if (this->chunks_[i] != that.chunks()[i])
326 return false;
327 return true;
328 }
329
operator !=(const extended_int & that) const330 bool operator!=(const extended_int& that) const {
331 return !(*this == that);
332 }
333
operator <(const extended_int & that) const334 bool operator<(const extended_int& that) const {
335 if (this->count_ != that.count())
336 return this->count_ < that.count();
337 std::size_t i = this->size();
338 if (!i)
339 return false;
340 do {
341 --i;
342 if (this->chunks_[i] != that.chunks()[i])
343 return (this->chunks_[i] < that.chunks()[i]) ^ (this->count_ < 0);
344 } while (i);
345 return false;
346 }
347
operator >(const extended_int & that) const348 bool operator>(const extended_int& that) const {
349 return that < *this;
350 }
351
operator <=(const extended_int & that) const352 bool operator<=(const extended_int& that) const {
353 return !(that < *this);
354 }
355
operator >=(const extended_int & that) const356 bool operator>=(const extended_int& that) const {
357 return !(*this < that);
358 }
359
operator -() const360 extended_int operator-() const {
361 extended_int ret_val = *this;
362 ret_val.neg();
363 return ret_val;
364 }
365
neg()366 void neg() {
367 this->count_ = -this->count_;
368 }
369
operator +(const extended_int & that) const370 extended_int operator+(const extended_int& that) const {
371 extended_int ret_val;
372 ret_val.add(*this, that);
373 return ret_val;
374 }
375
add(const extended_int & e1,const extended_int & e2)376 void add(const extended_int& e1, const extended_int& e2) {
377 if (!e1.count()) {
378 *this = e2;
379 return;
380 }
381 if (!e2.count()) {
382 *this = e1;
383 return;
384 }
385 if ((e1.count() > 0) ^ (e2.count() > 0)) {
386 dif(e1.chunks(), e1.size(), e2.chunks(), e2.size());
387 } else {
388 add(e1.chunks(), e1.size(), e2.chunks(), e2.size());
389 }
390 if (e1.count() < 0)
391 this->count_ = -this->count_;
392 }
393
operator -(const extended_int & that) const394 extended_int operator-(const extended_int& that) const {
395 extended_int ret_val;
396 ret_val.dif(*this, that);
397 return ret_val;
398 }
399
dif(const extended_int & e1,const extended_int & e2)400 void dif(const extended_int& e1, const extended_int& e2) {
401 if (!e1.count()) {
402 *this = e2;
403 this->count_ = -this->count_;
404 return;
405 }
406 if (!e2.count()) {
407 *this = e1;
408 return;
409 }
410 if ((e1.count() > 0) ^ (e2.count() > 0)) {
411 add(e1.chunks(), e1.size(), e2.chunks(), e2.size());
412 } else {
413 dif(e1.chunks(), e1.size(), e2.chunks(), e2.size());
414 }
415 if (e1.count() < 0)
416 this->count_ = -this->count_;
417 }
418
operator *(int32 that) const419 extended_int operator*(int32 that) const {
420 extended_int temp(that);
421 return (*this) * temp;
422 }
423
operator *(int64 that) const424 extended_int operator*(int64 that) const {
425 extended_int temp(that);
426 return (*this) * temp;
427 }
428
operator *(const extended_int & that) const429 extended_int operator*(const extended_int& that) const {
430 extended_int ret_val;
431 ret_val.mul(*this, that);
432 return ret_val;
433 }
434
mul(const extended_int & e1,const extended_int & e2)435 void mul(const extended_int& e1, const extended_int& e2) {
436 if (!e1.count() || !e2.count()) {
437 this->count_ = 0;
438 return;
439 }
440 mul(e1.chunks(), e1.size(), e2.chunks(), e2.size());
441 if ((e1.count() > 0) ^ (e2.count() > 0))
442 this->count_ = -this->count_;
443 }
444
chunks() const445 const uint32* chunks() const {
446 return chunks_;
447 }
448
count() const449 int32 count() const {
450 return count_;
451 }
452
size() const453 std::size_t size() const {
454 return (std::abs)(count_);
455 }
456
p() const457 std::pair<fpt64, int> p() const {
458 std::pair<fpt64, int> ret_val(0, 0);
459 std::size_t sz = this->size();
460 if (!sz) {
461 return ret_val;
462 } else {
463 if (sz == 1) {
464 ret_val.first = static_cast<fpt64>(this->chunks_[0]);
465 } else if (sz == 2) {
466 ret_val.first = static_cast<fpt64>(this->chunks_[1]) *
467 static_cast<fpt64>(0x100000000LL) +
468 static_cast<fpt64>(this->chunks_[0]);
469 } else {
470 for (std::size_t i = 1; i <= 3; ++i) {
471 ret_val.first *= static_cast<fpt64>(0x100000000LL);
472 ret_val.first += static_cast<fpt64>(this->chunks_[sz - i]);
473 }
474 ret_val.second = static_cast<int>((sz - 3) << 5);
475 }
476 }
477 if (this->count_ < 0)
478 ret_val.first = -ret_val.first;
479 return ret_val;
480 }
481
d() const482 fpt64 d() const {
483 std::pair<fpt64, int> p = this->p();
484 return std::ldexp(p.first, p.second);
485 }
486
487 private:
add(const uint32 * c1,std::size_t sz1,const uint32 * c2,std::size_t sz2)488 void add(const uint32* c1, std::size_t sz1,
489 const uint32* c2, std::size_t sz2) {
490 if (sz1 < sz2) {
491 add(c2, sz2, c1, sz1);
492 return;
493 }
494 this->count_ = static_cast<int32>(sz1);
495 uint64 temp = 0;
496 for (std::size_t i = 0; i < sz2; ++i) {
497 temp += static_cast<uint64>(c1[i]) + static_cast<uint64>(c2[i]);
498 this->chunks_[i] = static_cast<uint32>(temp);
499 temp >>= 32;
500 }
501 for (std::size_t i = sz2; i < sz1; ++i) {
502 temp += static_cast<uint64>(c1[i]);
503 this->chunks_[i] = static_cast<uint32>(temp);
504 temp >>= 32;
505 }
506 if (temp && (this->count_ != N)) {
507 this->chunks_[this->count_] = static_cast<uint32>(temp);
508 ++this->count_;
509 }
510 }
511
dif(const uint32 * c1,std::size_t sz1,const uint32 * c2,std::size_t sz2,bool rec=false)512 void dif(const uint32* c1, std::size_t sz1,
513 const uint32* c2, std::size_t sz2,
514 bool rec = false) {
515 if (sz1 < sz2) {
516 dif(c2, sz2, c1, sz1, true);
517 this->count_ = -this->count_;
518 return;
519 } else if ((sz1 == sz2) && !rec) {
520 do {
521 --sz1;
522 if (c1[sz1] < c2[sz1]) {
523 ++sz1;
524 dif(c2, sz1, c1, sz1, true);
525 this->count_ = -this->count_;
526 return;
527 } else if (c1[sz1] > c2[sz1]) {
528 ++sz1;
529 break;
530 }
531 } while (sz1);
532 if (!sz1) {
533 this->count_ = 0;
534 return;
535 }
536 sz2 = sz1;
537 }
538 this->count_ = static_cast<int32>(sz1-1);
539 bool flag = false;
540 for (std::size_t i = 0; i < sz2; ++i) {
541 this->chunks_[i] = c1[i] - c2[i] - (flag?1:0);
542 flag = (c1[i] < c2[i]) || ((c1[i] == c2[i]) && flag);
543 }
544 for (std::size_t i = sz2; i < sz1; ++i) {
545 this->chunks_[i] = c1[i] - (flag?1:0);
546 flag = !c1[i] && flag;
547 }
548 if (this->chunks_[this->count_])
549 ++this->count_;
550 }
551
mul(const uint32 * c1,std::size_t sz1,const uint32 * c2,std::size_t sz2)552 void mul(const uint32* c1, std::size_t sz1,
553 const uint32* c2, std::size_t sz2) {
554 uint64 cur = 0, nxt, tmp;
555 this->count_ = static_cast<int32>((std::min)(N, sz1 + sz2 - 1));
556 for (std::size_t shift = 0; shift < static_cast<std::size_t>(this->count_);
557 ++shift) {
558 nxt = 0;
559 for (std::size_t first = 0; first <= shift; ++first) {
560 if (first >= sz1)
561 break;
562 std::size_t second = shift - first;
563 if (second >= sz2)
564 continue;
565 tmp = static_cast<uint64>(c1[first]) * static_cast<uint64>(c2[second]);
566 cur += static_cast<uint32>(tmp);
567 nxt += tmp >> 32;
568 }
569 this->chunks_[shift] = static_cast<uint32>(cur);
570 cur = nxt + (cur >> 32);
571 }
572 if (cur && (this->count_ != N)) {
573 this->chunks_[this->count_] = static_cast<uint32>(cur);
574 ++this->count_;
575 }
576 }
577
578 uint32 chunks_[N];
579 int32 count_;
580 };
581
582 template <std::size_t N>
is_pos(const extended_int<N> & that)583 bool is_pos(const extended_int<N>& that) {
584 return that.count() > 0;
585 }
586
587 template <std::size_t N>
is_neg(const extended_int<N> & that)588 bool is_neg(const extended_int<N>& that) {
589 return that.count() < 0;
590 }
591
592 template <std::size_t N>
is_zero(const extended_int<N> & that)593 bool is_zero(const extended_int<N>& that) {
594 return !that.count();
595 }
596
597 struct type_converter_fpt {
598 template <typename T>
operator ()boost::polygon::detail::type_converter_fpt599 fpt64 operator()(const T& that) const {
600 return static_cast<fpt64>(that);
601 }
602
603 template <std::size_t N>
operator ()boost::polygon::detail::type_converter_fpt604 fpt64 operator()(const extended_int<N>& that) const {
605 return that.d();
606 }
607
operator ()boost::polygon::detail::type_converter_fpt608 fpt64 operator()(const extended_exponent_fpt<fpt64>& that) const {
609 return that.d();
610 }
611 };
612
613 struct type_converter_efpt {
614 template <std::size_t N>
operator ()boost::polygon::detail::type_converter_efpt615 extended_exponent_fpt<fpt64> operator()(const extended_int<N>& that) const {
616 std::pair<fpt64, int> p = that.p();
617 return extended_exponent_fpt<fpt64>(p.first, p.second);
618 }
619 };
620
621 // Voronoi coordinate type traits make it possible to extend algorithm
622 // input coordinate range to any user provided integer type and algorithm
623 // output coordinate range to any ieee-754 like floating point type.
624 template <typename T>
625 struct voronoi_ctype_traits;
626
627 template <>
628 struct voronoi_ctype_traits<int32> {
629 typedef int32 int_type;
630 typedef int64 int_x2_type;
631 typedef uint64 uint_x2_type;
632 typedef extended_int<64> big_int_type;
633 typedef fpt64 fpt_type;
634 typedef extended_exponent_fpt<fpt_type> efpt_type;
635 typedef ulp_comparison<fpt_type> ulp_cmp_type;
636 typedef type_converter_fpt to_fpt_converter_type;
637 typedef type_converter_efpt to_efpt_converter_type;
638 };
639 } // detail
640 } // polygon
641 } // boost
642
643 #endif // BOOST_POLYGON_DETAIL_VORONOI_CTYPES
644