ceph/src/boost/boost/polygon/detail/voronoi_ctypes.hpp

   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 <cmath>
  16 #include <cstring>
  17 #include <utility>
  18 #include <vector>
  19
  20 namespace boost {
  21 namespace polygon {
  22 namespace detail {
  23
  24 typedef boost::int32_t int32;
  25 typedef boost::int64_t int64;
  26 typedef boost::uint32_t uint32;
  27 typedef boost::uint64_t uint64;
  28 typedef double fpt64;
  29
  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;
  36
  37 template <>
  38 struct ulp_comparison<fpt64> {
  39   enum Result {
  40     LESS = -1,
  41     EQUAL = 0,
  42     MORE = 1
  43   };
  44
  45   Result operator()(fpt64 a, fpt64 b, unsigned int maxUlps) const {
  46     uint64 ll_a, ll_b;
  47
  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));
  51
  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;
  60
  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.
  64     if (ll_a > ll_b)
  65       return (ll_a - ll_b <= maxUlps) ? EQUAL : LESS;
  66     return (ll_b - ll_a <= maxUlps) ? EQUAL : MORE;
  67   }
  68 };
  69
  70 template <typename _fpt>
  71 struct extened_exponent_fpt_traits;
  72
  73 template <>
  74 struct extened_exponent_fpt_traits<fpt64> {
  75  public:
  76   typedef int exp_type;
  77   enum {
  78     MAX_SIGNIFICANT_EXP_DIF = 54
  79   };
  80 };
  81
  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
  84 // numbers or NaNs.
  85 template <typename _fpt, typename _traits = extened_exponent_fpt_traits<_fpt> >
  86 class extended_exponent_fpt {
  87  public:
  88   typedef _fpt fpt_type;
  89   typedef typename _traits::exp_type exp_type;
  90
  91   explicit extended_exponent_fpt(fpt_type val) {
  92     val_ = std::frexp(val, &exp_);
  93   }
  94
  95   extended_exponent_fpt(fpt_type val, exp_type exp) {
  96     val_ = std::frexp(val, &exp_);
  97     exp_ += exp;
  98   }
  99
 100   bool is_pos() const {
 101     return val_ > 0;
 102   }
 103
 104   bool is_neg() const {
 105     return val_ < 0;
 106   }
 107
 108   bool is_zero() const {
 109     return val_ == 0;
 110   }
 111
 112   extended_exponent_fpt operator-() const {
 113     return extended_exponent_fpt(-val_, exp_);
 114   }
 115
 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) {
 119       return that;
 120     }
 121     if (that.val_ == 0.0 ||
 122         this->exp_ > that.exp_ + _traits::MAX_SIGNIFICANT_EXP_DIF) {
 123       return *this;
 124     }
 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_);
 129     } else {
 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_);
 133     }
 134   }
 135
 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_);
 140     }
 141     if (that.val_ == 0.0 ||
 142         this->exp_ > that.exp_ + _traits::MAX_SIGNIFICANT_EXP_DIF) {
 143       return *this;
 144     }
 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_);
 149     } else {
 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_);
 153     }
 154   }
 155
 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);
 160   }
 161
 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);
 166   }
 167
 168   extended_exponent_fpt& operator+=(const extended_exponent_fpt& that) {
 169     return *this = *this + that;
 170   }
 171
 172   extended_exponent_fpt& operator-=(const extended_exponent_fpt& that) {
 173     return *this = *this - that;
 174   }
 175
 176   extended_exponent_fpt& operator*=(const extended_exponent_fpt& that) {
 177     return *this = *this * that;
 178   }
 179
 180   extended_exponent_fpt& operator/=(const extended_exponent_fpt& that) {
 181     return *this = *this / that;
 182   }
 183
 184   extended_exponent_fpt sqrt() const {
 185     fpt_type val = val_;
 186     exp_type exp = exp_;
 187     if (exp & 1) {
 188       val *= 2.0;
 189       --exp;
 190     }
 191     return extended_exponent_fpt(std::sqrt(val), exp >> 1);
 192   }
 193
 194   fpt_type d() const {
 195     return std::ldexp(val_, exp_);
 196   }
 197
 198  private:
 199   fpt_type val_;
 200   exp_type exp_;
 201 };
 202 typedef extended_exponent_fpt<double> efpt64;
 203
 204 template <typename _fpt>
 205 extended_exponent_fpt<_fpt> get_sqrt(const extended_exponent_fpt<_fpt>& that) {
 206   return that.sqrt();
 207 }
 208
 209 template <typename _fpt>
 210 bool is_pos(const extended_exponent_fpt<_fpt>& that) {
 211   return that.is_pos();
 212 }
 213
 214 template <typename _fpt>
 215 bool is_neg(const extended_exponent_fpt<_fpt>& that) {
 216   return that.is_neg();
 217 }
 218
 219 template <typename _fpt>
 220 bool is_zero(const extended_exponent_fpt<_fpt>& that) {
 221   return that.is_zero();
 222 }
 223
 224 // Very efficient stack allocated big integer class.
 225 // Supports next set of arithmetic operations: +, -, *.
 226 template<std::size_t N>
 227 class extended_int {
 228  public:
 229   extended_int() {}
 230
 231   extended_int(int32 that) {
 232     if (that > 0) {
 233       this->chunks_[0] = that;
 234       this->count_ = 1;
 235     } else if (that < 0) {
 236       this->chunks_[0] = -that;
 237       this->count_ = -1;
 238     } else {
 239       this->count_ = 0;
 240     }
 241   }
 242
 243   extended_int(int64 that) {
 244     if (that > 0) {
 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) {
 249       that = -that;
 250       this->chunks_[0] = static_cast<uint32>(that);
 251       this->chunks_[1] = that >> 32;
 252       this->count_ = this->chunks_[1] ? -2 : -1;
 253     } else {
 254       this->count_ = 0;
 255     }
 256   }
 257
 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];
 262     if (!plus)
 263       this->count_ = -this->count_;
 264   }
 265
 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));
 270   }
 271
 272   extended_int& operator=(int32 that) {
 273     if (that > 0) {
 274       this->chunks_[0] = that;
 275       this->count_ = 1;
 276     } else if (that < 0) {
 277       this->chunks_[0] = -that;
 278       this->count_ = -1;
 279     } else {
 280       this->count_ = 0;
 281     }
 282     return *this;
 283   }
 284
 285   extended_int& operator=(int64 that) {
 286     if (that > 0) {
 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) {
 291       that = -that;
 292       this->chunks_[0] = static_cast<uint32>(that);
 293       this->chunks_[1] = that >> 32;
 294       this->count_ = this->chunks_[1] ? -2 : -1;
 295     } else {
 296       this->count_ = 0;
 297     }
 298     return *this;
 299   }
 300
 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));
 305     return *this;
 306   }
 307
 308   bool is_pos() const {
 309     return this->count_ > 0;
 310   }
 311
 312   bool is_neg() const {
 313     return this->count_ < 0;
 314   }
 315
 316   bool is_zero() const {
 317     return this->count_ == 0;
 318   }
 319
 320   bool operator==(const extended_int& that) const {
 321     if (this->count_ != that.count())
 322       return false;
 323     for (std::size_t i = 0; i < this->size(); ++i)
 324       if (this->chunks_[i] != that.chunks()[i])
 325         return false;
 326     return true;
 327   }
 328
 329   bool operator!=(const extended_int& that) const {
 330     return !(*this == that);
 331   }
 332
 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();
 337     if (!i)
 338       return false;
 339     do {
 340       --i;
 341       if (this->chunks_[i] != that.chunks()[i])
 342         return (this->chunks_[i] < that.chunks()[i]) ^ (this->count_ < 0);
 343     } while (i);
 344     return false;
 345   }
 346
 347   bool operator>(const extended_int& that) const {
 348     return that < *this;
 349   }
 350
 351   bool operator<=(const extended_int& that) const {
 352     return !(that < *this);
 353   }
 354
 355   bool operator>=(const extended_int& that) const {
 356     return !(*this < that);
 357   }
 358
 359   extended_int operator-() const {
 360     extended_int ret_val = *this;
 361     ret_val.neg();
 362     return ret_val;
 363   }
 364
 365   void neg() {
 366     this->count_ = -this->count_;
 367   }
 368
 369   extended_int operator+(const extended_int& that) const {
 370     extended_int ret_val;
 371     ret_val.add(*this, that);
 372     return ret_val;
 373   }
 374
 375   void add(const extended_int& e1, const extended_int& e2) {
 376     if (!e1.count()) {
 377       *this = e2;
 378       return;
 379     }
 380     if (!e2.count()) {
 381       *this = e1;
 382       return;
 383     }
 384     if ((e1.count() > 0) ^ (e2.count() > 0)) {
 385       dif(e1.chunks(), e1.size(), e2.chunks(), e2.size());
 386     } else {
 387       add(e1.chunks(), e1.size(), e2.chunks(), e2.size());
 388     }
 389     if (e1.count() < 0)
 390       this->count_ = -this->count_;
 391   }
 392
 393   extended_int operator-(const extended_int& that) const {
 394     extended_int ret_val;
 395     ret_val.dif(*this, that);
 396     return ret_val;
 397   }
 398
 399   void dif(const extended_int& e1, const extended_int& e2) {
 400     if (!e1.count()) {
 401       *this = e2;
 402       this->count_ = -this->count_;
 403       return;
 404     }
 405     if (!e2.count()) {
 406       *this = e1;
 407       return;
 408     }
 409     if ((e1.count() > 0) ^ (e2.count() > 0)) {
 410       add(e1.chunks(), e1.size(), e2.chunks(), e2.size());
 411     } else {
 412       dif(e1.chunks(), e1.size(), e2.chunks(), e2.size());
 413     }
 414     if (e1.count() < 0)
 415       this->count_ = -this->count_;
 416   }
 417
 418   extended_int operator*(int32 that) const {
 419     extended_int temp(that);
 420     return (*this) * temp;
 421   }
 422
 423   extended_int operator*(int64 that) const {
 424     extended_int temp(that);
 425     return (*this) * temp;
 426   }
 427
 428   extended_int operator*(const extended_int& that) const {
 429     extended_int ret_val;
 430     ret_val.mul(*this, that);
 431     return ret_val;
 432   }
 433
 434   void mul(const extended_int& e1, const extended_int& e2) {
 435     if (!e1.count() || !e2.count()) {
 436       this->count_ = 0;
 437       return;
 438     }
 439     mul(e1.chunks(), e1.size(), e2.chunks(), e2.size());
 440     if ((e1.count() > 0) ^ (e2.count() > 0))
 441       this->count_ = -this->count_;
 442   }
 443
 444   const uint32* chunks() const {
 445     return chunks_;
 446   }
 447
 448   int32 count() const {
 449     return count_;
 450   }
 451
 452   std::size_t size() const {
 453     return (std::abs)(count_);
 454   }
 455
 456   std::pair<fpt64, int> p() const {
 457     std::pair<fpt64, int> ret_val(0, 0);
 458     std::size_t sz = this->size();
 459     if (!sz) {
 460       return ret_val;
 461     } else {
 462       if (sz == 1) {
 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]);
 468       } else {
 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]);
 472         }
 473         ret_val.second = static_cast<int>((sz - 3) << 5);
 474       }
 475     }
 476     if (this->count_ < 0)
 477       ret_val.first = -ret_val.first;
 478     return ret_val;
 479   }
 480
 481   fpt64 d() const {
 482     std::pair<fpt64, int> p = this->p();
 483     return std::ldexp(p.first, p.second);
 484   }
 485
 486  private:
 487   void add(const uint32* c1, std::size_t sz1,
 488            const uint32* c2, std::size_t sz2) {
 489     if (sz1 < sz2) {
 490       add(c2, sz2, c1, sz1);
 491       return;
 492     }
 493     this->count_ = static_cast<int32>(sz1);
 494     uint64 temp = 0;
 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);
 498       temp >>= 32;
 499     }
 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);
 503       temp >>= 32;
 504     }
 505     if (temp && (this->count_ != N)) {
 506       this->chunks_[this->count_] = static_cast<uint32>(temp);
 507       ++this->count_;
 508     }
 509   }
 510
 511   void dif(const uint32* c1, std::size_t sz1,
 512            const uint32* c2, std::size_t sz2,
 513            bool rec = false) {
 514     if (sz1 < sz2) {
 515       dif(c2, sz2, c1, sz1, true);
 516       this->count_ = -this->count_;
 517       return;
 518     } else if ((sz1 == sz2) && !rec) {
 519       do {
 520         --sz1;
 521         if (c1[sz1] < c2[sz1]) {
 522           ++sz1;
 523           dif(c2, sz1, c1, sz1, true);
 524           this->count_ = -this->count_;
 525           return;
 526         } else if (c1[sz1] > c2[sz1]) {
 527           ++sz1;
 528           break;
 529         }
 530       } while (sz1);
 531       if (!sz1) {
 532         this->count_ = 0;
 533         return;
 534       }
 535       sz2 = sz1;
 536     }
 537     this->count_ = static_cast<int32>(sz1-1);
 538     bool flag = false;
 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);
 542     }
 543     for (std::size_t i = sz2; i < sz1; ++i) {
 544       this->chunks_[i] = c1[i] - (flag?1:0);
 545       flag = !c1[i] && flag;
 546     }
 547     if (this->chunks_[this->count_])
 548       ++this->count_;
 549   }
 550
 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_);
 556          ++shift) {
 557       nxt = 0;
 558       for (std::size_t first = 0; first <= shift; ++first) {
 559         if (first >= sz1)
 560           break;
 561         std::size_t second = shift - first;
 562         if (second >= sz2)
 563           continue;
 564         tmp = static_cast<uint64>(c1[first]) * static_cast<uint64>(c2[second]);
 565         cur += static_cast<uint32>(tmp);
 566         nxt += tmp >> 32;
 567       }
 568       this->chunks_[shift] = static_cast<uint32>(cur);
 569       cur = nxt + (cur >> 32);
 570     }
 571     if (cur && (this->count_ != N)) {
 572       this->chunks_[this->count_] = static_cast<uint32>(cur);
 573       ++this->count_;
 574     }
 575   }
 576
 577   uint32 chunks_[N];
 578   int32 count_;
 579 };
 580
 581 template <std::size_t N>
 582 bool is_pos(const extended_int<N>& that) {
 583   return that.count() > 0;
 584 }
 585
 586 template <std::size_t N>
 587 bool is_neg(const extended_int<N>& that) {
 588   return that.count() < 0;
 589 }
 590
 591 template <std::size_t N>
 592 bool is_zero(const extended_int<N>& that) {
 593   return !that.count();
 594 }
 595
 596 struct type_converter_fpt {
 597   template <typename T>
 598   fpt64 operator()(const T& that) const {
 599     return static_cast<fpt64>(that);
 600   }
 601
 602   template <std::size_t N>
 603   fpt64 operator()(const extended_int<N>& that) const {
 604     return that.d();
 605   }
 606
 607   fpt64 operator()(const extended_exponent_fpt<fpt64>& that) const {
 608     return that.d();
 609   }
 610 };
 611
 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);
 617   }
 618 };
 619
 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;
 625
 626 template <>
 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;
 637 };
 638 }  // detail
 639 }  // polygon
 640 }  // boost
 641
 642 #endif  // BOOST_POLYGON_DETAIL_VORONOI_CTYPES