ceph/src/boost/libs/polygon/include/boost/polygon/detail/voronoi_predicates.hpp

   1 // Boost.Polygon library detail/voronoi_predicates.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_PREDICATES
  11 #define BOOST_POLYGON_DETAIL_VORONOI_PREDICATES
  12
  13 #include <utility>
  14
  15 #include "voronoi_robust_fpt.hpp"
  16
  17 namespace boost {
  18 namespace polygon {
  19 namespace detail {
  20
  21 // Predicate utilities. Operates with the coordinate types that could
  22 // be converted to the 32-bit signed integer without precision loss.
  23 template <typename CTYPE_TRAITS>
  24 class voronoi_predicates {
  25  public:
  26   typedef typename CTYPE_TRAITS::int_type int_type;
  27   typedef typename CTYPE_TRAITS::int_x2_type int_x2_type;
  28   typedef typename CTYPE_TRAITS::uint_x2_type uint_x2_type;
  29   typedef typename CTYPE_TRAITS::big_int_type big_int_type;
  30   typedef typename CTYPE_TRAITS::fpt_type fpt_type;
  31   typedef typename CTYPE_TRAITS::efpt_type efpt_type;
  32   typedef typename CTYPE_TRAITS::ulp_cmp_type ulp_cmp_type;
  33   typedef typename CTYPE_TRAITS::to_fpt_converter_type to_fpt_converter;
  34   typedef typename CTYPE_TRAITS::to_efpt_converter_type to_efpt_converter;
  35
  36   enum {
  37     ULPS = 64,
  38     ULPSx2 = 128
  39   };
  40
  41   template <typename Point>
  42   static bool is_vertical(const Point& point1, const Point& point2) {
  43     return point1.x() == point2.x();
  44   }
  45
  46   template <typename Site>
  47   static bool is_vertical(const Site& site) {
  48     return is_vertical(site.point0(), site.point1());
  49   }
  50
  51   // Compute robust cross_product: a1 * b2 - b1 * a2.
  52   // It was mathematically proven that the result is correct
  53   // with epsilon relative error equal to 1EPS.
  54   static fpt_type robust_cross_product(int_x2_type a1_,
  55                                        int_x2_type b1_,
  56                                        int_x2_type a2_,
  57                                        int_x2_type b2_) {
  58     static to_fpt_converter to_fpt;
  59     uint_x2_type a1 = static_cast<uint_x2_type>(is_neg(a1_) ? -a1_ : a1_);
  60     uint_x2_type b1 = static_cast<uint_x2_type>(is_neg(b1_) ? -b1_ : b1_);
  61     uint_x2_type a2 = static_cast<uint_x2_type>(is_neg(a2_) ? -a2_ : a2_);
  62     uint_x2_type b2 = static_cast<uint_x2_type>(is_neg(b2_) ? -b2_ : b2_);
  63
  64     uint_x2_type l = a1 * b2;
  65     uint_x2_type r = b1 * a2;
  66
  67     if (is_neg(a1_) ^ is_neg(b2_)) {
  68       if (is_neg(a2_) ^ is_neg(b1_))
  69         return (l > r) ? -to_fpt(l - r) : to_fpt(r - l);
  70       else
  71         return -to_fpt(l + r);
  72     } else {
  73       if (is_neg(a2_) ^ is_neg(b1_))
  74         return to_fpt(l + r);
  75       else
  76         return (l < r) ? -to_fpt(r - l) : to_fpt(l - r);
  77     }
  78   }
  79
  80   typedef struct orientation_test {
  81    public:
  82     // Represents orientation test result.
  83     enum Orientation {
  84       RIGHT = -1,
  85       COLLINEAR = 0,
  86       LEFT = 1
  87     };
  88
  89     // Value is a determinant of two vectors (e.g. x1 * y2 - x2 * y1).
  90     // Return orientation based on the sign of the determinant.
  91     template <typename T>
  92     static Orientation eval(T value) {
  93       if (is_zero(value)) return COLLINEAR;
  94       return (is_neg(value)) ? RIGHT : LEFT;
  95     }
  96
  97     static Orientation eval(int_x2_type dif_x1_,
  98                             int_x2_type dif_y1_,
  99                             int_x2_type dif_x2_,
 100                             int_x2_type dif_y2_) {
 101       return eval(robust_cross_product(dif_x1_, dif_y1_, dif_x2_, dif_y2_));
 102     }
 103
 104     template <typename Point>
 105     static Orientation eval(const Point& point1,
 106                             const Point& point2,
 107                             const Point& point3) {
 108       int_x2_type dx1 = static_cast<int_x2_type>(point1.x()) -
 109                         static_cast<int_x2_type>(point2.x());
 110       int_x2_type dx2 = static_cast<int_x2_type>(point2.x()) -
 111                         static_cast<int_x2_type>(point3.x());
 112       int_x2_type dy1 = static_cast<int_x2_type>(point1.y()) -
 113                         static_cast<int_x2_type>(point2.y());
 114       int_x2_type dy2 = static_cast<int_x2_type>(point2.y()) -
 115                         static_cast<int_x2_type>(point3.y());
 116       return eval(robust_cross_product(dx1, dy1, dx2, dy2));
 117     }
 118   } ot;
 119
 120   template <typename Point>
 121   class point_comparison_predicate {
 122    public:
 123     typedef Point point_type;
 124
 125     bool operator()(const point_type& lhs, const point_type& rhs) const {
 126       if (lhs.x() == rhs.x())
 127         return lhs.y() < rhs.y();
 128       return lhs.x() < rhs.x();
 129     }
 130   };
 131
 132   template <typename Site, typename Circle>
 133   class event_comparison_predicate {
 134    public:
 135     typedef Site site_type;
 136     typedef Circle circle_type;
 137
 138     bool operator()(const site_type& lhs, const site_type& rhs) const {
 139       if (lhs.x0() != rhs.x0())
 140         return lhs.x0() < rhs.x0();
 141       if (!lhs.is_segment()) {
 142         if (!rhs.is_segment())
 143           return lhs.y0() < rhs.y0();
 144         if (is_vertical(rhs))
 145           return lhs.y0() <= rhs.y0();
 146         return true;
 147       } else {
 148         if (is_vertical(rhs)) {
 149           if (is_vertical(lhs))
 150             return lhs.y0() < rhs.y0();
 151           return false;
 152         }
 153         if (is_vertical(lhs))
 154           return true;
 155         if (lhs.y0() != rhs.y0())
 156           return lhs.y0() < rhs.y0();
 157         return ot::eval(lhs.point1(), lhs.point0(), rhs.point1()) == ot::LEFT;
 158       }
 159     }
 160
 161     bool operator()(const site_type& lhs, const circle_type& rhs) const {
 162       typename ulp_cmp_type::Result xCmp =
 163           ulp_cmp(to_fpt(lhs.x0()), to_fpt(rhs.lower_x()), ULPS);
 164       return xCmp == ulp_cmp_type::LESS;
 165     }
 166
 167     bool operator()(const circle_type& lhs, const site_type& rhs) const {
 168       typename ulp_cmp_type::Result xCmp =
 169           ulp_cmp(to_fpt(lhs.lower_x()), to_fpt(rhs.x0()), ULPS);
 170       return xCmp == ulp_cmp_type::LESS;
 171     }
 172
 173     bool operator()(const circle_type& lhs, const circle_type& rhs) const {
 174       if (lhs.lower_x() != rhs.lower_x()) {
 175         return lhs.lower_x() < rhs.lower_x();
 176       }
 177       return lhs.y() < rhs.y();
 178     }
 179
 180    private:
 181     ulp_cmp_type ulp_cmp;
 182     to_fpt_converter to_fpt;
 183   };
 184
 185   template <typename Site>
 186   class distance_predicate {
 187    public:
 188     typedef Site site_type;
 189     typedef typename site_type::point_type point_type;
 190
 191     // Returns true if a horizontal line going through a new site intersects
 192     // right arc at first, else returns false. If horizontal line goes
 193     // through intersection point of the given two arcs returns false also.
 194     bool operator()(const site_type& left_site,
 195                     const site_type& right_site,
 196                     const point_type& new_point) const {
 197       if (!left_site.is_segment()) {
 198         if (!right_site.is_segment()) {
 199           return pp(left_site, right_site, new_point);
 200         } else {
 201           return ps(left_site, right_site, new_point, false);
 202         }
 203       } else {
 204         if (!right_site.is_segment()) {
 205           return ps(right_site, left_site, new_point, true);
 206         } else {
 207           return ss(left_site, right_site, new_point);
 208         }
 209       }
 210     }
 211
 212    private:
 213     // Represents the result of the epsilon robust predicate. If the
 214     // result is undefined some further processing is usually required.
 215     enum kPredicateResult {
 216       LESS = -1,
 217       UNDEFINED = 0,
 218       MORE = 1
 219     };
 220
 221     // Robust predicate, avoids using high-precision libraries.
 222     // Returns true if a horizontal line going through the new point site
 223     // intersects right arc at first, else returns false. If horizontal line
 224     // goes through intersection point of the given two arcs returns false.
 225     bool pp(const site_type& left_site,
 226             const site_type& right_site,
 227             const point_type& new_point) const {
 228       const point_type& left_point = left_site.point0();
 229       const point_type& right_point = right_site.point0();
 230       if (left_point.x() > right_point.x()) {
 231         if (new_point.y() <= left_point.y())
 232           return false;
 233       } else if (left_point.x() < right_point.x()) {
 234         if (new_point.y() >= right_point.y())
 235           return true;
 236       } else {
 237         return static_cast<int_x2_type>(left_point.y()) +
 238                static_cast<int_x2_type>(right_point.y()) <
 239                static_cast<int_x2_type>(new_point.y()) * 2;
 240       }
 241
 242       fpt_type dist1 = find_distance_to_point_arc(left_site, new_point);
 243       fpt_type dist2 = find_distance_to_point_arc(right_site, new_point);
 244
 245       // The undefined ulp range is equal to 3EPS + 3EPS <= 6ULP.
 246       return dist1 < dist2;
 247     }
 248
 249     bool ps(const site_type& left_site, const site_type& right_site,
 250             const point_type& new_point, bool reverse_order) const {
 251       kPredicateResult fast_res = fast_ps(
 252         left_site, right_site, new_point, reverse_order);
 253       if (fast_res != UNDEFINED) {
 254         return fast_res == LESS;
 255       }
 256
 257       fpt_type dist1 = find_distance_to_point_arc(left_site, new_point);
 258       fpt_type dist2 = find_distance_to_segment_arc(right_site, new_point);
 259
 260       // The undefined ulp range is equal to 3EPS + 7EPS <= 10ULP.
 261       return reverse_order ^ (dist1 < dist2);
 262     }
 263
 264     bool ss(const site_type& left_site,
 265             const site_type& right_site,
 266             const point_type& new_point) const {
 267       // Handle temporary segment sites.
 268       if (left_site.sorted_index() == right_site.sorted_index()) {
 269         return ot::eval(
 270             left_site.point0(), left_site.point1(), new_point) == ot::LEFT;
 271       }
 272
 273       fpt_type dist1 = find_distance_to_segment_arc(left_site, new_point);
 274       fpt_type dist2 = find_distance_to_segment_arc(right_site, new_point);
 275
 276       // The undefined ulp range is equal to 7EPS + 7EPS <= 14ULP.
 277       return dist1 < dist2;
 278     }
 279
 280     fpt_type find_distance_to_point_arc(
 281         const site_type& site, const point_type& point) const {
 282       fpt_type dx = to_fpt(site.x()) - to_fpt(point.x());
 283       fpt_type dy = to_fpt(site.y()) - to_fpt(point.y());
 284       // The relative error is at most 3EPS.
 285       return (dx * dx + dy * dy) / (to_fpt(2.0) * dx);
 286     }
 287
 288     fpt_type find_distance_to_segment_arc(
 289         const site_type& site, const point_type& point) const {
 290       if (is_vertical(site)) {
 291         return (to_fpt(site.x()) - to_fpt(point.x())) * to_fpt(0.5);
 292       } else {
 293         const point_type& segment0 = site.point0();
 294         const point_type& segment1 = site.point1();
 295         fpt_type a1 = to_fpt(segment1.x()) - to_fpt(segment0.x());
 296         fpt_type b1 = to_fpt(segment1.y()) - to_fpt(segment0.y());
 297         fpt_type k = get_sqrt(a1 * a1 + b1 * b1);
 298         // Avoid subtraction while computing k.
 299         if (!is_neg(b1)) {
 300           k = to_fpt(1.0) / (b1 + k);
 301         } else {
 302           k = (k - b1) / (a1 * a1);
 303         }
 304         // The relative error is at most 7EPS.
 305         return k * robust_cross_product(
 306             static_cast<int_x2_type>(segment1.x()) -
 307             static_cast<int_x2_type>(segment0.x()),
 308             static_cast<int_x2_type>(segment1.y()) -
 309             static_cast<int_x2_type>(segment0.y()),
 310             static_cast<int_x2_type>(point.x()) -
 311             static_cast<int_x2_type>(segment0.x()),
 312             static_cast<int_x2_type>(point.y()) -
 313             static_cast<int_x2_type>(segment0.y()));
 314       }
 315     }
 316
 317     kPredicateResult fast_ps(
 318         const site_type& left_site, const site_type& right_site,
 319         const point_type& new_point, bool reverse_order) const {
 320       const point_type& site_point = left_site.point0();
 321       const point_type& segment_start = right_site.point0();
 322       const point_type& segment_end = right_site.point1();
 323
 324       if (ot::eval(segment_start, segment_end, new_point) != ot::RIGHT)
 325         return (!right_site.is_inverse()) ? LESS : MORE;
 326
 327       fpt_type dif_x = to_fpt(new_point.x()) - to_fpt(site_point.x());
 328       fpt_type dif_y = to_fpt(new_point.y()) - to_fpt(site_point.y());
 329       fpt_type a = to_fpt(segment_end.x()) - to_fpt(segment_start.x());
 330       fpt_type b = to_fpt(segment_end.y()) - to_fpt(segment_start.y());
 331
 332       if (is_vertical(right_site)) {
 333         if (new_point.y() < site_point.y() && !reverse_order)
 334           return MORE;
 335         else if (new_point.y() > site_point.y() && reverse_order)
 336           return LESS;
 337         return UNDEFINED;
 338       } else {
 339         typename ot::Orientation orientation = ot::eval(
 340             static_cast<int_x2_type>(segment_end.x()) -
 341             static_cast<int_x2_type>(segment_start.x()),
 342             static_cast<int_x2_type>(segment_end.y()) -
 343             static_cast<int_x2_type>(segment_start.y()),
 344             static_cast<int_x2_type>(new_point.x()) -
 345             static_cast<int_x2_type>(site_point.x()),
 346             static_cast<int_x2_type>(new_point.y()) -
 347             static_cast<int_x2_type>(site_point.y()));
 348         if (orientation == ot::LEFT) {
 349           if (!right_site.is_inverse())
 350             return reverse_order ? LESS : UNDEFINED;
 351           return reverse_order ? UNDEFINED : MORE;
 352         }
 353       }
 354
 355       fpt_type fast_left_expr = a * (dif_y + dif_x) * (dif_y - dif_x);
 356       fpt_type fast_right_expr = (to_fpt(2.0) * b) * dif_x * dif_y;
 357       typename ulp_cmp_type::Result expr_cmp =
 358           ulp_cmp(fast_left_expr, fast_right_expr, 4);
 359       if (expr_cmp != ulp_cmp_type::EQUAL) {
 360         if ((expr_cmp == ulp_cmp_type::MORE) ^ reverse_order)
 361           return reverse_order ? LESS : MORE;
 362         return UNDEFINED;
 363       }
 364       return UNDEFINED;
 365     }
 366
 367    private:
 368     ulp_cmp_type ulp_cmp;
 369     to_fpt_converter to_fpt;
 370   };
 371
 372   template <typename Node>
 373   class node_comparison_predicate {
 374    public:
 375     typedef Node node_type;
 376     typedef typename Node::site_type site_type;
 377     typedef typename site_type::point_type point_type;
 378     typedef typename point_type::coordinate_type coordinate_type;
 379     typedef point_comparison_predicate<point_type> point_comparison_type;
 380     typedef distance_predicate<site_type> distance_predicate_type;
 381
 382     // Compares nodes in the balanced binary search tree. Nodes are
 383     // compared based on the y coordinates of the arcs intersection points.
 384     // Nodes with less y coordinate of the intersection point go first.
 385     // Comparison is only called during the new site events processing.
 386     // That's why one of the nodes will always lie on the sweepline and may
 387     // be represented as a straight horizontal line.
 388     bool operator() (const node_type& node1,
 389                      const node_type& node2) const {
 390       // Get x coordinate of the rightmost site from both nodes.
 391       const site_type& site1 = get_comparison_site(node1);
 392       const site_type& site2 = get_comparison_site(node2);
 393       const point_type& point1 = get_comparison_point(site1);
 394       const point_type& point2 = get_comparison_point(site2);
 395
 396       if (point1.x() < point2.x()) {
 397         // The second node contains a new site.
 398         return distance_predicate_(
 399             node1.left_site(), node1.right_site(), point2);
 400       } else if (point1.x() > point2.x()) {
 401         // The first node contains a new site.
 402         return !distance_predicate_(
 403             node2.left_site(), node2.right_site(), point1);
 404       } else {
 405         // This checks were evaluated experimentally.
 406         if (site1.sorted_index() == site2.sorted_index()) {
 407           // Both nodes are new (inserted during same site event processing).
 408           return get_comparison_y(node1) < get_comparison_y(node2);
 409         } else if (site1.sorted_index() < site2.sorted_index()) {
 410           std::pair<coordinate_type, int> y1 = get_comparison_y(node1, false);
 411           std::pair<coordinate_type, int> y2 = get_comparison_y(node2, true);
 412           if (y1.first != y2.first) return y1.first < y2.first;
 413           return (!site1.is_segment()) ? (y1.second < 0) : false;
 414         } else {
 415           std::pair<coordinate_type, int> y1 = get_comparison_y(node1, true);
 416           std::pair<coordinate_type, int> y2 = get_comparison_y(node2, false);
 417           if (y1.first != y2.first) return y1.first < y2.first;
 418           return (!site2.is_segment()) ? (y2.second > 0) : true;
 419         }
 420       }
 421     }
 422
 423    private:
 424     // Get the newer site.
 425     const site_type& get_comparison_site(const node_type& node) const {
 426       if (node.left_site().sorted_index() > node.right_site().sorted_index()) {
 427         return node.left_site();
 428       }
 429       return node.right_site();
 430     }
 431
 432     const point_type& get_comparison_point(const site_type& site) const {
 433       return point_comparison_(site.point0(), site.point1()) ?
 434           site.point0() : site.point1();
 435     }
 436
 437     // Get comparison pair: y coordinate and direction of the newer site.
 438     std::pair<coordinate_type, int> get_comparison_y(
 439       const node_type& node, bool is_new_node = true) const {
 440       if (node.left_site().sorted_index() ==
 441           node.right_site().sorted_index()) {
 442         return std::make_pair(node.left_site().y0(), 0);
 443       }
 444       if (node.left_site().sorted_index() > node.right_site().sorted_index()) {
 445         if (!is_new_node &&
 446             node.left_site().is_segment() &&
 447             is_vertical(node.left_site())) {
 448           return std::make_pair(node.left_site().y0(), 1);
 449         }
 450         return std::make_pair(node.left_site().y1(), 1);
 451       }
 452       return std::make_pair(node.right_site().y0(), -1);
 453     }
 454
 455     point_comparison_type point_comparison_;
 456     distance_predicate_type distance_predicate_;
 457   };
 458
 459   template <typename Site>
 460   class circle_existence_predicate {
 461    public:
 462     typedef typename Site::point_type point_type;
 463     typedef Site site_type;
 464
 465     bool ppp(const site_type& site1,
 466              const site_type& site2,
 467              const site_type& site3) const {
 468       return ot::eval(site1.point0(),
 469                       site2.point0(),
 470                       site3.point0()) == ot::RIGHT;
 471     }
 472
 473     bool pps(const site_type& site1,
 474              const site_type& site2,
 475              const site_type& site3,
 476              int segment_index) const {
 477       if (segment_index != 2) {
 478         typename ot::Orientation orient1 = ot::eval(
 479             site1.point0(), site2.point0(), site3.point0());
 480         typename ot::Orientation orient2 = ot::eval(
 481             site1.point0(), site2.point0(), site3.point1());
 482         if (segment_index == 1 && site1.x0() >= site2.x0()) {
 483           if (orient1 != ot::RIGHT)
 484             return false;
 485         } else if (segment_index == 3 && site2.x0() >= site1.x0()) {
 486           if (orient2 != ot::RIGHT)
 487             return false;
 488         } else if (orient1 != ot::RIGHT && orient2 != ot::RIGHT) {
 489           return false;
 490         }
 491       } else {
 492         return (site3.point0() != site1.point0()) ||
 493                (site3.point1() != site2.point0());
 494       }
 495       return true;
 496     }
 497
 498     bool pss(const site_type& site1,
 499              const site_type& site2,
 500              const site_type& site3,
 501              int point_index) const {
 502       if (site2.sorted_index() == site3.sorted_index()) {
 503         return false;
 504       }
 505       if (point_index == 2) {
 506         if (!site2.is_inverse() && site3.is_inverse())
 507           return false;
 508         if (site2.is_inverse() == site3.is_inverse() &&
 509             ot::eval(site2.point0(),
 510                      site1.point0(),
 511                      site3.point1()) != ot::RIGHT)
 512           return false;
 513       }
 514       return true;
 515     }
 516
 517     bool sss(const site_type& site1,
 518              const site_type& site2,
 519              const site_type& site3) const {
 520       return (site1.sorted_index() != site2.sorted_index()) &&
 521              (site2.sorted_index() != site3.sorted_index());
 522     }
 523   };
 524
 525   template <typename Site, typename Circle>
 526   class mp_circle_formation_functor {
 527    public:
 528     typedef typename Site::point_type point_type;
 529     typedef Site site_type;
 530     typedef Circle circle_type;
 531     typedef robust_sqrt_expr<big_int_type, efpt_type, to_efpt_converter>
 532         robust_sqrt_expr_type;
 533
 534     void ppp(const site_type& site1,
 535              const site_type& site2,
 536              const site_type& site3,
 537              circle_type& circle,
 538              bool recompute_c_x = true,
 539              bool recompute_c_y = true,
 540              bool recompute_lower_x = true) {
 541       big_int_type dif_x[3], dif_y[3], sum_x[2], sum_y[2];
 542       dif_x[0] = static_cast<int_x2_type>(site1.x()) -
 543                  static_cast<int_x2_type>(site2.x());
 544       dif_x[1] = static_cast<int_x2_type>(site2.x()) -
 545                  static_cast<int_x2_type>(site3.x());
 546       dif_x[2] = static_cast<int_x2_type>(site1.x()) -
 547                  static_cast<int_x2_type>(site3.x());
 548       dif_y[0] = static_cast<int_x2_type>(site1.y()) -
 549                  static_cast<int_x2_type>(site2.y());
 550       dif_y[1] = static_cast<int_x2_type>(site2.y()) -
 551                  static_cast<int_x2_type>(site3.y());
 552       dif_y[2] = static_cast<int_x2_type>(site1.y()) -
 553                  static_cast<int_x2_type>(site3.y());
 554       sum_x[0] = static_cast<int_x2_type>(site1.x()) +
 555                  static_cast<int_x2_type>(site2.x());
 556       sum_x[1] = static_cast<int_x2_type>(site2.x()) +
 557                  static_cast<int_x2_type>(site3.x());
 558       sum_y[0] = static_cast<int_x2_type>(site1.y()) +
 559                  static_cast<int_x2_type>(site2.y());
 560       sum_y[1] = static_cast<int_x2_type>(site2.y()) +
 561                  static_cast<int_x2_type>(site3.y());
 562       fpt_type inv_denom = to_fpt(0.5) / to_fpt(static_cast<big_int_type>(
 563           dif_x[0] * dif_y[1] - dif_x[1] * dif_y[0]));
 564       big_int_type numer1 = dif_x[0] * sum_x[0] + dif_y[0] * sum_y[0];
 565       big_int_type numer2 = dif_x[1] * sum_x[1] + dif_y[1] * sum_y[1];
 566
 567       if (recompute_c_x || recompute_lower_x) {
 568         big_int_type c_x = numer1 * dif_y[1] - numer2 * dif_y[0];
 569         circle.x(to_fpt(c_x) * inv_denom);
 570
 571         if (recompute_lower_x) {
 572           // Evaluate radius of the circle.
 573           big_int_type sqr_r = (dif_x[0] * dif_x[0] + dif_y[0] * dif_y[0]) *
 574                                (dif_x[1] * dif_x[1] + dif_y[1] * dif_y[1]) *
 575                                (dif_x[2] * dif_x[2] + dif_y[2] * dif_y[2]);
 576           fpt_type r = get_sqrt(to_fpt(sqr_r));
 577
 578           // If c_x >= 0 then lower_x = c_x + r,
 579           // else lower_x = (c_x * c_x - r * r) / (c_x - r).
 580           // To guarantee epsilon relative error.
 581           if (!is_neg(circle.x())) {
 582             if (!is_neg(inv_denom)) {
 583               circle.lower_x(circle.x() + r * inv_denom);
 584             } else {
 585               circle.lower_x(circle.x() - r * inv_denom);
 586             }
 587           } else {
 588             big_int_type numer = c_x * c_x - sqr_r;
 589             fpt_type lower_x = to_fpt(numer) * inv_denom / (to_fpt(c_x) + r);
 590             circle.lower_x(lower_x);
 591           }
 592         }
 593       }
 594
 595       if (recompute_c_y) {
 596         big_int_type c_y = numer2 * dif_x[0] - numer1 * dif_x[1];
 597         circle.y(to_fpt(c_y) * inv_denom);
 598       }
 599     }
 600
 601     // Recompute parameters of the circle event using high-precision library.
 602     void pps(const site_type& site1,
 603              const site_type& site2,
 604              const site_type& site3,
 605              int segment_index,
 606              circle_type& c_event,
 607              bool recompute_c_x = true,
 608              bool recompute_c_y = true,
 609              bool recompute_lower_x = true) {
 610       big_int_type cA[4], cB[4];
 611       big_int_type line_a = static_cast<int_x2_type>(site3.y1()) -
 612                             static_cast<int_x2_type>(site3.y0());
 613       big_int_type line_b = static_cast<int_x2_type>(site3.x0()) -
 614                             static_cast<int_x2_type>(site3.x1());
 615       big_int_type segm_len = line_a * line_a + line_b * line_b;
 616       big_int_type vec_x = static_cast<int_x2_type>(site2.y()) -
 617                            static_cast<int_x2_type>(site1.y());
 618       big_int_type vec_y = static_cast<int_x2_type>(site1.x()) -
 619                            static_cast<int_x2_type>(site2.x());
 620       big_int_type sum_x = static_cast<int_x2_type>(site1.x()) +
 621                            static_cast<int_x2_type>(site2.x());
 622       big_int_type sum_y = static_cast<int_x2_type>(site1.y()) +
 623                            static_cast<int_x2_type>(site2.y());
 624       big_int_type teta = line_a * vec_x + line_b * vec_y;
 625       big_int_type denom = vec_x * line_b - vec_y * line_a;
 626
 627       big_int_type dif0 = static_cast<int_x2_type>(site3.y1()) -
 628                           static_cast<int_x2_type>(site1.y());
 629       big_int_type dif1 = static_cast<int_x2_type>(site1.x()) -
 630                           static_cast<int_x2_type>(site3.x1());
 631       big_int_type A = line_a * dif1 - line_b * dif0;
 632       dif0 = static_cast<int_x2_type>(site3.y1()) -
 633              static_cast<int_x2_type>(site2.y());
 634       dif1 = static_cast<int_x2_type>(site2.x()) -
 635              static_cast<int_x2_type>(site3.x1());
 636       big_int_type B = line_a * dif1 - line_b * dif0;
 637       big_int_type sum_AB = A + B;
 638
 639       if (is_zero(denom)) {
 640         big_int_type numer = teta * teta - sum_AB * sum_AB;
 641         denom = teta * sum_AB;
 642         cA[0] = denom * sum_x * 2 + numer * vec_x;
 643         cB[0] = segm_len;
 644         cA[1] = denom * sum_AB * 2 + numer * teta;
 645         cB[1] = 1;
 646         cA[2] = denom * sum_y * 2 + numer * vec_y;
 647         fpt_type inv_denom = to_fpt(1.0) / to_fpt(denom);
 648         if (recompute_c_x)
 649           c_event.x(to_fpt(0.25) * to_fpt(cA[0]) * inv_denom);
 650         if (recompute_c_y)
 651           c_event.y(to_fpt(0.25) * to_fpt(cA[2]) * inv_denom);
 652         if (recompute_lower_x) {
 653           c_event.lower_x(to_fpt(0.25) * to_fpt(sqrt_expr_.eval2(cA, cB)) *
 654               inv_denom / get_sqrt(to_fpt(segm_len)));
 655         }
 656         return;
 657       }
 658
 659       big_int_type det = (teta * teta + denom * denom) * A * B * 4;
 660       fpt_type inv_denom_sqr = to_fpt(1.0) / to_fpt(denom);
 661       inv_denom_sqr *= inv_denom_sqr;
 662
 663       if (recompute_c_x || recompute_lower_x) {
 664         cA[0] = sum_x * denom * denom + teta * sum_AB * vec_x;
 665         cB[0] = 1;
 666         cA[1] = (segment_index == 2) ? -vec_x : vec_x;
 667         cB[1] = det;
 668         if (recompute_c_x) {
 669           c_event.x(to_fpt(0.5) * to_fpt(sqrt_expr_.eval2(cA, cB)) *
 670               inv_denom_sqr);
 671         }
 672       }
 673
 674       if (recompute_c_y || recompute_lower_x) {
 675         cA[2] = sum_y * denom * denom + teta * sum_AB * vec_y;
 676         cB[2] = 1;
 677         cA[3] = (segment_index == 2) ? -vec_y : vec_y;
 678         cB[3] = det;
 679         if (recompute_c_y) {
 680           c_event.y(to_fpt(0.5) * to_fpt(sqrt_expr_.eval2(&cA[2], &cB[2])) *
 681                     inv_denom_sqr);
 682         }
 683       }
 684
 685       if (recompute_lower_x) {
 686         cB[0] = cB[0] * segm_len;
 687         cB[1] = cB[1] * segm_len;
 688         cA[2] = sum_AB * (denom * denom + teta * teta);
 689         cB[2] = 1;
 690         cA[3] = (segment_index == 2) ? -teta : teta;
 691         cB[3] = det;
 692         c_event.lower_x(to_fpt(0.5) * to_fpt(sqrt_expr_.eval4(cA, cB)) *
 693             inv_denom_sqr / get_sqrt(to_fpt(segm_len)));
 694       }
 695     }
 696
 697     // Recompute parameters of the circle event using high-precision library.
 698     void pss(const site_type& site1,
 699              const site_type& site2,
 700              const site_type& site3,
 701              int point_index,
 702              circle_type& c_event,
 703              bool recompute_c_x = true,
 704              bool recompute_c_y = true,
 705              bool recompute_lower_x = true) {
 706       big_int_type a[2], b[2], c[2], cA[4], cB[4];
 707       const point_type& segm_start1 = site2.point1();
 708       const point_type& segm_end1 = site2.point0();
 709       const point_type& segm_start2 = site3.point0();
 710       const point_type& segm_end2 = site3.point1();
 711       a[0] = static_cast<int_x2_type>(segm_end1.x()) -
 712              static_cast<int_x2_type>(segm_start1.x());
 713       b[0] = static_cast<int_x2_type>(segm_end1.y()) -
 714              static_cast<int_x2_type>(segm_start1.y());
 715       a[1] = static_cast<int_x2_type>(segm_end2.x()) -
 716              static_cast<int_x2_type>(segm_start2.x());
 717       b[1] = static_cast<int_x2_type>(segm_end2.y()) -
 718              static_cast<int_x2_type>(segm_start2.y());
 719       big_int_type orientation = a[1] * b[0] - a[0] * b[1];
 720       if (is_zero(orientation)) {
 721         fpt_type denom = to_fpt(2.0) * to_fpt(
 722             static_cast<big_int_type>(a[0] * a[0] + b[0] * b[0]));
 723         c[0] = b[0] * (static_cast<int_x2_type>(segm_start2.x()) -
 724                        static_cast<int_x2_type>(segm_start1.x())) -
 725                a[0] * (static_cast<int_x2_type>(segm_start2.y()) -
 726                        static_cast<int_x2_type>(segm_start1.y()));
 727         big_int_type dx = a[0] * (static_cast<int_x2_type>(site1.y()) -
 728                                   static_cast<int_x2_type>(segm_start1.y())) -
 729                           b[0] * (static_cast<int_x2_type>(site1.x()) -
 730                                   static_cast<int_x2_type>(segm_start1.x()));
 731         big_int_type dy = b[0] * (static_cast<int_x2_type>(site1.x()) -
 732                                   static_cast<int_x2_type>(segm_start2.x())) -
 733                           a[0] * (static_cast<int_x2_type>(site1.y()) -
 734                                   static_cast<int_x2_type>(segm_start2.y()));
 735         cB[0] = dx * dy;
 736         cB[1] = 1;
 737
 738         if (recompute_c_y) {
 739           cA[0] = b[0] * ((point_index == 2) ? 2 : -2);
 740           cA[1] = a[0] * a[0] * (static_cast<int_x2_type>(segm_start1.y()) +
 741                                  static_cast<int_x2_type>(segm_start2.y())) -
 742                   a[0] * b[0] * (static_cast<int_x2_type>(segm_start1.x()) +
 743                                  static_cast<int_x2_type>(segm_start2.x()) -
 744                                  static_cast<int_x2_type>(site1.x()) * 2) +
 745                   b[0] * b[0] * (static_cast<int_x2_type>(site1.y()) * 2);
 746           fpt_type c_y = to_fpt(sqrt_expr_.eval2(cA, cB));
 747           c_event.y(c_y / denom);
 748         }
 749
 750         if (recompute_c_x || recompute_lower_x) {
 751           cA[0] = a[0] * ((point_index == 2) ? 2 : -2);
 752           cA[1] = b[0] * b[0] * (static_cast<int_x2_type>(segm_start1.x()) +
 753                                  static_cast<int_x2_type>(segm_start2.x())) -
 754                   a[0] * b[0] * (static_cast<int_x2_type>(segm_start1.y()) +
 755                                  static_cast<int_x2_type>(segm_start2.y()) -
 756                                  static_cast<int_x2_type>(site1.y()) * 2) +
 757                   a[0] * a[0] * (static_cast<int_x2_type>(site1.x()) * 2);
 758
 759           if (recompute_c_x) {
 760             fpt_type c_x = to_fpt(sqrt_expr_.eval2(cA, cB));
 761             c_event.x(c_x / denom);
 762           }
 763
 764           if (recompute_lower_x) {
 765             cA[2] = is_neg(c[0]) ? -c[0] : c[0];
 766             cB[2] = a[0] * a[0] + b[0] * b[0];
 767             fpt_type lower_x = to_fpt(sqrt_expr_.eval3(cA, cB));
 768             c_event.lower_x(lower_x / denom);
 769           }
 770         }
 771         return;
 772       }
 773       c[0] = b[0] * segm_end1.x() - a[0] * segm_end1.y();
 774       c[1] = a[1] * segm_end2.y() - b[1] * segm_end2.x();
 775       big_int_type ix = a[0] * c[1] + a[1] * c[0];
 776       big_int_type iy = b[0] * c[1] + b[1] * c[0];
 777       big_int_type dx = ix - orientation * site1.x();
 778       big_int_type dy = iy - orientation * site1.y();
 779       if (is_zero(dx) && is_zero(dy)) {
 780         fpt_type denom = to_fpt(orientation);
 781         fpt_type c_x = to_fpt(ix) / denom;
 782         fpt_type c_y = to_fpt(iy) / denom;
 783         c_event = circle_type(c_x, c_y, c_x);
 784         return;
 785       }
 786
 787       big_int_type sign = ((point_index == 2) ? 1 : -1) *
 788                           (is_neg(orientation) ? 1 : -1);
 789       cA[0] = a[1] * -dx + b[1] * -dy;
 790       cA[1] = a[0] * -dx + b[0] * -dy;
 791       cA[2] = sign;
 792       cA[3] = 0;
 793       cB[0] = a[0] * a[0] + b[0] * b[0];
 794       cB[1] = a[1] * a[1] + b[1] * b[1];
 795       cB[2] = a[0] * a[1] + b[0] * b[1];
 796       cB[3] = (a[0] * dy - b[0] * dx) * (a[1] * dy - b[1] * dx) * -2;
 797       fpt_type temp = to_fpt(
 798           sqrt_expr_evaluator_pss4<big_int_type, efpt_type>(cA, cB));
 799       fpt_type denom = temp * to_fpt(orientation);
 800
 801       if (recompute_c_y) {
 802         cA[0] = b[1] * (dx * dx + dy * dy) - iy * (dx * a[1] + dy * b[1]);
 803         cA[1] = b[0] * (dx * dx + dy * dy) - iy * (dx * a[0] + dy * b[0]);
 804         cA[2] = iy * sign;
 805         fpt_type cy = to_fpt(
 806             sqrt_expr_evaluator_pss4<big_int_type, efpt_type>(cA, cB));
 807         c_event.y(cy / denom);
 808       }
 809
 810       if (recompute_c_x || recompute_lower_x) {
 811         cA[0] = a[1] * (dx * dx + dy * dy) - ix * (dx * a[1] + dy * b[1]);
 812         cA[1] = a[0] * (dx * dx + dy * dy) - ix * (dx * a[0] + dy * b[0]);
 813         cA[2] = ix * sign;
 814
 815         if (recompute_c_x) {
 816           fpt_type cx = to_fpt(
 817               sqrt_expr_evaluator_pss4<big_int_type, efpt_type>(cA, cB));
 818           c_event.x(cx / denom);
 819         }
 820
 821         if (recompute_lower_x) {
 822           cA[3] = orientation * (dx * dx + dy * dy) * (is_neg(temp) ? -1 : 1);
 823           fpt_type lower_x = to_fpt(
 824               sqrt_expr_evaluator_pss4<big_int_type, efpt_type>(cA, cB));
 825           c_event.lower_x(lower_x / denom);
 826         }
 827       }
 828     }
 829
 830     // Recompute parameters of the circle event using high-precision library.
 831     void sss(const site_type& site1,
 832              const site_type& site2,
 833              const site_type& site3,
 834              circle_type& c_event,
 835              bool recompute_c_x = true,
 836              bool recompute_c_y = true,
 837              bool recompute_lower_x = true) {
 838       big_int_type a[3], b[3], c[3], cA[4], cB[4];
 839       // cA - corresponds to the cross product.
 840       // cB - corresponds to the squared length.
 841       a[0] = static_cast<int_x2_type>(site1.x1()) -
 842              static_cast<int_x2_type>(site1.x0());
 843       a[1] = static_cast<int_x2_type>(site2.x1()) -
 844              static_cast<int_x2_type>(site2.x0());
 845       a[2] = static_cast<int_x2_type>(site3.x1()) -
 846              static_cast<int_x2_type>(site3.x0());
 847
 848       b[0] = static_cast<int_x2_type>(site1.y1()) -
 849              static_cast<int_x2_type>(site1.y0());
 850       b[1] = static_cast<int_x2_type>(site2.y1()) -
 851              static_cast<int_x2_type>(site2.y0());
 852       b[2] = static_cast<int_x2_type>(site3.y1()) -
 853              static_cast<int_x2_type>(site3.y0());
 854
 855       c[0] = static_cast<int_x2_type>(site1.x0()) *
 856              static_cast<int_x2_type>(site1.y1()) -
 857              static_cast<int_x2_type>(site1.y0()) *
 858              static_cast<int_x2_type>(site1.x1());
 859       c[1] = static_cast<int_x2_type>(site2.x0()) *
 860              static_cast<int_x2_type>(site2.y1()) -
 861              static_cast<int_x2_type>(site2.y0()) *
 862              static_cast<int_x2_type>(site2.x1());
 863       c[2] = static_cast<int_x2_type>(site3.x0()) *
 864              static_cast<int_x2_type>(site3.y1()) -
 865              static_cast<int_x2_type>(site3.y0()) *
 866              static_cast<int_x2_type>(site3.x1());
 867
 868       for (int i = 0; i < 3; ++i)
 869         cB[i] = a[i] * a[i] + b[i] * b[i];
 870
 871       for (int i = 0; i < 3; ++i) {
 872         int j = (i+1) % 3;
 873         int k = (i+2) % 3;
 874         cA[i] = a[j] * b[k] - a[k] * b[j];
 875       }
 876       fpt_type denom = to_fpt(sqrt_expr_.eval3(cA, cB));
 877
 878       if (recompute_c_y) {
 879         for (int i = 0; i < 3; ++i) {
 880           int j = (i+1) % 3;
 881           int k = (i+2) % 3;
 882           cA[i] = b[j] * c[k] - b[k] * c[j];
 883         }
 884         fpt_type c_y = to_fpt(sqrt_expr_.eval3(cA, cB));
 885         c_event.y(c_y / denom);
 886       }
 887
 888       if (recompute_c_x || recompute_lower_x) {
 889         cA[3] = 0;
 890         for (int i = 0; i < 3; ++i) {
 891           int j = (i+1) % 3;
 892           int k = (i+2) % 3;
 893           cA[i] = a[j] * c[k] - a[k] * c[j];
 894           if (recompute_lower_x) {
 895             cA[3] = cA[3] + cA[i] * b[i];
 896           }
 897         }
 898
 899         if (recompute_c_x) {
 900           fpt_type c_x = to_fpt(sqrt_expr_.eval3(cA, cB));
 901           c_event.x(c_x / denom);
 902         }
 903
 904         if (recompute_lower_x) {
 905           cB[3] = 1;
 906           fpt_type lower_x = to_fpt(sqrt_expr_.eval4(cA, cB));
 907           c_event.lower_x(lower_x / denom);
 908         }
 909       }
 910     }
 911
 912    private:
 913     // Evaluates A[3] + A[0] * sqrt(B[0]) + A[1] * sqrt(B[1]) +
 914     //           A[2] * sqrt(B[3] * (sqrt(B[0] * B[1]) + B[2])).
 915     template <typename _int, typename _fpt>
 916     _fpt sqrt_expr_evaluator_pss4(_int *A, _int *B) {
 917       _int cA[4], cB[4];
 918       if (is_zero(A[3])) {
 919         _fpt lh = sqrt_expr_.eval2(A, B);
 920         cA[0] = 1;
 921         cB[0] = B[0] * B[1];
 922         cA[1] = B[2];
 923         cB[1] = 1;
 924         _fpt rh = sqrt_expr_.eval1(A+2, B+3) *
 925             get_sqrt(sqrt_expr_.eval2(cA, cB));
 926         if ((!is_neg(lh) && !is_neg(rh)) || (!is_pos(lh) && !is_pos(rh)))
 927           return lh + rh;
 928         cA[0] = A[0] * A[0] * B[0] + A[1] * A[1] * B[1] -
 929                 A[2] * A[2] * B[3] * B[2];
 930         cB[0] = 1;
 931         cA[1] = A[0] * A[1] * 2 - A[2] * A[2] * B[3];
 932         cB[1] = B[0] * B[1];
 933         _fpt numer = sqrt_expr_.eval2(cA, cB);
 934         return numer / (lh - rh);
 935       }
 936       cA[0] = 1;
 937       cB[0] = B[0] * B[1];
 938       cA[1] = B[2];
 939       cB[1] = 1;
 940       _fpt rh = sqrt_expr_.eval1(A+2, B+3) * get_sqrt(sqrt_expr_.eval2(cA, cB));
 941       cA[0] = A[0];
 942       cB[0] = B[0];
 943       cA[1] = A[1];
 944       cB[1] = B[1];
 945       cA[2] = A[3];
 946       cB[2] = 1;
 947       _fpt lh = sqrt_expr_.eval3(cA, cB);
 948       if ((!is_neg(lh) && !is_neg(rh)) || (!is_pos(lh) && !is_pos(rh)))
 949         return lh + rh;
 950       cA[0] = A[3] * A[0] * 2;
 951       cA[1] = A[3] * A[1] * 2;
 952       cA[2] = A[0] * A[0] * B[0] + A[1] * A[1] * B[1] +
 953               A[3] * A[3] - A[2] * A[2] * B[2] * B[3];
 954       cA[3] = A[0] * A[1] * 2 - A[2] * A[2] * B[3];
 955       cB[3] = B[0] * B[1];
 956       _fpt numer = sqrt_expr_evaluator_pss3<_int, _fpt>(cA, cB);
 957       return numer / (lh - rh);
 958     }
 959
 960     template <typename _int, typename _fpt>
 961     // Evaluates A[0] * sqrt(B[0]) + A[1] * sqrt(B[1]) +
 962     //           A[2] + A[3] * sqrt(B[0] * B[1]).
 963     // B[3] = B[0] * B[1].
 964     _fpt sqrt_expr_evaluator_pss3(_int *A, _int *B) {
 965       _int cA[2], cB[2];
 966       _fpt lh = sqrt_expr_.eval2(A, B);
 967       _fpt rh = sqrt_expr_.eval2(A+2, B+2);
 968       if ((!is_neg(lh) && !is_neg(rh)) || (!is_pos(lh) && !is_pos(rh)))
 969         return lh + rh;
 970       cA[0] = A[0] * A[0] * B[0] + A[1] * A[1] * B[1] -
 971               A[2] * A[2] - A[3] * A[3] * B[0] * B[1];
 972       cB[0] = 1;
 973       cA[1] = (A[0] * A[1] - A[2] * A[3]) * 2;
 974       cB[1] = B[3];
 975       _fpt numer = sqrt_expr_.eval2(cA, cB);
 976       return numer / (lh - rh);
 977     }
 978
 979     robust_sqrt_expr_type sqrt_expr_;
 980     to_fpt_converter to_fpt;
 981   };
 982
 983   template <typename Site, typename Circle>
 984   class lazy_circle_formation_functor {
 985    public:
 986     typedef robust_fpt<fpt_type> robust_fpt_type;
 987     typedef robust_dif<robust_fpt_type> robust_dif_type;
 988     typedef typename Site::point_type point_type;
 989     typedef Site site_type;
 990     typedef Circle circle_type;
 991     typedef mp_circle_formation_functor<site_type, circle_type>
 992         exact_circle_formation_functor_type;
 993
 994     void ppp(const site_type& site1,
 995              const site_type& site2,
 996              const site_type& site3,
 997              circle_type& c_event) {
 998       fpt_type dif_x1 = to_fpt(site1.x()) - to_fpt(site2.x());
 999       fpt_type dif_x2 = to_fpt(site2.x()) - to_fpt(site3.x());
1000       fpt_type dif_y1 = to_fpt(site1.y()) - to_fpt(site2.y());
1001       fpt_type dif_y2 = to_fpt(site2.y()) - to_fpt(site3.y());
1002       fpt_type orientation = robust_cross_product(
1003           static_cast<int_x2_type>(site1.x()) -
1004           static_cast<int_x2_type>(site2.x()),
1005           static_cast<int_x2_type>(site2.x()) -
1006           static_cast<int_x2_type>(site3.x()),
1007           static_cast<int_x2_type>(site1.y()) -
1008           static_cast<int_x2_type>(site2.y()),
1009           static_cast<int_x2_type>(site2.y()) -
1010           static_cast<int_x2_type>(site3.y()));
1011       robust_fpt_type inv_orientation(to_fpt(0.5) / orientation, to_fpt(2.0));
1012       fpt_type sum_x1 = to_fpt(site1.x()) + to_fpt(site2.x());
1013       fpt_type sum_x2 = to_fpt(site2.x()) + to_fpt(site3.x());
1014       fpt_type sum_y1 = to_fpt(site1.y()) + to_fpt(site2.y());
1015       fpt_type sum_y2 = to_fpt(site2.y()) + to_fpt(site3.y());
1016       fpt_type dif_x3 = to_fpt(site1.x()) - to_fpt(site3.x());
1017       fpt_type dif_y3 = to_fpt(site1.y()) - to_fpt(site3.y());
1018       robust_dif_type c_x, c_y;
1019       c_x += robust_fpt_type(dif_x1 * sum_x1 * dif_y2, to_fpt(2.0));
1020       c_x += robust_fpt_type(dif_y1 * sum_y1 * dif_y2, to_fpt(2.0));
1021       c_x -= robust_fpt_type(dif_x2 * sum_x2 * dif_y1, to_fpt(2.0));
1022       c_x -= robust_fpt_type(dif_y2 * sum_y2 * dif_y1, to_fpt(2.0));
1023       c_y += robust_fpt_type(dif_x2 * sum_x2 * dif_x1, to_fpt(2.0));
1024       c_y += robust_fpt_type(dif_y2 * sum_y2 * dif_x1, to_fpt(2.0));
1025       c_y -= robust_fpt_type(dif_x1 * sum_x1 * dif_x2, to_fpt(2.0));
1026       c_y -= robust_fpt_type(dif_y1 * sum_y1 * dif_x2, to_fpt(2.0));
1027       robust_dif_type lower_x(c_x);
1028       lower_x -= robust_fpt_type(get_sqrt(
1029           (dif_x1 * dif_x1 + dif_y1 * dif_y1) *
1030           (dif_x2 * dif_x2 + dif_y2 * dif_y2) *
1031           (dif_x3 * dif_x3 + dif_y3 * dif_y3)), to_fpt(5.0));
1032       c_event = circle_type(
1033           c_x.dif().fpv() * inv_orientation.fpv(),
1034           c_y.dif().fpv() * inv_orientation.fpv(),
1035           lower_x.dif().fpv() * inv_orientation.fpv());
1036       bool recompute_c_x = c_x.dif().ulp() > ULPS;
1037       bool recompute_c_y = c_y.dif().ulp() > ULPS;
1038       bool recompute_lower_x = lower_x.dif().ulp() > ULPS;
1039       if (recompute_c_x || recompute_c_y || recompute_lower_x) {
1040         exact_circle_formation_functor_.ppp(
1041             site1, site2, site3, c_event,
1042             recompute_c_x, recompute_c_y, recompute_lower_x);
1043       }
1044     }
1045
1046     void pps(const site_type& site1,
1047              const site_type& site2,
1048              const site_type& site3,
1049              int segment_index,
1050              circle_type& c_event) {
1051       fpt_type line_a = to_fpt(site3.y1()) - to_fpt(site3.y0());
1052       fpt_type line_b = to_fpt(site3.x0()) - to_fpt(site3.x1());
1053       fpt_type vec_x = to_fpt(site2.y()) - to_fpt(site1.y());
1054       fpt_type vec_y = to_fpt(site1.x()) - to_fpt(site2.x());
1055       robust_fpt_type teta(robust_cross_product(
1056           static_cast<int_x2_type>(site3.y1()) -
1057           static_cast<int_x2_type>(site3.y0()),
1058           static_cast<int_x2_type>(site3.x0()) -
1059           static_cast<int_x2_type>(site3.x1()),
1060           static_cast<int_x2_type>(site2.x()) -
1061           static_cast<int_x2_type>(site1.x()),
1062           static_cast<int_x2_type>(site2.y()) -
1063           static_cast<int_x2_type>(site1.y())), to_fpt(1.0));
1064       robust_fpt_type A(robust_cross_product(
1065           static_cast<int_x2_type>(site3.y0()) -
1066           static_cast<int_x2_type>(site3.y1()),
1067           static_cast<int_x2_type>(site3.x0()) -
1068           static_cast<int_x2_type>(site3.x1()),
1069           static_cast<int_x2_type>(site3.y1()) -
1070           static_cast<int_x2_type>(site1.y()),
1071           static_cast<int_x2_type>(site3.x1()) -
1072           static_cast<int_x2_type>(site1.x())), to_fpt(1.0));
1073       robust_fpt_type B(robust_cross_product(
1074           static_cast<int_x2_type>(site3.y0()) -
1075           static_cast<int_x2_type>(site3.y1()),
1076           static_cast<int_x2_type>(site3.x0()) -
1077           static_cast<int_x2_type>(site3.x1()),
1078           static_cast<int_x2_type>(site3.y1()) -
1079           static_cast<int_x2_type>(site2.y()),
1080           static_cast<int_x2_type>(site3.x1()) -
1081           static_cast<int_x2_type>(site2.x())), to_fpt(1.0));
1082       robust_fpt_type denom(robust_cross_product(
1083           static_cast<int_x2_type>(site1.y()) -
1084           static_cast<int_x2_type>(site2.y()),
1085           static_cast<int_x2_type>(site1.x()) -
1086           static_cast<int_x2_type>(site2.x()),
1087           static_cast<int_x2_type>(site3.y1()) -
1088           static_cast<int_x2_type>(site3.y0()),
1089           static_cast<int_x2_type>(site3.x1()) -
1090           static_cast<int_x2_type>(site3.x0())), to_fpt(1.0));
1091       robust_fpt_type inv_segm_len(to_fpt(1.0) /
1092           get_sqrt(line_a * line_a + line_b * line_b), to_fpt(3.0));
1093       robust_dif_type t;
1094       if (ot::eval(denom) == ot::COLLINEAR) {
1095         t += teta / (robust_fpt_type(to_fpt(8.0)) * A);
1096         t -= A / (robust_fpt_type(to_fpt(2.0)) * teta);
1097       } else {
1098         robust_fpt_type det = ((teta * teta + denom * denom) * A * B).sqrt();
1099         if (segment_index == 2) {
1100           t -= det / (denom * denom);
1101         } else {
1102           t += det / (denom * denom);
1103         }
1104         t += teta * (A + B) / (robust_fpt_type(to_fpt(2.0)) * denom * denom);
1105       }
1106       robust_dif_type c_x, c_y;
1107       c_x += robust_fpt_type(to_fpt(0.5) *
1108           (to_fpt(site1.x()) + to_fpt(site2.x())));
1109       c_x += robust_fpt_type(vec_x) * t;
1110       c_y += robust_fpt_type(to_fpt(0.5) *
1111           (to_fpt(site1.y()) + to_fpt(site2.y())));
1112       c_y += robust_fpt_type(vec_y) * t;
1113       robust_dif_type r, lower_x(c_x);
1114       r -= robust_fpt_type(line_a) * robust_fpt_type(site3.x0());
1115       r -= robust_fpt_type(line_b) * robust_fpt_type(site3.y0());
1116       r += robust_fpt_type(line_a) * c_x;
1117       r += robust_fpt_type(line_b) * c_y;
1118       if (r.pos().fpv() < r.neg().fpv())
1119         r = -r;
1120       lower_x += r * inv_segm_len;
1121       c_event = circle_type(
1122           c_x.dif().fpv(), c_y.dif().fpv(), lower_x.dif().fpv());
1123       bool recompute_c_x = c_x.dif().ulp() > ULPS;
1124       bool recompute_c_y = c_y.dif().ulp() > ULPS;
1125       bool recompute_lower_x = lower_x.dif().ulp() > ULPS;
1126       if (recompute_c_x || recompute_c_y || recompute_lower_x) {
1127         exact_circle_formation_functor_.pps(
1128             site1, site2, site3, segment_index, c_event,
1129             recompute_c_x, recompute_c_y, recompute_lower_x);
1130       }
1131     }
1132
1133     void pss(const site_type& site1,
1134              const site_type& site2,
1135              const site_type& site3,
1136              int point_index,
1137              circle_type& c_event) {
1138       const point_type& segm_start1 = site2.point1();
1139       const point_type& segm_end1 = site2.point0();
1140       const point_type& segm_start2 = site3.point0();
1141       const point_type& segm_end2 = site3.point1();
1142       fpt_type a1 = to_fpt(segm_end1.x()) - to_fpt(segm_start1.x());
1143       fpt_type b1 = to_fpt(segm_end1.y()) - to_fpt(segm_start1.y());
1144       fpt_type a2 = to_fpt(segm_end2.x()) - to_fpt(segm_start2.x());
1145       fpt_type b2 = to_fpt(segm_end2.y()) - to_fpt(segm_start2.y());
1146       bool recompute_c_x, recompute_c_y, recompute_lower_x;
1147       robust_fpt_type orientation(robust_cross_product(
1148         static_cast<int_x2_type>(segm_end1.y()) -
1149         static_cast<int_x2_type>(segm_start1.y()),
1150         static_cast<int_x2_type>(segm_end1.x()) -
1151         static_cast<int_x2_type>(segm_start1.x()),
1152         static_cast<int_x2_type>(segm_end2.y()) -
1153         static_cast<int_x2_type>(segm_start2.y()),
1154         static_cast<int_x2_type>(segm_end2.x()) -
1155         static_cast<int_x2_type>(segm_start2.x())), to_fpt(1.0));
1156       if (ot::eval(orientation) == ot::COLLINEAR) {
1157         robust_fpt_type a(a1 * a1 + b1 * b1, to_fpt(2.0));
1158         robust_fpt_type c(robust_cross_product(
1159             static_cast<int_x2_type>(segm_end1.y()) -
1160             static_cast<int_x2_type>(segm_start1.y()),
1161             static_cast<int_x2_type>(segm_end1.x()) -
1162             static_cast<int_x2_type>(segm_start1.x()),
1163             static_cast<int_x2_type>(segm_start2.y()) -
1164             static_cast<int_x2_type>(segm_start1.y()),
1165             static_cast<int_x2_type>(segm_start2.x()) -
1166             static_cast<int_x2_type>(segm_start1.x())), to_fpt(1.0));
1167         robust_fpt_type det(
1168             robust_cross_product(
1169                 static_cast<int_x2_type>(segm_end1.x()) -
1170                 static_cast<int_x2_type>(segm_start1.x()),
1171                 static_cast<int_x2_type>(segm_end1.y()) -
1172                 static_cast<int_x2_type>(segm_start1.y()),
1173                 static_cast<int_x2_type>(site1.x()) -
1174                 static_cast<int_x2_type>(segm_start1.x()),
1175                 static_cast<int_x2_type>(site1.y()) -
1176                 static_cast<int_x2_type>(segm_start1.y())) *
1177             robust_cross_product(
1178                 static_cast<int_x2_type>(segm_end1.y()) -
1179                 static_cast<int_x2_type>(segm_start1.y()),
1180                 static_cast<int_x2_type>(segm_end1.x()) -
1181                 static_cast<int_x2_type>(segm_start1.x()),
1182                 static_cast<int_x2_type>(site1.y()) -
1183                 static_cast<int_x2_type>(segm_start2.y()),
1184                 static_cast<int_x2_type>(site1.x()) -
1185                 static_cast<int_x2_type>(segm_start2.x())),
1186             to_fpt(3.0));
1187         robust_dif_type t;
1188         t -= robust_fpt_type(a1) * robust_fpt_type((
1189              to_fpt(segm_start1.x()) + to_fpt(segm_start2.x())) * to_fpt(0.5) -
1190              to_fpt(site1.x()));
1191         t -= robust_fpt_type(b1) * robust_fpt_type((
1192              to_fpt(segm_start1.y()) + to_fpt(segm_start2.y())) * to_fpt(0.5) -
1193              to_fpt(site1.y()));
1194         if (point_index == 2) {
1195           t += det.sqrt();
1196         } else {
1197           t -= det.sqrt();
1198         }
1199         t /= a;
1200         robust_dif_type c_x, c_y;
1201         c_x += robust_fpt_type(to_fpt(0.5) * (
1202             to_fpt(segm_start1.x()) + to_fpt(segm_start2.x())));
1203         c_x += robust_fpt_type(a1) * t;
1204         c_y += robust_fpt_type(to_fpt(0.5) * (
1205             to_fpt(segm_start1.y()) + to_fpt(segm_start2.y())));
1206         c_y += robust_fpt_type(b1) * t;
1207         robust_dif_type lower_x(c_x);
1208         if (is_neg(c)) {
1209           lower_x -= robust_fpt_type(to_fpt(0.5)) * c / a.sqrt();
1210         } else {
1211           lower_x += robust_fpt_type(to_fpt(0.5)) * c / a.sqrt();
1212         }
1213         recompute_c_x = c_x.dif().ulp() > ULPS;
1214         recompute_c_y = c_y.dif().ulp() > ULPS;
1215         recompute_lower_x = lower_x.dif().ulp() > ULPS;
1216         c_event =
1217             circle_type(c_x.dif().fpv(), c_y.dif().fpv(), lower_x.dif().fpv());
1218       } else {
1219         robust_fpt_type sqr_sum1(get_sqrt(a1 * a1 + b1 * b1), to_fpt(2.0));
1220         robust_fpt_type sqr_sum2(get_sqrt(a2 * a2 + b2 * b2), to_fpt(2.0));
1221         robust_fpt_type a(robust_cross_product(
1222           static_cast<int_x2_type>(segm_end1.x()) -
1223           static_cast<int_x2_type>(segm_start1.x()),
1224           static_cast<int_x2_type>(segm_end1.y()) -
1225           static_cast<int_x2_type>(segm_start1.y()),
1226           static_cast<int_x2_type>(segm_start2.y()) -
1227           static_cast<int_x2_type>(segm_end2.y()),
1228           static_cast<int_x2_type>(segm_end2.x()) -
1229           static_cast<int_x2_type>(segm_start2.x())), to_fpt(1.0));
1230         if (!is_neg(a)) {
1231           a += sqr_sum1 * sqr_sum2;
1232         } else {
1233           a = (orientation * orientation) / (sqr_sum1 * sqr_sum2 - a);
1234         }
1235         robust_fpt_type or1(robust_cross_product(
1236             static_cast<int_x2_type>(segm_end1.y()) -
1237             static_cast<int_x2_type>(segm_start1.y()),
1238             static_cast<int_x2_type>(segm_end1.x()) -
1239             static_cast<int_x2_type>(segm_start1.x()),
1240             static_cast<int_x2_type>(segm_end1.y()) -
1241             static_cast<int_x2_type>(site1.y()),
1242             static_cast<int_x2_type>(segm_end1.x()) -
1243             static_cast<int_x2_type>(site1.x())), to_fpt(1.0));
1244         robust_fpt_type or2(robust_cross_product(
1245             static_cast<int_x2_type>(segm_end2.x()) -
1246             static_cast<int_x2_type>(segm_start2.x()),
1247             static_cast<int_x2_type>(segm_end2.y()) -
1248             static_cast<int_x2_type>(segm_start2.y()),
1249             static_cast<int_x2_type>(segm_end2.x()) -
1250             static_cast<int_x2_type>(site1.x()),
1251             static_cast<int_x2_type>(segm_end2.y()) -
1252             static_cast<int_x2_type>(site1.y())), to_fpt(1.0));
1253         robust_fpt_type det = robust_fpt_type(to_fpt(2.0)) * a * or1 * or2;
1254         robust_fpt_type c1(robust_cross_product(
1255             static_cast<int_x2_type>(segm_end1.y()) -
1256             static_cast<int_x2_type>(segm_start1.y()),
1257             static_cast<int_x2_type>(segm_end1.x()) -
1258             static_cast<int_x2_type>(segm_start1.x()),
1259             static_cast<int_x2_type>(segm_end1.y()),
1260             static_cast<int_x2_type>(segm_end1.x())), to_fpt(1.0));
1261         robust_fpt_type c2(robust_cross_product(
1262             static_cast<int_x2_type>(segm_end2.x()) -
1263             static_cast<int_x2_type>(segm_start2.x()),
1264             static_cast<int_x2_type>(segm_end2.y()) -
1265             static_cast<int_x2_type>(segm_start2.y()),
1266             static_cast<int_x2_type>(segm_end2.x()),
1267             static_cast<int_x2_type>(segm_end2.y())), to_fpt(1.0));
1268         robust_fpt_type inv_orientation =
1269             robust_fpt_type(to_fpt(1.0)) / orientation;
1270         robust_dif_type t, b, ix, iy;
1271         ix += robust_fpt_type(a2) * c1 * inv_orientation;
1272         ix += robust_fpt_type(a1) * c2 * inv_orientation;
1273         iy += robust_fpt_type(b1) * c2 * inv_orientation;
1274         iy += robust_fpt_type(b2) * c1 * inv_orientation;
1275
1276         b += ix * (robust_fpt_type(a1) * sqr_sum2);
1277         b += ix * (robust_fpt_type(a2) * sqr_sum1);
1278         b += iy * (robust_fpt_type(b1) * sqr_sum2);
1279         b += iy * (robust_fpt_type(b2) * sqr_sum1);
1280         b -= sqr_sum1 * robust_fpt_type(robust_cross_product(
1281             static_cast<int_x2_type>(segm_end2.x()) -
1282             static_cast<int_x2_type>(segm_start2.x()),
1283             static_cast<int_x2_type>(segm_end2.y()) -
1284             static_cast<int_x2_type>(segm_start2.y()),
1285             static_cast<int_x2_type>(-site1.y()),
1286             static_cast<int_x2_type>(site1.x())), to_fpt(1.0));
1287         b -= sqr_sum2 * robust_fpt_type(robust_cross_product(
1288             static_cast<int_x2_type>(segm_end1.x()) -
1289             static_cast<int_x2_type>(segm_start1.x()),
1290             static_cast<int_x2_type>(segm_end1.y()) -
1291             static_cast<int_x2_type>(segm_start1.y()),
1292             static_cast<int_x2_type>(-site1.y()),
1293             static_cast<int_x2_type>(site1.x())), to_fpt(1.0));
1294         t -= b;
1295         if (point_index == 2) {
1296           t += det.sqrt();
1297         } else {
1298           t -= det.sqrt();
1299         }
1300         t /= (a * a);
1301         robust_dif_type c_x(ix), c_y(iy);
1302         c_x += t * (robust_fpt_type(a1) * sqr_sum2);
1303         c_x += t * (robust_fpt_type(a2) * sqr_sum1);
1304         c_y += t * (robust_fpt_type(b1) * sqr_sum2);
1305         c_y += t * (robust_fpt_type(b2) * sqr_sum1);
1306         if (t.pos().fpv() < t.neg().fpv()) {
1307           t = -t;
1308         }
1309         robust_dif_type lower_x(c_x);
1310         if (is_neg(orientation)) {
1311           lower_x -= t * orientation;
1312         } else {
1313           lower_x += t * orientation;
1314         }
1315         recompute_c_x = c_x.dif().ulp() > ULPS;
1316         recompute_c_y = c_y.dif().ulp() > ULPS;
1317         recompute_lower_x = lower_x.dif().ulp() > ULPS;
1318         c_event = circle_type(
1319             c_x.dif().fpv(), c_y.dif().fpv(), lower_x.dif().fpv());
1320       }
1321       if (recompute_c_x || recompute_c_y || recompute_lower_x) {
1322           exact_circle_formation_functor_.pss(
1323               site1, site2, site3, point_index, c_event,
1324         recompute_c_x, recompute_c_y, recompute_lower_x);
1325       }
1326     }
1327
1328     void sss(const site_type& site1,
1329              const site_type& site2,
1330              const site_type& site3,
1331              circle_type& c_event) {
1332       robust_fpt_type a1(to_fpt(site1.x1()) - to_fpt(site1.x0()));
1333       robust_fpt_type b1(to_fpt(site1.y1()) - to_fpt(site1.y0()));
1334       robust_fpt_type c1(robust_cross_product(
1335           site1.x0(), site1.y0(),
1336           site1.x1(), site1.y1()), to_fpt(1.0));
1337
1338       robust_fpt_type a2(to_fpt(site2.x1()) - to_fpt(site2.x0()));
1339       robust_fpt_type b2(to_fpt(site2.y1()) - to_fpt(site2.y0()));
1340       robust_fpt_type c2(robust_cross_product(
1341           site2.x0(), site2.y0(),
1342           site2.x1(), site2.y1()), to_fpt(1.0));
1343
1344       robust_fpt_type a3(to_fpt(site3.x1()) - to_fpt(site3.x0()));
1345       robust_fpt_type b3(to_fpt(site3.y1()) - to_fpt(site3.y0()));
1346       robust_fpt_type c3(robust_cross_product(
1347           site3.x0(), site3.y0(),
1348           site3.x1(), site3.y1()), to_fpt(1.0));
1349
1350       robust_fpt_type len1 = (a1 * a1 + b1 * b1).sqrt();
1351       robust_fpt_type len2 = (a2 * a2 + b2 * b2).sqrt();
1352       robust_fpt_type len3 = (a3 * a3 + b3 * b3).sqrt();
1353       robust_fpt_type cross_12(robust_cross_product(
1354           static_cast<int_x2_type>(site1.x1()) -
1355           static_cast<int_x2_type>(site1.x0()),
1356           static_cast<int_x2_type>(site1.y1()) -
1357           static_cast<int_x2_type>(site1.y0()),
1358           static_cast<int_x2_type>(site2.x1()) -
1359           static_cast<int_x2_type>(site2.x0()),
1360           static_cast<int_x2_type>(site2.y1()) -
1361           static_cast<int_x2_type>(site2.y0())), to_fpt(1.0));
1362       robust_fpt_type cross_23(robust_cross_product(
1363           static_cast<int_x2_type>(site2.x1()) -
1364           static_cast<int_x2_type>(site2.x0()),
1365           static_cast<int_x2_type>(site2.y1()) -
1366           static_cast<int_x2_type>(site2.y0()),
1367           static_cast<int_x2_type>(site3.x1()) -
1368           static_cast<int_x2_type>(site3.x0()),
1369           static_cast<int_x2_type>(site3.y1()) -
1370           static_cast<int_x2_type>(site3.y0())), to_fpt(1.0));
1371       robust_fpt_type cross_31(robust_cross_product(
1372           static_cast<int_x2_type>(site3.x1()) -
1373           static_cast<int_x2_type>(site3.x0()),
1374           static_cast<int_x2_type>(site3.y1()) -
1375           static_cast<int_x2_type>(site3.y0()),
1376           static_cast<int_x2_type>(site1.x1()) -
1377           static_cast<int_x2_type>(site1.x0()),
1378           static_cast<int_x2_type>(site1.y1()) -
1379           static_cast<int_x2_type>(site1.y0())), to_fpt(1.0));
1380
1381       // denom = cross_12 * len3 + cross_23 * len1 + cross_31 * len2.
1382       robust_dif_type denom;
1383       denom += cross_12 * len3;
1384       denom += cross_23 * len1;
1385       denom += cross_31 * len2;
1386
1387       // denom * r = (b2 * c_x - a2 * c_y - c2 * denom) / len2.
1388       robust_dif_type r;
1389       r -= cross_12 * c3;
1390       r -= cross_23 * c1;
1391       r -= cross_31 * c2;
1392
1393       robust_dif_type c_x;
1394       c_x += a1 * c2 * len3;
1395       c_x -= a2 * c1 * len3;
1396       c_x += a2 * c3 * len1;
1397       c_x -= a3 * c2 * len1;
1398       c_x += a3 * c1 * len2;
1399       c_x -= a1 * c3 * len2;
1400
1401       robust_dif_type c_y;
1402       c_y += b1 * c2 * len3;
1403       c_y -= b2 * c1 * len3;
1404       c_y += b2 * c3 * len1;
1405       c_y -= b3 * c2 * len1;
1406       c_y += b3 * c1 * len2;
1407       c_y -= b1 * c3 * len2;
1408
1409       robust_dif_type lower_x = c_x + r;
1410
1411       robust_fpt_type denom_dif = denom.dif();
1412       robust_fpt_type c_x_dif = c_x.dif() / denom_dif;
1413       robust_fpt_type c_y_dif = c_y.dif() / denom_dif;
1414       robust_fpt_type lower_x_dif = lower_x.dif() / denom_dif;
1415
1416       bool recompute_c_x = c_x_dif.ulp() > ULPS;
1417       bool recompute_c_y = c_y_dif.ulp() > ULPS;
1418       bool recompute_lower_x = lower_x_dif.ulp() > ULPS;
1419       c_event = circle_type(c_x_dif.fpv(), c_y_dif.fpv(), lower_x_dif.fpv());
1420       if (recompute_c_x || recompute_c_y || recompute_lower_x) {
1421         exact_circle_formation_functor_.sss(
1422             site1, site2, site3, c_event,
1423             recompute_c_x, recompute_c_y, recompute_lower_x);
1424       }
1425     }
1426
1427    private:
1428     exact_circle_formation_functor_type exact_circle_formation_functor_;
1429     to_fpt_converter to_fpt;
1430   };
1431
1432   template <typename Site,
1433             typename Circle,
1434             typename CEP = circle_existence_predicate<Site>,
1435             typename CFF = lazy_circle_formation_functor<Site, Circle> >
1436   class circle_formation_predicate {
1437    public:
1438     typedef Site site_type;
1439     typedef Circle circle_type;
1440     typedef CEP circle_existence_predicate_type;
1441     typedef CFF circle_formation_functor_type;
1442
1443     // Create a circle event from the given three sites.
1444     // Returns true if the circle event exists, else false.
1445     // If exists circle event is saved into the c_event variable.
1446     bool operator()(const site_type& site1, const site_type& site2,
1447                     const site_type& site3, circle_type& circle) {
1448       if (!site1.is_segment()) {
1449         if (!site2.is_segment()) {
1450           if (!site3.is_segment()) {
1451             // (point, point, point) sites.
1452             if (!circle_existence_predicate_.ppp(site1, site2, site3))
1453               return false;
1454             circle_formation_functor_.ppp(site1, site2, site3, circle);
1455           } else {
1456             // (point, point, segment) sites.
1457             if (!circle_existence_predicate_.pps(site1, site2, site3, 3))
1458               return false;
1459             circle_formation_functor_.pps(site1, site2, site3, 3, circle);
1460           }
1461         } else {
1462           if (!site3.is_segment()) {
1463             // (point, segment, point) sites.
1464             if (!circle_existence_predicate_.pps(site1, site3, site2, 2))
1465               return false;
1466             circle_formation_functor_.pps(site1, site3, site2, 2, circle);
1467           } else {
1468             // (point, segment, segment) sites.
1469             if (!circle_existence_predicate_.pss(site1, site2, site3, 1))
1470               return false;
1471             circle_formation_functor_.pss(site1, site2, site3, 1, circle);
1472           }
1473         }
1474       } else {
1475         if (!site2.is_segment()) {
1476           if (!site3.is_segment()) {
1477             // (segment, point, point) sites.
1478             if (!circle_existence_predicate_.pps(site2, site3, site1, 1))
1479               return false;
1480             circle_formation_functor_.pps(site2, site3, site1, 1, circle);
1481           } else {
1482             // (segment, point, segment) sites.
1483             if (!circle_existence_predicate_.pss(site2, site1, site3, 2))
1484               return false;
1485             circle_formation_functor_.pss(site2, site1, site3, 2, circle);
1486           }
1487         } else {
1488           if (!site3.is_segment()) {
1489             // (segment, segment, point) sites.
1490             if (!circle_existence_predicate_.pss(site3, site1, site2, 3))
1491               return false;
1492             circle_formation_functor_.pss(site3, site1, site2, 3, circle);
1493           } else {
1494             // (segment, segment, segment) sites.
1495             if (!circle_existence_predicate_.sss(site1, site2, site3))
1496               return false;
1497             circle_formation_functor_.sss(site1, site2, site3, circle);
1498           }
1499         }
1500       }
1501       if (lies_outside_vertical_segment(circle, site1) ||
1502           lies_outside_vertical_segment(circle, site2) ||
1503           lies_outside_vertical_segment(circle, site3)) {
1504         return false;
1505       }
1506       return true;
1507     }
1508
1509    private:
1510     bool lies_outside_vertical_segment(
1511         const circle_type& c, const site_type& s) {
1512       if (!s.is_segment() || !is_vertical(s)) {
1513         return false;
1514       }
1515       fpt_type y0 = to_fpt(s.is_inverse() ? s.y1() : s.y0());
1516       fpt_type y1 = to_fpt(s.is_inverse() ? s.y0() : s.y1());
1517       return ulp_cmp(c.y(), y0, ULPS) == ulp_cmp_type::LESS ||
1518              ulp_cmp(c.y(), y1, ULPS) == ulp_cmp_type::MORE;
1519     }
1520
1521    private:
1522     to_fpt_converter to_fpt;
1523     ulp_cmp_type ulp_cmp;
1524     circle_existence_predicate_type circle_existence_predicate_;
1525     circle_formation_functor_type circle_formation_functor_;
1526   };
1527 };
1528 }  // detail
1529 }  // polygon
1530 }  // boost
1531
1532 #endif  // BOOST_POLYGON_DETAIL_VORONOI_PREDICATES