Boost-Commit :

Date view	Thread view	Subject view	Author view

Subject: [Boost-commit] svn:boost r77544 - sandbox/gtl/boost/polygon/detail
From: sydorchuk.andriy_at_[hidden]
Date: 2012-03-25 15:08:16

Next message: vicente.botet_at_[hidden]: "[Boost-commit] svn:boost r77545 - trunk/boost/chrono/detail/inlined/win"
Previous message: vicente.botet_at_[hidden]: "[Boost-commit] svn:boost r77543 - in trunk/libs/utility/doc: . html"

Author: asydorchuk
Date: 2012-03-25 15:08:13 EDT (Sun, 25 Mar 2012)
New Revision: 77544
URL: http://svn.boost.org/trac/boost/changeset/77544

Log:
Code styling.
Text files modified:
   sandbox/gtl/boost/polygon/detail/voronoi_predicates.hpp | 2587 ++++++++++++++++++++-------------------
   sandbox/gtl/boost/polygon/detail/voronoi_robust_fpt.hpp | 2
   2 files changed, 1298 insertions(+), 1291 deletions(-)

Modified: sandbox/gtl/boost/polygon/detail/voronoi_predicates.hpp
==============================================================================
--- sandbox/gtl/boost/polygon/detail/voronoi_predicates.hpp (original)
+++ sandbox/gtl/boost/polygon/detail/voronoi_predicates.hpp 2012-03-25 15:08:13 EDT (Sun, 25 Mar 2012)
@@ -21,1356 +21,1363 @@
template <typename CTYPE_TRAITS>
class voronoi_predicates {
public:
- typedef typename CTYPE_TRAITS::int_type int_type;
- typedef typename CTYPE_TRAITS::int_x2_type int_x2_type;
- typedef typename CTYPE_TRAITS::uint_x2_type uint_x2_type;
- typedef typename CTYPE_TRAITS::big_int_type big_int_type;
- typedef typename CTYPE_TRAITS::fpt_type fpt_type;
- typedef typename CTYPE_TRAITS::efpt_type efpt_type;
- typedef typename CTYPE_TRAITS::ulp_cmp_type ulp_cmp_type;
- typedef typename CTYPE_TRAITS::to_fpt_converter_type to_fpt_converter;
- typedef typename CTYPE_TRAITS::to_efpt_converter_type to_efpt_converter;
-
- enum {
- ULPS = CTYPE_TRAITS::ULPS,
- ULPSx2 = ULPS * 2
- };
-
- template <typename Point>
- static bool is_vertical(const Point &point1, const Point &point2) {
- return point1.x() == point2.x();
+ typedef typename CTYPE_TRAITS::int_type int_type;
+ typedef typename CTYPE_TRAITS::int_x2_type int_x2_type;
+ typedef typename CTYPE_TRAITS::uint_x2_type uint_x2_type;
+ typedef typename CTYPE_TRAITS::big_int_type big_int_type;
+ typedef typename CTYPE_TRAITS::fpt_type fpt_type;
+ typedef typename CTYPE_TRAITS::efpt_type efpt_type;
+ typedef typename CTYPE_TRAITS::ulp_cmp_type ulp_cmp_type;
+ typedef typename CTYPE_TRAITS::to_fpt_converter_type to_fpt_converter;
+ typedef typename CTYPE_TRAITS::to_efpt_converter_type to_efpt_converter;
+
+ enum {
+ ULPS = CTYPE_TRAITS::ULPS,
+ ULPSx2 = ULPS * 2
+ };
+
+ template <typename Point>
+ static bool is_vertical(const Point &point1, const Point &point2) {
+ return point1.x() == point2.x();
+ }
+
+ template <typename Site>
+ static bool is_vertical(const Site &site) {
+ return is_vertical(site.point0(), site.point1());
+ }
+
+ // Compute robust cross_product: a1 * b2 - b1 * a2.
+ // It was mathematically proven that the result is correct
+ // with epsilon relative error equal to 1EPS.
+ template <typename T>
+ static fpt_type robust_cross_product(T a1_, T b1_, T a2_, T b2_) {
+ static to_fpt_converter to_fpt;
+ uint_x2_type a1 = static_cast<uint_x2_type>(is_neg(a1_) ? -a1_ : a1_);
+ uint_x2_type b1 = static_cast<uint_x2_type>(is_neg(b1_) ? -b1_ : b1_);
+ uint_x2_type a2 = static_cast<uint_x2_type>(is_neg(a2_) ? -a2_ : a2_);
+ uint_x2_type b2 = static_cast<uint_x2_type>(is_neg(b2_) ? -b2_ : b2_);
+
+ uint_x2_type l = a1 * b2;
+ uint_x2_type r = b1 * a2;
+
+ if (is_neg(a1_) ^ is_neg(b2_)) {
+ if (is_neg(a2_) ^ is_neg(b1_))
+ return (l > r) ? -to_fpt(l - r) : to_fpt(r - l);
+ else
+ return -to_fpt(l + r);
+ } else {
+ if (is_neg(a2_) ^ is_neg(b1_))
+ return to_fpt(l + r);
+ else
+ return (l < r) ? -to_fpt(r - l) : to_fpt(l - r);
     }
+ }

- template <typename Site>
- static bool is_vertical(const Site &site) {
- return is_vertical(site.point0(), site.point1());
- }
+ typedef struct orientation_test {
+ public:
+ // Represents orientation test result.
+ enum Orientation {
+ RIGHT = -1,
+ COLLINEAR = 0,
+ LEFT = 1
+ };

- // Compute robust cross_product: a1 * b2 - b1 * a2.
- // It was mathematically proven that the result is correct
- // with epsilon relative error equal to 1EPS.
+ // Value is a determinant of two vectors (e.g. x1 * y2 - x2 * y1).
+ // Return orientation based on the sign of the determinant.
     template <typename T>
- static fpt_type robust_cross_product(T a1_, T b1_, T a2_, T b2_) {
- static to_fpt_converter to_fpt;
- uint_x2_type a1 = static_cast<uint_x2_type>(is_neg(a1_) ? -a1_ : a1_);
- uint_x2_type b1 = static_cast<uint_x2_type>(is_neg(b1_) ? -b1_ : b1_);
- uint_x2_type a2 = static_cast<uint_x2_type>(is_neg(a2_) ? -a2_ : a2_);
- uint_x2_type b2 = static_cast<uint_x2_type>(is_neg(b2_) ? -b2_ : b2_);
-
- uint_x2_type l = a1 * b2;
- uint_x2_type r = b1 * a2;
-
- if (is_neg(a1_) ^ is_neg(b2_)) {
- if (is_neg(a2_) ^ is_neg(b1_))
- return (l > r) ? -to_fpt(l - r) : to_fpt(r - l);
- else
- return -to_fpt(l + r);
- } else {
- if (is_neg(a2_) ^ is_neg(b1_))
- return to_fpt(l + r);
- else
- return (l < r) ? -to_fpt(r - l) : to_fpt(l - r);
- }
+ static Orientation eval(T value) {
+ if (is_zero(value)) return COLLINEAR;
+ return (is_neg(value)) ? RIGHT : LEFT;
     }

- typedef struct orientation_test {
- public:
- // Represents orientation test result.
- enum Orientation {
- RIGHT = -1,
- COLLINEAR = 0,
- LEFT = 1
- };
-
- // Value is a determinant of two vectors (e.g. x1 * y2 - x2 * y1).
- // Return orientation based on the sign of the determinant.
- template <typename T>
- static Orientation eval(T value) {
- if (is_zero(value)) return COLLINEAR;
- return (is_neg(value)) ? RIGHT : LEFT;
- }
-
- template <typename T>
- static Orientation eval(T dif_x1_, T dif_y1_, T dif_x2_, T dif_y2_) {
- return eval(robust_cross_product(dif_x1_, dif_y1_,
- dif_x2_, dif_y2_));
- }
-
- template <typename Point>
- static Orientation eval(const Point &point1,
- const Point &point2,
- const Point &point3) {
- int_x2_type dx1 = static_cast<int_x2_type>(point1.x()) -
- static_cast<int_x2_type>(point2.x());
- int_x2_type dx2 = static_cast<int_x2_type>(point2.x()) -
- static_cast<int_x2_type>(point3.x());
- int_x2_type dy1 = static_cast<int_x2_type>(point1.y()) -
- static_cast<int_x2_type>(point2.y());
- int_x2_type dy2 = static_cast<int_x2_type>(point2.y()) -
- static_cast<int_x2_type>(point3.y());
- return eval(robust_cross_product(dx1, dy1, dx2, dy2));
- }
- } ot;
+ template <typename T>
+ static Orientation eval(T dif_x1_, T dif_y1_, T dif_x2_, T dif_y2_) {
+ return eval(robust_cross_product(dif_x1_, dif_y1_, dif_x2_, dif_y2_));
+ }

     template <typename Point>
- class point_comparison_predicate {
- public:
- typedef Point point_type;
-
- bool operator()(const point_type &lhs, const point_type &rhs) const {
- if (lhs.x() == rhs.x()) {
- return lhs.y() < rhs.y();
- }
- return lhs.x() < rhs.x();
- }
- };
+ static Orientation eval(const Point &point1,
+ const Point &point2,
+ const Point &point3) {
+ int_x2_type dx1 = static_cast<int_x2_type>(point1.x()) -
+ static_cast<int_x2_type>(point2.x());
+ int_x2_type dx2 = static_cast<int_x2_type>(point2.x()) -
+ static_cast<int_x2_type>(point3.x());
+ int_x2_type dy1 = static_cast<int_x2_type>(point1.y()) -
+ static_cast<int_x2_type>(point2.y());
+ int_x2_type dy2 = static_cast<int_x2_type>(point2.y()) -
+ static_cast<int_x2_type>(point3.y());
+ return eval(robust_cross_product(dx1, dy1, dx2, dy2));
+ }
+ } ot;

- template <typename Site, typename Circle>
- class event_comparison_predicate {
- public:
- typedef Site site_type;
- typedef Circle circle_type;
-
- bool operator()(const site_type &lhs, const site_type &rhs) const {
- if (lhs.x0() != rhs.x0()) {
- return lhs.x0() < rhs.x0();
- }
- if (!lhs.is_segment()) {
- if (!rhs.is_segment()) {
- return lhs.y0() < rhs.y0();
- }
- if (is_vertical(rhs)) {
- return lhs.y0() <= rhs.y0();
- }
- return true;
- } else {
- if (is_vertical(rhs)) {
- if(is_vertical(lhs)) {
- return lhs.y0() < rhs.y0();
- }
- return false;
- }
- if (is_vertical(lhs)) {
- return true;
- }
- if (lhs.y0() != rhs.y0()) {
- return lhs.y0() < rhs.y0();
- }
- return ot::eval(lhs.point1(), lhs.point0(), rhs.point1()) == ot::LEFT;
- }
- }
+ template <typename Point>
+ class point_comparison_predicate {
+ public:
+ typedef Point point_type;
+
+ bool operator()(const point_type &lhs, const point_type &rhs) const {
+ if (lhs.x() == rhs.x())
+ return lhs.y() < rhs.y();
+ return lhs.x() < rhs.x();
+ }
+ };

- bool operator()(const site_type &lhs, const circle_type &rhs) const {
- typename ulp_cmp_type::Result xCmp =
- ulp_cmp(to_fpt(lhs.x()), to_fpt(rhs.lower_x()), ULPS);
- if (xCmp != ulp_cmp_type::EQUAL) {
- return xCmp == ulp_cmp_type::LESS;
- }
- typename ulp_cmp_type::Result yCmp =
- ulp_cmp(to_fpt(lhs.y()), to_fpt(rhs.lower_y()), ULPS);
- return yCmp == ulp_cmp_type::LESS;
- }
+ template <typename Site, typename Circle>
+ class event_comparison_predicate {
+ public:
+ typedef Site site_type;
+ typedef Circle circle_type;
+
+ bool operator()(const site_type &lhs, const site_type &rhs) const {
+ if (lhs.x0() != rhs.x0())
+ return lhs.x0() < rhs.x0();
+ if (!lhs.is_segment()) {
+ if (!rhs.is_segment())
+ return lhs.y0() < rhs.y0();
+ if (is_vertical(rhs))
+ return lhs.y0() <= rhs.y0();
+ return true;
+ } else {
+ if (is_vertical(rhs)) {
+ if(is_vertical(lhs))
+ return lhs.y0() < rhs.y0();
+ return false;
+ }
+ if (is_vertical(lhs))
+ return true;
+ if (lhs.y0() != rhs.y0())
+ return lhs.y0() < rhs.y0();
+ return ot::eval(lhs.point1(), lhs.point0(), rhs.point1()) == ot::LEFT;
+ }
+ }

- bool operator()(const circle_type &lhs, const site_type &rhs) const {
- typename ulp_cmp_type::Result xCmp =
- ulp_cmp(to_fpt(lhs.lower_x()), to_fpt(rhs.x()), ULPS);
- if (xCmp != ulp_cmp_type::EQUAL) {
- return xCmp == ulp_cmp_type::LESS;
- }
- typename ulp_cmp_type::Result yCmp =
- ulp_cmp(to_fpt(lhs.lower_y()), to_fpt(rhs.y()), ULPS);
- return yCmp == ulp_cmp_type::LESS;
- }
+ bool operator()(const site_type &lhs, const circle_type &rhs) const {
+ typename ulp_cmp_type::Result xCmp =
+ ulp_cmp(to_fpt(lhs.x()), to_fpt(rhs.lower_x()), ULPS);
+ if (xCmp != ulp_cmp_type::EQUAL)
+ return xCmp == ulp_cmp_type::LESS;
+ typename ulp_cmp_type::Result yCmp =
+ ulp_cmp(to_fpt(lhs.y()), to_fpt(rhs.lower_y()), ULPS);
+ return yCmp == ulp_cmp_type::LESS;
+ }

- bool operator()(const circle_type &lhs, const circle_type &rhs) const {
- typename ulp_cmp_type::Result xCmp =
- ulp_cmp(to_fpt(lhs.lower_x()), to_fpt(rhs.lower_x()), ULPSx2);
- if (xCmp != ulp_cmp_type::EQUAL) {
- return xCmp == ulp_cmp_type::LESS;
- }
- typename ulp_cmp_type::Result yCmp =
- ulp_cmp(to_fpt(lhs.lower_y()), to_fpt(rhs.lower_y()), ULPSx2);
- return yCmp == ulp_cmp_type::LESS;
- }
+ bool operator()(const circle_type &lhs, const site_type &rhs) const {
+ typename ulp_cmp_type::Result xCmp =
+ ulp_cmp(to_fpt(lhs.lower_x()), to_fpt(rhs.x()), ULPS);
+ if (xCmp != ulp_cmp_type::EQUAL)
+ return xCmp == ulp_cmp_type::LESS;
+ typename ulp_cmp_type::Result yCmp =
+ ulp_cmp(to_fpt(lhs.lower_y()), to_fpt(rhs.y()), ULPS);
+ return yCmp == ulp_cmp_type::LESS;
+ }

- private:
- ulp_cmp_type ulp_cmp;
- to_fpt_converter to_fpt;
- };
+ bool operator()(const circle_type &lhs, const circle_type &rhs) const {
+ typename ulp_cmp_type::Result xCmp =
+ ulp_cmp(to_fpt(lhs.lower_x()), to_fpt(rhs.lower_x()), ULPSx2);
+ if (xCmp != ulp_cmp_type::EQUAL)
+ return xCmp == ulp_cmp_type::LESS;
+ typename ulp_cmp_type::Result yCmp =
+ ulp_cmp(to_fpt(lhs.lower_y()), to_fpt(rhs.lower_y()), ULPSx2);
+ return yCmp == ulp_cmp_type::LESS;
+ }

- template <typename Site>
- class distance_predicate {
- public:
- typedef Site site_type;
-
- // Returns true if a horizontal line going through a new site intersects
- // right arc at first, else returns false. If horizontal line goes
- // through intersection point of the given two arcs returns false also.
- bool operator()(const site_type &left_site,
- const site_type &right_site,
- const site_type &new_site) const {
- if (!left_site.is_segment()) {
- if (!right_site.is_segment()) {
- return pp(left_site, right_site, new_site);
- } else {
- return ps(left_site, right_site, new_site, false);
- }
- } else {
- if (!right_site.is_segment()) {
- return ps(right_site, left_site, new_site, true);
- } else {
- return ss(left_site, right_site, new_site);
- }
- }
+ private:
+ ulp_cmp_type ulp_cmp;
+ to_fpt_converter to_fpt;
+ };
+
+ template <typename Site>
+ class distance_predicate {
+ public:
+ typedef Site site_type;
+
+ // Returns true if a horizontal line going through a new site intersects
+ // right arc at first, else returns false. If horizontal line goes
+ // through intersection point of the given two arcs returns false also.
+ bool operator()(const site_type &left_site,
+ const site_type &right_site,
+ const site_type &new_site) const {
+ if (!left_site.is_segment()) {
+ if (!right_site.is_segment()) {
+ return pp(left_site, right_site, new_site);
+ } else {
+ return ps(left_site, right_site, new_site, false);
         }
-
- private:
- // Represents the result of the epsilon robust predicate. If the
- // result is undefined some further processing is usually required.
- enum kPredicateResult {
- LESS = -1,
- UNDEFINED = 0,
- MORE = 1
- };
-
- typedef typename Site::point_type point_type;
-
- // Robust predicate, avoids using high-precision libraries.
- // Returns true if a horizontal line going through the new point site
- // intersects right arc at first, else returns false. If horizontal line
- // goes through intersection point of the given two arcs returns false.
- bool pp(const site_type &left_site,
- const site_type &right_site,
- const site_type &new_site) const {
- const point_type &left_point = left_site.point0();
- const point_type &right_point = right_site.point0();
- const point_type &new_point = new_site.point0();
- if (left_point.x() > right_point.x()) {
- if (new_point.y() <= left_point.y())
- return false;
- } else if (left_point.x() < right_point.x()) {
- if (new_point.y() >= right_point.y())
- return true;
- } else {
- return static_cast<int_x2_type>(left_point.y()) +
- static_cast<int_x2_type>(right_point.y()) <
- static_cast<int_x2_type>(new_point.y()) * 2;
- }
-
- fpt_type dist1 = find_distance_to_point_arc(left_site, new_point);
- fpt_type dist2 = find_distance_to_point_arc(right_site, new_point);
-
- // The undefined ulp range is equal to 3EPS + 3EPS <= 6ULP.
- return dist1 < dist2;
+ } else {
+ if (!right_site.is_segment()) {
+ return ps(right_site, left_site, new_site, true);
+ } else {
+ return ss(left_site, right_site, new_site);
         }
+ }
+ }

- bool ps(const site_type &left_site, const site_type &right_site,
- const site_type &new_site, bool reverse_order) const {
- kPredicateResult fast_res = fast_ps(
- left_site, right_site, new_site, reverse_order);
- if (fast_res != UNDEFINED) {
- return (fast_res == LESS);
- }
-
- fpt_type dist1 = find_distance_to_point_arc(left_site, new_site.point0());
- fpt_type dist2 = find_distance_to_segment_arc(right_site, new_site.point0());
-
- // The undefined ulp range is equal to 3EPS + 7EPS <= 10ULP.
- return reverse_order ^ (dist1 < dist2);
- }
+ private:
+ // Represents the result of the epsilon robust predicate. If the
+ // result is undefined some further processing is usually required.
+ enum kPredicateResult {
+ LESS = -1,
+ UNDEFINED = 0,
+ MORE = 1
+ };

- bool ss(const site_type &left_site,
- const site_type &right_site,
- const site_type &new_site) const {
- // Handle temporary segment sites.
- if (left_site.point0() == right_site.point0() &&
- left_site.point1() == right_site.point1()) {
- return ot::eval(left_site.point0(),
- left_site.point1(),
- new_site.point0()) == ot::LEFT;
- }
+ typedef typename Site::point_type point_type;

- fpt_type dist1 = find_distance_to_segment_arc(left_site, new_site.point0());
- fpt_type dist2 = find_distance_to_segment_arc(right_site, new_site.point0());
+ // Robust predicate, avoids using high-precision libraries.
+ // Returns true if a horizontal line going through the new point site
+ // intersects right arc at first, else returns false. If horizontal line
+ // goes through intersection point of the given two arcs returns false.
+ bool pp(const site_type &left_site,
+ const site_type &right_site,
+ const site_type &new_site) const {
+ const point_type &left_point = left_site.point0();
+ const point_type &right_point = right_site.point0();
+ const point_type &new_point = new_site.point0();
+ if (left_point.x() > right_point.x()) {
+ if (new_point.y() <= left_point.y())
+ return false;
+ } else if (left_point.x() < right_point.x()) {
+ if (new_point.y() >= right_point.y())
+ return true;
+ } else {
+ return static_cast<int_x2_type>(left_point.y()) +
+ static_cast<int_x2_type>(right_point.y()) <
+ static_cast<int_x2_type>(new_point.y()) * 2;
+ }

- // The undefined ulp range is equal to 7EPS + 7EPS <= 14ULP.
- return dist1 < dist2;
- }
+ fpt_type dist1 = find_distance_to_point_arc(left_site, new_point);
+ fpt_type dist2 = find_distance_to_point_arc(right_site, new_point);

- fpt_type find_distance_to_point_arc(const site_type &site,
- const point_type &point) const {
- fpt_type dx = to_fpt(site.x()) - to_fpt(point.x());
- fpt_type dy = to_fpt(site.y()) - to_fpt(point.y());
- // The relative error is at most 3EPS.
- return (dx * dx + dy * dy) / (to_fpt(2.0) * dx);
- }
+ // The undefined ulp range is equal to 3EPS + 3EPS <= 6ULP.
+ return dist1 < dist2;
+ }

- fpt_type find_distance_to_segment_arc(const site_type &site,
- const point_type &point) const {
- if (is_vertical(site)) {
- return (to_fpt(site.x()) - to_fpt(point.x())) * to_fpt(0.5);
- } else {
- const point_type &segment0 = site.point0(true);
- const point_type &segment1 = site.point1(true);
- fpt_type a1 = to_fpt(segment1.x()) - to_fpt(segment0.x());
- fpt_type b1 = to_fpt(segment1.y()) - to_fpt(segment0.y());
- fpt_type a3 = to_fpt(point.x()) - to_fpt(segment0.x());
- fpt_type b3 = to_fpt(point.y()) - to_fpt(segment0.y());
- fpt_type k = get_sqrt(a1 * a1 + b1 * b1);
- // Avoid subtraction while computing k.
- if (!is_neg(b1)) {
- k = to_fpt(1.0) / (b1 + k);
- } else {
- k = (k - b1) / (a1 * a1);
- }
- // The relative error is at most 7EPS.
- return robust_cross_product(a1, b1, a3, b3) * k;
- }
- }
+ bool ps(const site_type &left_site, const site_type &right_site,
+ const site_type &new_site, bool reverse_order) const {
+ kPredicateResult fast_res = fast_ps(
+ left_site, right_site, new_site, reverse_order);
+ if (fast_res != UNDEFINED)
+ return (fast_res == LESS);
+
+ fpt_type dist1 = find_distance_to_point_arc(
+ left_site, new_site.point0());
+ fpt_type dist2 = find_distance_to_segment_arc(
+ right_site, new_site.point0());

- kPredicateResult fast_ps(const site_type &left_site, const site_type &right_site,
- const site_type &new_site, bool reverse_order) const {
- const point_type &site_point = left_site.point0();
- const point_type &segment_start = right_site.point0(true);
- const point_type &segment_end = right_site.point1(true);
- const point_type &new_point = new_site.point0();
- if (ot::eval(segment_start, segment_end, new_point) != ot::RIGHT) {
- return (!right_site.is_inverse()) ? LESS : MORE;
- }
+ // The undefined ulp range is equal to 3EPS + 7EPS <= 10ULP.
+ return reverse_order ^ (dist1 < dist2);
+ }

- fpt_type dif_x = to_fpt(new_point.x()) - to_fpt(site_point.x());
- fpt_type dif_y = to_fpt(new_point.y()) - to_fpt(site_point.y());
- fpt_type a = to_fpt(segment_end.x()) - to_fpt(segment_start.x());
- fpt_type b = to_fpt(segment_end.y()) - to_fpt(segment_start.y());
-
- if (is_vertical(right_site)) {
- if (new_point.y() < site_point.y() && !reverse_order)
- return MORE;
- else if (new_point.y() > site_point.y() && reverse_order)
- return LESS;
- return UNDEFINED;
- } else {
- typename ot::Orientation orientation = ot::eval(a, b, dif_x, dif_y);
- if (orientation == ot::LEFT) {
- if (!right_site.is_inverse())
- return reverse_order ? LESS : UNDEFINED;
- return reverse_order ? UNDEFINED : MORE;
- }
- }
+ bool ss(const site_type &left_site,
+ const site_type &right_site,
+ const site_type &new_site) const {
+ // Handle temporary segment sites.
+ if (left_site.point0() == right_site.point0() &&
+ left_site.point1() == right_site.point1()) {
+ return ot::eval(left_site.point0(),
+ left_site.point1(),
+ new_site.point0()) == ot::LEFT;
+ }
+
+ fpt_type dist1 = find_distance_to_segment_arc(
+ left_site, new_site.point0());
+ fpt_type dist2 = find_distance_to_segment_arc(
+ right_site, new_site.point0());

- fpt_type fast_left_expr = a * (dif_y + dif_x) * (dif_y - dif_x);
- fpt_type fast_right_expr = (to_fpt(2.0) * b) * dif_x * dif_y;
- typename ulp_cmp_type::Result expr_cmp = ulp_cmp(fast_left_expr, fast_right_expr, 4);
- if (expr_cmp != ulp_cmp_type::EQUAL) {
- if ((expr_cmp == ulp_cmp_type::MORE) ^ reverse_order)
- return reverse_order ? LESS : MORE;
- return UNDEFINED;
- }
- return UNDEFINED;
- }
+ // The undefined ulp range is equal to 7EPS + 7EPS <= 14ULP.
+ return dist1 < dist2;
+ }

- private:
- ulp_cmp_type ulp_cmp;
- to_fpt_converter to_fpt;
- };
+ fpt_type find_distance_to_point_arc(
+ const site_type &site, const point_type &point) const {
+ fpt_type dx = to_fpt(site.x()) - to_fpt(point.x());
+ fpt_type dy = to_fpt(site.y()) - to_fpt(point.y());
+ // The relative error is at most 3EPS.
+ return (dx * dx + dy * dy) / (to_fpt(2.0) * dx);
+ }

- template <typename Node>
- class node_comparison_predicate {
- public:
- typedef Node node_type;
- typedef typename Node::site_type site_type;
- typedef typename site_type::coordinate_type coordinate_type;
- typedef distance_predicate<site_type> distance_predicate_type;
-
- // Compares nodes in the balanced binary search tree. Nodes are
- // compared based on the y coordinates of the arcs intersection points.
- // Nodes with less y coordinate of the intersection point go first.
- // Comparison is only called during the new site events processing.
- // That's why one of the nodes will always lie on the sweepline and may
- // be represented as a straight horizontal line.
- bool operator() (const node_type &node1,
- const node_type &node2) const {
- // Get x coordinate of the rightmost site from both nodes.
- const site_type &site1 = get_comparison_site(node1);
- const site_type &site2 = get_comparison_site(node2);
-
- if (site1.x() < site2.x()) {
- // The second node contains a new site.
- return predicate_(node1.left_site(), node1.right_site(), site2);
- } else if (site1.x() > site2.x()) {
- // The first node contains a new site.
- return !predicate_(node2.left_site(), node2.right_site(), site1);
- } else {
- // This checks were evaluated experimentally.
- if (site1.index() == site2.index()) {
- // Both nodes are new (inserted during the same site event processing).
- return get_comparison_y(node1) < get_comparison_y(node2);
- } else if (site1.index() < site2.index()) {
- std::pair<coordinate_type, int> y1 = get_comparison_y(node1, false);
- std::pair<coordinate_type, int> y2 = get_comparison_y(node2, true);
- if (y1.first != y2.first) return y1.first < y2.first;
- return (!site1.is_segment()) ? (y1.second < 0) : false;
- } else {
- std::pair<coordinate_type, int> y1 = get_comparison_y(node1, true);
- std::pair<coordinate_type, int> y2 = get_comparison_y(node2, false);
- if (y1.first != y2.first) return y1.first < y2.first;
- return (!site2.is_segment()) ? (y2.second > 0) : true;
- }
- }
+ fpt_type find_distance_to_segment_arc(
+ const site_type &site, const point_type &point) const {
+ if (is_vertical(site)) {
+ return (to_fpt(site.x()) - to_fpt(point.x())) * to_fpt(0.5);
+ } else {
+ const point_type &segment0 = site.point0(true);
+ const point_type &segment1 = site.point1(true);
+ fpt_type a1 = to_fpt(segment1.x()) - to_fpt(segment0.x());
+ fpt_type b1 = to_fpt(segment1.y()) - to_fpt(segment0.y());
+ fpt_type a3 = to_fpt(point.x()) - to_fpt(segment0.x());
+ fpt_type b3 = to_fpt(point.y()) - to_fpt(segment0.y());
+ fpt_type k = get_sqrt(a1 * a1 + b1 * b1);
+ // Avoid subtraction while computing k.
+ if (!is_neg(b1)) {
+ k = to_fpt(1.0) / (b1 + k);
+ } else {
+ k = (k - b1) / (a1 * a1);
         }
+ // The relative error is at most 7EPS.
+ return robust_cross_product(a1, b1, a3, b3) * k;
+ }
+ }

- private:
- // Get the newer site.
- const site_type &get_comparison_site(const node_type &node) const {
- if (node.left_site().index() > node.right_site().index()) {
- return node.left_site();
- }
- return node.right_site();
- }
+ kPredicateResult fast_ps(
+ const site_type &left_site, const site_type &right_site,
+ const site_type &new_site, bool reverse_order) const {
+ const point_type &site_point = left_site.point0();
+ const point_type &segment_start = right_site.point0(true);
+ const point_type &segment_end = right_site.point1(true);
+ const point_type &new_point = new_site.point0();
+
+ if (ot::eval(segment_start, segment_end, new_point) != ot::RIGHT)
+ return (!right_site.is_inverse()) ? LESS : MORE;
+
+ fpt_type dif_x = to_fpt(new_point.x()) - to_fpt(site_point.x());
+ fpt_type dif_y = to_fpt(new_point.y()) - to_fpt(site_point.y());
+ fpt_type a = to_fpt(segment_end.x()) - to_fpt(segment_start.x());
+ fpt_type b = to_fpt(segment_end.y()) - to_fpt(segment_start.y());
+
+ if (is_vertical(right_site)) {
+ if (new_point.y() < site_point.y() && !reverse_order)
+ return MORE;
+ else if (new_point.y() > site_point.y() && reverse_order)
+ return LESS;
+ return UNDEFINED;
+ } else {
+ typename ot::Orientation orientation = ot::eval(a, b, dif_x, dif_y);
+ if (orientation == ot::LEFT) {
+ if (!right_site.is_inverse())
+ return reverse_order ? LESS : UNDEFINED;
+ return reverse_order ? UNDEFINED : MORE;
+ }
+ }
+
+ fpt_type fast_left_expr = a * (dif_y + dif_x) * (dif_y - dif_x);
+ fpt_type fast_right_expr = (to_fpt(2.0) * b) * dif_x * dif_y;
+ typename ulp_cmp_type::Result expr_cmp =
+ ulp_cmp(fast_left_expr, fast_right_expr, 4);
+ if (expr_cmp != ulp_cmp_type::EQUAL) {
+ if ((expr_cmp == ulp_cmp_type::MORE) ^ reverse_order)
+ return reverse_order ? LESS : MORE;
+ return UNDEFINED;
+ }
+ return UNDEFINED;
+ }

- // Get comparison pair: y coordinate and direction of the newer site.
- std::pair<coordinate_type, int> get_comparison_y(
- const node_type &node, bool is_new_node = true) const {
- if (node.left_site().index() == node.right_site().index()) {
- return std::make_pair(node.left_site().y(), 0);
- }
- if (node.left_site().index() > node.right_site().index()) {
- if (!is_new_node &&
- node.left_site().is_segment() &&
- is_vertical(node.left_site())) {
- return std::make_pair(node.left_site().y1(), 1);
- }
- return std::make_pair(node.left_site().y(), 1);
- }
- return std::make_pair(node.right_site().y(), -1);
+ private:
+ ulp_cmp_type ulp_cmp;
+ to_fpt_converter to_fpt;
+ };
+
+ template <typename Node>
+ class node_comparison_predicate {
+ public:
+ typedef Node node_type;
+ typedef typename Node::site_type site_type;
+ typedef typename site_type::coordinate_type coordinate_type;
+ typedef distance_predicate<site_type> distance_predicate_type;
+
+ // Compares nodes in the balanced binary search tree. Nodes are
+ // compared based on the y coordinates of the arcs intersection points.
+ // Nodes with less y coordinate of the intersection point go first.
+ // Comparison is only called during the new site events processing.
+ // That's why one of the nodes will always lie on the sweepline and may
+ // be represented as a straight horizontal line.
+ bool operator() (const node_type &node1,
+ const node_type &node2) const {
+ // Get x coordinate of the rightmost site from both nodes.
+ const site_type &site1 = get_comparison_site(node1);
+ const site_type &site2 = get_comparison_site(node2);
+
+ if (site1.x() < site2.x()) {
+ // The second node contains a new site.
+ return predicate_(node1.left_site(), node1.right_site(), site2);
+ } else if (site1.x() > site2.x()) {
+ // The first node contains a new site.
+ return !predicate_(node2.left_site(), node2.right_site(), site1);
+ } else {
+ // This checks were evaluated experimentally.
+ if (site1.index() == site2.index()) {
+ // Both nodes are new (inserted during same site event processing).
+ return get_comparison_y(node1) < get_comparison_y(node2);
+ } else if (site1.index() < site2.index()) {
+ std::pair<coordinate_type, int> y1 = get_comparison_y(node1, false);
+ std::pair<coordinate_type, int> y2 = get_comparison_y(node2, true);
+ if (y1.first != y2.first) return y1.first < y2.first;
+ return (!site1.is_segment()) ? (y1.second < 0) : false;
+ } else {
+ std::pair<coordinate_type, int> y1 = get_comparison_y(node1, true);
+ std::pair<coordinate_type, int> y2 = get_comparison_y(node2, false);
+ if (y1.first != y2.first) return y1.first < y2.first;
+ return (!site2.is_segment()) ? (y2.second > 0) : true;
         }
+ }
+ }

- distance_predicate_type predicate_;
- };
+ private:
+ // Get the newer site.
+ const site_type &get_comparison_site(const node_type &node) const {
+ if (node.left_site().index() > node.right_site().index()) {
+ return node.left_site();
+ }
+ return node.right_site();
+ }

- template <typename Site>
- class circle_existence_predicate {
- public:
- typedef typename Site::point_type point_type;
- typedef Site site_type;
-
- bool ppp(const site_type &site1,
- const site_type &site2,
- const site_type &site3) const {
- return ot::eval(site1.point0(), site2.point0(), site3.point0()) == ot::RIGHT;
- }
-
- bool pps(const site_type &site1,
- const site_type &site2,
- const site_type &site3,
- int segment_index) const {
- if (segment_index != 2) {
- typename ot::Orientation orient1 = ot::eval(site1.point0(),
- site2.point0(), site3.point0(true));
- typename ot::Orientation orient2 = ot::eval(site1.point0(),
- site2.point0(), site3.point1(true));
- if (segment_index == 1 && site1.x0() >= site2.x0()) {
- if (orient1 != ot::RIGHT)
- return false;
- } else if (segment_index == 3 && site2.x0() >= site1.x0()) {
- if (orient2 != ot::RIGHT)
- return false;
- } else if (orient1 != ot::RIGHT && orient2 != ot::RIGHT) {
- return false;
- }
- } else {
- if (site3.point0(true) == site1.point0() &&
- site3.point1(true) == site2.point0())
- return false;
- }
- return true;
- }
+ // Get comparison pair: y coordinate and direction of the newer site.
+ std::pair<coordinate_type, int> get_comparison_y(
+ const node_type &node, bool is_new_node = true) const {
+ if (node.left_site().index() == node.right_site().index()) {
+ return std::make_pair(node.left_site().y(), 0);
+ }
+ if (node.left_site().index() > node.right_site().index()) {
+ if (!is_new_node &&
+ node.left_site().is_segment() &&
+ is_vertical(node.left_site())) {
+ return std::make_pair(node.left_site().y1(), 1);
+ }
+ return std::make_pair(node.left_site().y(), 1);
+ }
+ return std::make_pair(node.right_site().y(), -1);
+ }

- bool pss(const site_type &site1,
- const site_type &site2,
- const site_type &site3,
- int point_index) const {
- if (site2.point0() == site3.point0() &&
- site2.point1() == site3.point1()) {
- return false;
- }
- if (point_index == 2) {
- if (!site2.is_inverse() && site3.is_inverse())
- return false;
- if (site2.is_inverse() == site3.is_inverse() &&
- ot::eval(site2.point0(true),
- site1.point0(),
- site3.point1(true)) != ot::RIGHT)
- return false;
- }
- return true;
- }
+ distance_predicate_type predicate_;
+ };

- bool sss(const site_type &site1,
- const site_type &site2,
- const site_type &site3) const {
- if (site1.point0() == site2.point0() &&
- site1.point1() == site2.point1())
- return false;
- if (site2.point0() == site3.point0() &&
- site2.point1() == site3.point1())
- return false;
- return true;
- }
- };
+ template <typename Site>
+ class circle_existence_predicate {
+ public:
+ typedef typename Site::point_type point_type;
+ typedef Site site_type;
+
+ bool ppp(const site_type &site1,
+ const site_type &site2,
+ const site_type &site3) const {
+ return ot::eval(site1.point0(), site2.point0(), site3.point0()) ==
+ ot::RIGHT;
+ }

- template <typename Site, typename Circle>
- class mp_circle_formation_functor {
- public:
- typedef typename Site::point_type point_type;
- typedef Site site_type;
- typedef Circle circle_type;
- typedef robust_sqrt_expr<big_int_type, efpt_type, to_efpt_converter> robust_sqrt_expr_type;
-
- void ppp(const site_type &site1,
- const site_type &site2,
- const site_type &site3,
- circle_type &circle,
- bool recompute_c_x = true,
- bool recompute_c_y = true,
- bool recompute_lower_x = true) {
- big_int_type dif_x[3], dif_y[3], sum_x[2], sum_y[2];
- dif_x[0] = static_cast<int_x2_type>(site1.x()) -
- static_cast<int_x2_type>(site2.x());
- dif_x[1] = static_cast<int_x2_type>(site2.x()) -
- static_cast<int_x2_type>(site3.x());
- dif_x[2] = static_cast<int_x2_type>(site1.x()) -
- static_cast<int_x2_type>(site3.x());
- dif_y[0] = static_cast<int_x2_type>(site1.y()) -
- static_cast<int_x2_type>(site2.y());
- dif_y[1] = static_cast<int_x2_type>(site2.y()) -
- static_cast<int_x2_type>(site3.y());
- dif_y[2] = static_cast<int_x2_type>(site1.y()) -
- static_cast<int_x2_type>(site3.y());
- sum_x[0] = static_cast<int_x2_type>(site1.x()) +
- static_cast<int_x2_type>(site2.x());
- sum_x[1] = static_cast<int_x2_type>(site2.x()) +
- static_cast<int_x2_type>(site3.x());
- sum_y[0] = static_cast<int_x2_type>(site1.y()) +
- static_cast<int_x2_type>(site2.y());
- sum_y[1] = static_cast<int_x2_type>(site2.y()) +
- static_cast<int_x2_type>(site3.y());
- fpt_type inv_denom = to_fpt(0.5) / to_fpt(static_cast<big_int_type>(
- dif_x[0] * dif_y[1] - dif_x[1] * dif_y[0]));
- big_int_type numer1 = dif_x[0] * sum_x[0] + dif_y[0] * sum_y[0];
- big_int_type numer2 = dif_x[1] * sum_x[1] + dif_y[1] * sum_y[1];
-
- if (recompute_c_x || recompute_lower_x) {
- big_int_type c_x = numer1 * dif_y[1] - numer2 * dif_y[0];
- circle.x(to_fpt(c_x) * inv_denom);
-
- if (recompute_lower_x) {
- // Evaluate radius of the circle.
- big_int_type sqr_r = (dif_x[0] * dif_x[0] + dif_y[0] * dif_y[0]) *
- (dif_x[1] * dif_x[1] + dif_y[1] * dif_y[1]) *
- (dif_x[2] * dif_x[2] + dif_y[2] * dif_y[2]);
- fpt_type r = get_sqrt(to_fpt(sqr_r));
-
- // If c_x >= 0 then lower_x = c_x + r,
- // else lower_x = (c_x * c_x - r * r) / (c_x - r).
- // To guarantee epsilon relative error.
- if (!is_neg(circle.x())) {
- if (!is_neg(inv_denom)) {
- circle.lower_x(circle.x() + r * inv_denom);
- } else {
- circle.lower_x(circle.x() - r * inv_denom);
- }
- } else {
- big_int_type numer = c_x * c_x - sqr_r;
- fpt_type lower_x = to_fpt(numer) * inv_denom /
- (to_fpt(c_x) + r);
- circle.lower_x(lower_x);
- }
- }
- }
+ bool pps(const site_type &site1,
+ const site_type &site2,
+ const site_type &site3,
+ int segment_index) const {
+ if (segment_index != 2) {
+ typename ot::Orientation orient1 = ot::eval(site1.point0(),
+ site2.point0(), site3.point0(true));
+ typename ot::Orientation orient2 = ot::eval(site1.point0(),
+ site2.point0(), site3.point1(true));
+ if (segment_index == 1 && site1.x0() >= site2.x0()) {
+ if (orient1 != ot::RIGHT)
+ return false;
+ } else if (segment_index == 3 && site2.x0() >= site1.x0()) {
+ if (orient2 != ot::RIGHT)
+ return false;
+ } else if (orient1 != ot::RIGHT && orient2 != ot::RIGHT) {
+ return false;
+ }
+ } else {
+ if (site3.point0(true) == site1.point0() &&
+ site3.point1(true) == site2.point0())
+ return false;
+ }
+ return true;
+ }

- if (recompute_c_y) {
- big_int_type c_y = numer2 * dif_x[0] - numer1 * dif_x[1];
- circle.y(to_fpt(c_y) * inv_denom);
- }
- }
+ bool pss(const site_type &site1,
+ const site_type &site2,
+ const site_type &site3,
+ int point_index) const {
+ if (site2.point0() == site3.point0() &&
+ site2.point1() == site3.point1()) {
+ return false;
+ }
+ if (point_index == 2) {
+ if (!site2.is_inverse() && site3.is_inverse())
+ return false;
+ if (site2.is_inverse() == site3.is_inverse() &&
+ ot::eval(site2.point0(true),
+ site1.point0(),
+ site3.point1(true)) != ot::RIGHT)
+ return false;
+ }
+ return true;
+ }

- // Recompute parameters of the circle event using high-precision library.
- void pps(const site_type &site1,
- const site_type &site2,
- const site_type &site3,
- int segment_index,
- circle_type &c_event,
- bool recompute_c_x = true,
- bool recompute_c_y = true,
- bool recompute_lower_x = true) {
- big_int_type cA[4], cB[4];
- big_int_type line_a = static_cast<int_x2_type>(site3.point1(true).y()) -
- static_cast<int_x2_type>(site3.point0(true).y());
- big_int_type line_b = static_cast<int_x2_type>(site3.point0(true).x()) -
- static_cast<int_x2_type>(site3.point1(true).x());
- big_int_type segm_len = line_a * line_a + line_b * line_b;
- big_int_type vec_x = static_cast<int_x2_type>(site2.y()) -
- static_cast<int_x2_type>(site1.y());
- big_int_type vec_y = static_cast<int_x2_type>(site1.x()) -
- static_cast<int_x2_type>(site2.x());
- big_int_type sum_x = static_cast<int_x2_type>(site1.x()) +
- static_cast<int_x2_type>(site2.x());
- big_int_type sum_y = static_cast<int_x2_type>(site1.y()) +
- static_cast<int_x2_type>(site2.y());
- big_int_type teta = line_a * vec_x + line_b * vec_y;
- big_int_type denom = vec_x * line_b - vec_y * line_a;
-
- big_int_type dif0 = static_cast<int_x2_type>(site3.point1().y()) -
- static_cast<int_x2_type>(site1.y());
- big_int_type dif1 = static_cast<int_x2_type>(site1.x()) -
- static_cast<int_x2_type>(site3.point1().x());
- big_int_type A = line_a * dif1 - line_b * dif0;
- dif0 = static_cast<int_x2_type>(site3.point1().y()) -
- static_cast<int_x2_type>(site2.y());
- dif1 = static_cast<int_x2_type>(site2.x()) -
- static_cast<int_x2_type>(site3.point1().x());
- big_int_type B = line_a * dif1 - line_b * dif0;
- big_int_type sum_AB = A + B;
-
- if (is_zero(denom)) {
- big_int_type numer = teta * teta - sum_AB * sum_AB;
- big_int_type denom = teta * sum_AB;
- cA[0] = denom * sum_x * 2 + numer * vec_x;
- cB[0] = segm_len;
- cA[1] = denom * sum_AB * 2 + numer * teta;
- cB[1] = 1;
- cA[2] = denom * sum_y * 2 + numer * vec_y;
- fpt_type inv_denom = to_fpt(1.0) / to_fpt(denom);
- if (recompute_c_x) {
- c_event.x(to_fpt(0.25) * to_fpt(cA[0]) * inv_denom);
- }
- if (recompute_c_y) {
- c_event.y(to_fpt(0.25) * to_fpt(cA[2]) * inv_denom);
- }
- if (recompute_lower_x) {
- c_event.lower_x(to_fpt(0.25) * to_fpt(sqrt_expr_.eval2(cA, cB)) * inv_denom /
- get_sqrt(to_fpt(segm_len)));
- }
- return;
- }
+ bool sss(const site_type &site1,
+ const site_type &site2,
+ const site_type &site3) const {
+ if (site1.point0() == site2.point0() && site1.point1() == site2.point1())
+ return false;
+ if (site2.point0() == site3.point0() && site2.point1() == site3.point1())
+ return false;
+ return true;
+ }
+ };

- big_int_type det = (teta * teta + denom * denom) * A * B * 4;
- fpt_type inv_denom_sqr = to_fpt(1.0) / to_fpt(denom);
- inv_denom_sqr *= inv_denom_sqr;
-
- if (recompute_c_x || recompute_lower_x) {
- cA[0] = sum_x * denom * denom + teta * sum_AB * vec_x;
- cB[0] = 1;
- cA[1] = (segment_index == 2) ? -vec_x : vec_x;
- cB[1] = det;
- if (recompute_c_x) {
- c_event.x(to_fpt(0.5) * to_fpt(sqrt_expr_.eval2(cA, cB)) * inv_denom_sqr);
- }
+ template <typename Site, typename Circle>
+ class mp_circle_formation_functor {
+ public:
+ typedef typename Site::point_type point_type;
+ typedef Site site_type;
+ typedef Circle circle_type;
+ typedef robust_sqrt_expr<big_int_type, efpt_type, to_efpt_converter>
+ robust_sqrt_expr_type;
+
+ void ppp(const site_type &site1,
+ const site_type &site2,
+ const site_type &site3,
+ circle_type &circle,
+ bool recompute_c_x = true,
+ bool recompute_c_y = true,
+ bool recompute_lower_x = true) {
+ big_int_type dif_x[3], dif_y[3], sum_x[2], sum_y[2];
+ dif_x[0] = static_cast<int_x2_type>(site1.x()) -
+ static_cast<int_x2_type>(site2.x());
+ dif_x[1] = static_cast<int_x2_type>(site2.x()) -
+ static_cast<int_x2_type>(site3.x());
+ dif_x[2] = static_cast<int_x2_type>(site1.x()) -
+ static_cast<int_x2_type>(site3.x());
+ dif_y[0] = static_cast<int_x2_type>(site1.y()) -
+ static_cast<int_x2_type>(site2.y());
+ dif_y[1] = static_cast<int_x2_type>(site2.y()) -
+ static_cast<int_x2_type>(site3.y());
+ dif_y[2] = static_cast<int_x2_type>(site1.y()) -
+ static_cast<int_x2_type>(site3.y());
+ sum_x[0] = static_cast<int_x2_type>(site1.x()) +
+ static_cast<int_x2_type>(site2.x());
+ sum_x[1] = static_cast<int_x2_type>(site2.x()) +
+ static_cast<int_x2_type>(site3.x());
+ sum_y[0] = static_cast<int_x2_type>(site1.y()) +
+ static_cast<int_x2_type>(site2.y());
+ sum_y[1] = static_cast<int_x2_type>(site2.y()) +
+ static_cast<int_x2_type>(site3.y());
+ fpt_type inv_denom = to_fpt(0.5) / to_fpt(static_cast<big_int_type>(
+ dif_x[0] * dif_y[1] - dif_x[1] * dif_y[0]));
+ big_int_type numer1 = dif_x[0] * sum_x[0] + dif_y[0] * sum_y[0];
+ big_int_type numer2 = dif_x[1] * sum_x[1] + dif_y[1] * sum_y[1];
+
+ if (recompute_c_x || recompute_lower_x) {
+ big_int_type c_x = numer1 * dif_y[1] - numer2 * dif_y[0];
+ circle.x(to_fpt(c_x) * inv_denom);
+
+ if (recompute_lower_x) {
+ // Evaluate radius of the circle.
+ big_int_type sqr_r = (dif_x[0] * dif_x[0] + dif_y[0] * dif_y[0]) *
+ (dif_x[1] * dif_x[1] + dif_y[1] * dif_y[1]) *
+ (dif_x[2] * dif_x[2] + dif_y[2] * dif_y[2]);
+ fpt_type r = get_sqrt(to_fpt(sqr_r));
+
+ // If c_x >= 0 then lower_x = c_x + r,
+ // else lower_x = (c_x * c_x - r * r) / (c_x - r).
+ // To guarantee epsilon relative error.
+ if (!is_neg(circle.x())) {
+ if (!is_neg(inv_denom)) {
+ circle.lower_x(circle.x() + r * inv_denom);
+ } else {
+ circle.lower_x(circle.x() - r * inv_denom);
             }
+ } else {
+ big_int_type numer = c_x * c_x - sqr_r;
+ fpt_type lower_x = to_fpt(numer) * inv_denom / (to_fpt(c_x) + r);
+ circle.lower_x(lower_x);
+ }
+ }
+ }
+
+ if (recompute_c_y) {
+ big_int_type c_y = numer2 * dif_x[0] - numer1 * dif_x[1];
+ circle.y(to_fpt(c_y) * inv_denom);
+ }
+ }

- if (recompute_c_y || recompute_lower_x) {
- cA[2] = sum_y * denom * denom + teta * sum_AB * vec_y;
- cB[2] = 1;
- cA[3] = (segment_index == 2) ? -vec_y : vec_y;
- cB[3] = det;
- if (recompute_c_y) {
- c_event.y(to_fpt(0.5) * to_fpt(sqrt_expr_.eval2(&cA[2], &cB[2])) *
- inv_denom_sqr);
- }
- }
+ // Recompute parameters of the circle event using high-precision library.
+ void pps(const site_type &site1,
+ const site_type &site2,
+ const site_type &site3,
+ int segment_index,
+ circle_type &c_event,
+ bool recompute_c_x = true,
+ bool recompute_c_y = true,
+ bool recompute_lower_x = true) {
+ big_int_type cA[4], cB[4];
+ big_int_type line_a = static_cast<int_x2_type>(site3.point1(true).y()) -
+ static_cast<int_x2_type>(site3.point0(true).y());
+ big_int_type line_b = static_cast<int_x2_type>(site3.point0(true).x()) -
+ static_cast<int_x2_type>(site3.point1(true).x());
+ big_int_type segm_len = line_a * line_a + line_b * line_b;
+ big_int_type vec_x = static_cast<int_x2_type>(site2.y()) -
+ static_cast<int_x2_type>(site1.y());
+ big_int_type vec_y = static_cast<int_x2_type>(site1.x()) -
+ static_cast<int_x2_type>(site2.x());
+ big_int_type sum_x = static_cast<int_x2_type>(site1.x()) +
+ static_cast<int_x2_type>(site2.x());
+ big_int_type sum_y = static_cast<int_x2_type>(site1.y()) +
+ static_cast<int_x2_type>(site2.y());
+ big_int_type teta = line_a * vec_x + line_b * vec_y;
+ big_int_type denom = vec_x * line_b - vec_y * line_a;
+
+ big_int_type dif0 = static_cast<int_x2_type>(site3.point1().y()) -
+ static_cast<int_x2_type>(site1.y());
+ big_int_type dif1 = static_cast<int_x2_type>(site1.x()) -
+ static_cast<int_x2_type>(site3.point1().x());
+ big_int_type A = line_a * dif1 - line_b * dif0;
+ dif0 = static_cast<int_x2_type>(site3.point1().y()) -
+ static_cast<int_x2_type>(site2.y());
+ dif1 = static_cast<int_x2_type>(site2.x()) -
+ static_cast<int_x2_type>(site3.point1().x());
+ big_int_type B = line_a * dif1 - line_b * dif0;
+ big_int_type sum_AB = A + B;
+
+ if (is_zero(denom)) {
+ big_int_type numer = teta * teta - sum_AB * sum_AB;
+ big_int_type denom = teta * sum_AB;
+ cA[0] = denom * sum_x * 2 + numer * vec_x;
+ cB[0] = segm_len;
+ cA[1] = denom * sum_AB * 2 + numer * teta;
+ cB[1] = 1;
+ cA[2] = denom * sum_y * 2 + numer * vec_y;
+ fpt_type inv_denom = to_fpt(1.0) / to_fpt(denom);
+ if (recompute_c_x)
+ c_event.x(to_fpt(0.25) * to_fpt(cA[0]) * inv_denom);
+ if (recompute_c_y)
+ c_event.y(to_fpt(0.25) * to_fpt(cA[2]) * inv_denom);
+ if (recompute_lower_x) {
+ c_event.lower_x(to_fpt(0.25) * to_fpt(sqrt_expr_.eval2(cA, cB)) *
+ inv_denom / get_sqrt(to_fpt(segm_len)));
+ }
+ return;
+ }
+
+ big_int_type det = (teta * teta + denom * denom) * A * B * 4;
+ fpt_type inv_denom_sqr = to_fpt(1.0) / to_fpt(denom);
+ inv_denom_sqr *= inv_denom_sqr;
+
+ if (recompute_c_x || recompute_lower_x) {
+ cA[0] = sum_x * denom * denom + teta * sum_AB * vec_x;
+ cB[0] = 1;
+ cA[1] = (segment_index == 2) ? -vec_x : vec_x;
+ cB[1] = det;
+ if (recompute_c_x) {
+ c_event.x(to_fpt(0.5) * to_fpt(sqrt_expr_.eval2(cA, cB)) *
+ inv_denom_sqr);
+ }
+ }
+
+ if (recompute_c_y || recompute_lower_x) {
+ cA[2] = sum_y * denom * denom + teta * sum_AB * vec_y;
+ cB[2] = 1;
+ cA[3] = (segment_index == 2) ? -vec_y : vec_y;
+ cB[3] = det;
+ if (recompute_c_y) {
+ c_event.y(to_fpt(0.5) * to_fpt(sqrt_expr_.eval2(&cA[2], &cB[2])) *
+ inv_denom_sqr);
+ }
+ }
+
+ if (recompute_lower_x) {
+ cB[0] = cB[0] * segm_len;
+ cB[1] = cB[1] * segm_len;
+ cA[2] = sum_AB * (denom * denom + teta * teta);
+ cB[2] = 1;
+ cA[3] = (segment_index == 2) ? -teta : teta;
+ cB[3] = det;
+ c_event.lower_x(to_fpt(0.5) * to_fpt(sqrt_expr_.eval4(cA, cB)) *
+ inv_denom_sqr / get_sqrt(to_fpt(segm_len)));
+ }
+ }

- if (recompute_lower_x) {
- cB[0] = cB[0] * segm_len;
- cB[1] = cB[1] * segm_len;
- cA[2] = sum_AB * (denom * denom + teta * teta);
- cB[2] = 1;
- cA[3] = (segment_index == 2) ? -teta : teta;
- cB[3] = det;
- c_event.lower_x(to_fpt(0.5) * to_fpt(sqrt_expr_.eval4(cA, cB)) * inv_denom_sqr /
- get_sqrt(to_fpt(segm_len)));
- }
+ // Recompute parameters of the circle event using high-precision library.
+ void pss(const site_type &site1,
+ const site_type &site2,
+ const site_type &site3,
+ int point_index,
+ circle_type &c_event,
+ bool recompute_c_x = true,
+ bool recompute_c_y = true,
+ bool recompute_lower_x = true) {
+ big_int_type a[2], b[2], c[2], cA[4], cB[4];
+ const point_type &segm_start1 = site2.point1(true);
+ const point_type &segm_end1 = site2.point0(true);
+ const point_type &segm_start2 = site3.point0(true);
+ const point_type &segm_end2 = site3.point1(true);
+ a[0] = static_cast<int_x2_type>(segm_end1.x()) -
+ static_cast<int_x2_type>(segm_start1.x());
+ b[0] = static_cast<int_x2_type>(segm_end1.y()) -
+ static_cast<int_x2_type>(segm_start1.y());
+ a[1] = static_cast<int_x2_type>(segm_end2.x()) -
+ static_cast<int_x2_type>(segm_start2.x());
+ b[1] = static_cast<int_x2_type>(segm_end2.y()) -
+ static_cast<int_x2_type>(segm_start2.y());
+ big_int_type orientation = a[1] * b[0] - a[0] * b[1];
+ if (is_zero(orientation)) {
+ fpt_type denom = to_fpt(2.0) * to_fpt(
+ static_cast<big_int_type>(a[0] * a[0] + b[0] * b[0]));
+ c[0] = b[0] * (static_cast<int_x2_type>(segm_start2.x()) -
+ static_cast<int_x2_type>(segm_start1.x())) -
+ a[0] * (static_cast<int_x2_type>(segm_start2.y()) -
+ static_cast<int_x2_type>(segm_start1.y()));
+ big_int_type dx = a[0] * (static_cast<int_x2_type>(site1.y()) -
+ static_cast<int_x2_type>(segm_start1.y())) -
+ b[0] * (static_cast<int_x2_type>(site1.x()) -
+ static_cast<int_x2_type>(segm_start1.x()));
+ big_int_type dy = b[0] * (static_cast<int_x2_type>(site1.x()) -
+ static_cast<int_x2_type>(segm_start2.x())) -
+ a[0] * (static_cast<int_x2_type>(site1.y()) -
+ static_cast<int_x2_type>(segm_start2.y()));
+ cB[0] = dx * dy;
+ cB[1] = 1;
+
+ if (recompute_c_y) {
+ cA[0] = b[0] * ((point_index == 2) ? 2 : -2);
+ cA[1] = a[0] * a[0] * (static_cast<int_x2_type>(segm_start1.y()) +
+ static_cast<int_x2_type>(segm_start2.y())) -
+ a[0] * b[0] * (static_cast<int_x2_type>(segm_start1.x()) +
+ static_cast<int_x2_type>(segm_start2.x()) -
+ static_cast<int_x2_type>(site1.x()) * 2) +
+ b[0] * b[0] * (static_cast<int_x2_type>(site1.y()) * 2);
+ fpt_type c_y = to_fpt(sqrt_expr_.eval2(cA, cB));
+ c_event.y(c_y / denom);
+ }
+
+ if (recompute_c_x || recompute_lower_x) {
+ cA[0] = a[0] * ((point_index == 2) ? 2 : -2);
+ cA[1] = b[0] * b[0] * (static_cast<int_x2_type>(segm_start1.x()) +
+ static_cast<int_x2_type>(segm_start2.x())) -
+ a[0] * b[0] * (static_cast<int_x2_type>(segm_start1.y()) +
+ static_cast<int_x2_type>(segm_start2.y()) -
+ static_cast<int_x2_type>(site1.y()) * 2) +
+ a[0] * a[0] * (static_cast<int_x2_type>(site1.x()) * 2);
+
+ if (recompute_c_x) {
+ fpt_type c_x = to_fpt(sqrt_expr_.eval2(cA, cB));
+ c_event.x(c_x / denom);
+ }
+
+ if (recompute_lower_x) {
+ cA[2] = is_neg(c[0]) ? -c[0] : c[0];
+ cB[2] = a[0] * a[0] + b[0] * b[0];
+ fpt_type lower_x = to_fpt(sqrt_expr_.eval3(cA, cB));
+ c_event.lower_x(lower_x / denom);
+ }
+ }
+ return;
+ }
+ c[0] = b[0] * segm_end1.x() - a[0] * segm_end1.y();
+ c[1] = a[1] * segm_end2.y() - b[1] * segm_end2.x();
+ big_int_type ix = a[0] * c[1] + a[1] * c[0];
+ big_int_type iy = b[0] * c[1] + b[1] * c[0];
+ big_int_type dx = ix - orientation * site1.x();
+ big_int_type dy = iy - orientation * site1.y();
+ if (is_zero(dx) && is_zero(dy)) {
+ fpt_type denom = to_fpt(orientation);
+ fpt_type c_x = to_fpt(ix) / denom;
+ fpt_type c_y = to_fpt(iy) / denom;
+ c_event = circle_type(c_x, c_y, c_x);
+ return;
+ }
+
+ big_int_type sign = ((point_index == 2) ? 1 : -1) *
+ (is_neg(orientation) ? 1 : -1);
+ cA[0] = a[1] * -dx + b[1] * -dy;
+ cA[1] = a[0] * -dx + b[0] * -dy;
+ cA[2] = sign;
+ cA[3] = 0;
+ cB[0] = a[0] * a[0] + b[0] * b[0];
+ cB[1] = a[1] * a[1] + b[1] * b[1];
+ cB[2] = a[0] * a[1] + b[0] * b[1];
+ cB[3] = (a[0] * dy - b[0] * dx) * (a[1] * dy - b[1] * dx) * -2;
+ fpt_type temp = to_fpt(
+ sqrt_expr_evaluator_pss4<big_int_type, efpt_type>(cA, cB));
+ fpt_type denom = temp * to_fpt(orientation);
+
+ if (recompute_c_y) {
+ cA[0] = b[1] * (dx * dx + dy * dy) - iy * (dx * a[1] + dy * b[1]);
+ cA[1] = b[0] * (dx * dx + dy * dy) - iy * (dx * a[0] + dy * b[0]);
+ cA[2] = iy * sign;
+ fpt_type cy = to_fpt(
+ sqrt_expr_evaluator_pss4<big_int_type, efpt_type>(cA, cB));
+ c_event.y(cy / denom);
+ }
+
+ if (recompute_c_x || recompute_lower_x) {
+ cA[0] = a[1] * (dx * dx + dy * dy) - ix * (dx * a[1] + dy * b[1]);
+ cA[1] = a[0] * (dx * dx + dy * dy) - ix * (dx * a[0] + dy * b[0]);
+ cA[2] = ix * sign;
+
+ if (recompute_c_x) {
+ fpt_type cx = to_fpt(
+ sqrt_expr_evaluator_pss4<big_int_type, efpt_type>(cA, cB));
+ c_event.x(cx / denom);
+ }
+
+ if (recompute_lower_x) {
+ cA[3] = orientation * (dx * dx + dy * dy) * (is_neg(temp) ? -1 : 1);
+ fpt_type lower_x = to_fpt(
+ sqrt_expr_evaluator_pss4<big_int_type, efpt_type>(cA, cB));
+ c_event.lower_x(lower_x / denom);
         }
+ }
+ }

- // Recompute parameters of the circle event using high-precision library.
- void pss(const site_type &site1,
- const site_type &site2,
- const site_type &site3,
- int point_index,
- circle_type &c_event,
- bool recompute_c_x = true,
- bool recompute_c_y = true,
- bool recompute_lower_x = true) {
- big_int_type a[2], b[2], c[2], cA[4], cB[4];
- const point_type &segm_start1 = site2.point1(true);
- const point_type &segm_end1 = site2.point0(true);
- const point_type &segm_start2 = site3.point0(true);
- const point_type &segm_end2 = site3.point1(true);
- a[0] = static_cast<int_x2_type>(segm_end1.x()) -
- static_cast<int_x2_type>(segm_start1.x());
- b[0] = static_cast<int_x2_type>(segm_end1.y()) -
- static_cast<int_x2_type>(segm_start1.y());
- a[1] = static_cast<int_x2_type>(segm_end2.x()) -
- static_cast<int_x2_type>(segm_start2.x());
- b[1] = static_cast<int_x2_type>(segm_end2.y()) -
- static_cast<int_x2_type>(segm_start2.y());
- big_int_type orientation = a[1] * b[0] - a[0] * b[1];
- if (is_zero(orientation)) {
- fpt_type denom = to_fpt(2.0) * to_fpt(
- static_cast<big_int_type>(a[0] * a[0] + b[0] * b[0]));
- c[0] = b[0] * (static_cast<int_x2_type>(segm_start2.x()) -
- static_cast<int_x2_type>(segm_start1.x())) -
- a[0] * (static_cast<int_x2_type>(segm_start2.y()) -
- static_cast<int_x2_type>(segm_start1.y()));
- big_int_type dx = a[0] * (static_cast<int_x2_type>(site1.y()) -
- static_cast<int_x2_type>(segm_start1.y())) -
- b[0] * (static_cast<int_x2_type>(site1.x()) -
- static_cast<int_x2_type>(segm_start1.x()));
- big_int_type dy = b[0] * (static_cast<int_x2_type>(site1.x()) -
- static_cast<int_x2_type>(segm_start2.x())) -
- a[0] * (static_cast<int_x2_type>(site1.y()) -
- static_cast<int_x2_type>(segm_start2.y()));
- cB[0] = dx * dy;
- cB[1] = 1;
-
- if (recompute_c_y) {
- cA[0] = b[0] * ((point_index == 2) ? 2 : -2);
- cA[1] = a[0] * a[0] * (static_cast<int_x2_type>(segm_start1.y()) +
- static_cast<int_x2_type>(segm_start2.y())) -
- a[0] * b[0] * (static_cast<int_x2_type>(segm_start1.x()) +
- static_cast<int_x2_type>(segm_start2.x()) -
- static_cast<int_x2_type>(site1.x()) * 2) +
- b[0] * b[0] * (static_cast<int_x2_type>(site1.y()) * 2);
- fpt_type c_y = to_fpt(sqrt_expr_.eval2(cA, cB));
- c_event.y(c_y / denom);
- }
-
- if (recompute_c_x || recompute_lower_x) {
- cA[0] = a[0] * ((point_index == 2) ? 2 : -2);
- cA[1] = b[0] * b[0] * (static_cast<int_x2_type>(segm_start1.x()) +
- static_cast<int_x2_type>(segm_start2.x())) -
- a[0] * b[0] * (static_cast<int_x2_type>(segm_start1.y()) +
- static_cast<int_x2_type>(segm_start2.y()) -
- static_cast<int_x2_type>(site1.y()) * 2) +
- a[0] * a[0] * (static_cast<int_x2_type>(site1.x()) * 2);
-
- if (recompute_c_x) {
- fpt_type c_x = to_fpt(sqrt_expr_.eval2(cA, cB));
- c_event.x(c_x / denom);
- }
-
- if (recompute_lower_x) {
- cA[2] = is_neg(c[0]) ? -c[0] : c[0];
- cB[2] = a[0] * a[0] + b[0] * b[0];
- fpt_type lower_x = to_fpt(sqrt_expr_.eval3(cA, cB));
- c_event.lower_x(lower_x / denom);
- }
- }
- return;
- }
- c[0] = b[0] * segm_end1.x() -
- a[0] * segm_end1.y();
- c[1] = a[1] * segm_end2.y() -
- b[1] * segm_end2.x();
- big_int_type ix = a[0] * c[1] + a[1] * c[0];
- big_int_type iy = b[0] * c[1] + b[1] * c[0];
- big_int_type dx = ix - orientation * site1.x();
- big_int_type dy = iy - orientation * site1.y();
- if (is_zero(dx) && is_zero(dy)) {
- fpt_type denom = to_fpt(orientation);
- fpt_type c_x = to_fpt(ix) / denom;
- fpt_type c_y = to_fpt(iy) / denom;
- c_event = circle_type(c_x, c_y, c_x);
- return;
- }
-
- big_int_type sign = ((point_index == 2) ? 1 : -1) * (is_neg(orientation) ? 1 : -1);
- cA[0] = a[1] * -dx + b[1] * -dy;
- cA[1] = a[0] * -dx + b[0] * -dy;
- cA[2] = sign;
- cA[3] = 0;
- cB[0] = a[0] * a[0] + b[0] * b[0];
- cB[1] = a[1] * a[1] + b[1] * b[1];
- cB[2] = a[0] * a[1] + b[0] * b[1];
- cB[3] = (a[0] * dy - b[0] * dx) * (a[1] * dy - b[1] * dx) * -2;
- fpt_type temp = to_fpt(sqrt_expr_evaluator_pss4<big_int_type, efpt_type>(cA, cB));
- fpt_type denom = temp * to_fpt(orientation);
-
- if (recompute_c_y) {
- cA[0] = b[1] * (dx * dx + dy * dy) - iy * (dx * a[1] + dy * b[1]);
- cA[1] = b[0] * (dx * dx + dy * dy) - iy * (dx * a[0] + dy * b[0]);
- cA[2] = iy * sign;
- fpt_type cy = to_fpt(sqrt_expr_evaluator_pss4<big_int_type, efpt_type>(cA, cB));
- c_event.y(cy / denom);
- }
-
- if (recompute_c_x || recompute_lower_x) {
- cA[0] = a[1] * (dx * dx + dy * dy) - ix * (dx * a[1] + dy * b[1]);
- cA[1] = a[0] * (dx * dx + dy * dy) - ix * (dx * a[0] + dy * b[0]);
- cA[2] = ix * sign;
-
- if (recompute_c_x) {
- fpt_type cx = to_fpt(sqrt_expr_evaluator_pss4<big_int_type, efpt_type>(cA, cB));
- c_event.x(cx / denom);
- }
-
- if (recompute_lower_x) {
- cA[3] = orientation * (dx * dx + dy * dy) * (is_neg(temp) ? -1 : 1);
- fpt_type lower_x = to_fpt(sqrt_expr_evaluator_pss4<big_int_type, efpt_type>(cA, cB));
- c_event.lower_x(lower_x / denom);
- }
- }
+ // Recompute parameters of the circle event using high-precision library.
+ void sss(const site_type &site1,
+ const site_type &site2,
+ const site_type &site3,
+ circle_type &c_event,
+ bool recompute_c_x = true,
+ bool recompute_c_y = true,
+ bool recompute_lower_x = true) {
+ big_int_type a[3], b[3], c[3], cA[4], cB[4];
+ // cA - corresponds to the cross product.
+ // cB - corresponds to the squared length.
+ a[0] = static_cast<int_x2_type>(site1.x1(true)) -
+ static_cast<int_x2_type>(site1.x0(true));
+ a[1] = static_cast<int_x2_type>(site2.x1(true)) -
+ static_cast<int_x2_type>(site2.x0(true));
+ a[2] = static_cast<int_x2_type>(site3.x1(true)) -
+ static_cast<int_x2_type>(site3.x0(true));
+
+ b[0] = static_cast<int_x2_type>(site1.y1(true)) -
+ static_cast<int_x2_type>(site1.y0(true));
+ b[1] = static_cast<int_x2_type>(site2.y1(true)) -
+ static_cast<int_x2_type>(site2.y0(true));
+ b[2] = static_cast<int_x2_type>(site3.y1(true)) -
+ static_cast<int_x2_type>(site3.y0(true));
+
+ c[0] = static_cast<int_x2_type>(site1.x0(true)) *
+ static_cast<int_x2_type>(site1.y1(true)) -
+ static_cast<int_x2_type>(site1.y0(true)) *
+ static_cast<int_x2_type>(site1.x1(true));
+ c[1] = static_cast<int_x2_type>(site2.x0(true)) *
+ static_cast<int_x2_type>(site2.y1(true)) -
+ static_cast<int_x2_type>(site2.y0(true)) *
+ static_cast<int_x2_type>(site2.x1(true));
+ c[2] = static_cast<int_x2_type>(site3.x0(true)) *
+ static_cast<int_x2_type>(site3.y1(true)) -
+ static_cast<int_x2_type>(site3.y0(true)) *
+ static_cast<int_x2_type>(site3.x1(true));
+
+ for (int i = 0; i < 3; ++i)
+ cB[i] = a[i] * a[i] + b[i] * b[i];
+
+ for (int i = 0; i < 3; ++i) {
+ int j = (i+1) % 3;
+ int k = (i+2) % 3;
+ cA[i] = a[j] * b[k] - a[k] * b[j];
+ }
+ fpt_type denom = to_fpt(sqrt_expr_.eval3(cA, cB));
+
+ if (recompute_c_y) {
+ for (int i = 0; i < 3; ++i) {
+ int j = (i+1) % 3;
+ int k = (i+2) % 3;
+ cA[i] = b[j] * c[k] - b[k] * c[j];
+ }
+ fpt_type c_y = to_fpt(sqrt_expr_.eval3(cA, cB));
+ c_event.y(c_y / denom);
+ }
+
+ if (recompute_c_x || recompute_lower_x) {
+ cA[3] = 0;
+ for (int i = 0; i < 3; ++i) {
+ int j = (i+1) % 3;
+ int k = (i+2) % 3;
+ cA[i] = a[j] * c[k] - a[k] * c[j];
+ if (recompute_lower_x) {
+ cA[3] = cA[3] + cA[i] * b[i];
+ }
+ }
+
+ if (recompute_c_x) {
+ fpt_type c_x = to_fpt(sqrt_expr_.eval3(cA, cB));
+ c_event.x(c_x / denom);
+ }
+
+ if (recompute_lower_x) {
+ cB[3] = 1;
+ fpt_type lower_x = to_fpt(sqrt_expr_.eval4(cA, cB));
+ c_event.lower_x(lower_x / denom);
         }
+ }
+ }

- // Recompute parameters of the circle event using high-precision library.
- void sss(const site_type &site1,
- const site_type &site2,
- const site_type &site3,
- circle_type &c_event,
- bool recompute_c_x = true,
- bool recompute_c_y = true,
- bool recompute_lower_x = true) {
- big_int_type a[3], b[3], c[3], cA[4], cB[4];
- // cA - corresponds to the cross product.
- // cB - corresponds to the squared length.
- a[0] = static_cast<int_x2_type>(site1.x1(true)) -
- static_cast<int_x2_type>(site1.x0(true));
- a[1] = static_cast<int_x2_type>(site2.x1(true)) -
- static_cast<int_x2_type>(site2.x0(true));
- a[2] = static_cast<int_x2_type>(site3.x1(true)) -
- static_cast<int_x2_type>(site3.x0(true));
-
- b[0] = static_cast<int_x2_type>(site1.y1(true)) -
- static_cast<int_x2_type>(site1.y0(true));
- b[1] = static_cast<int_x2_type>(site2.y1(true)) -
- static_cast<int_x2_type>(site2.y0(true));
- b[2] = static_cast<int_x2_type>(site3.y1(true)) -
- static_cast<int_x2_type>(site3.y0(true));
-
- c[0] = static_cast<int_x2_type>(site1.x0(true)) *
- static_cast<int_x2_type>(site1.y1(true)) -
- static_cast<int_x2_type>(site1.y0(true)) *
- static_cast<int_x2_type>(site1.x1(true));
- c[1] = static_cast<int_x2_type>(site2.x0(true)) *
- static_cast<int_x2_type>(site2.y1(true)) -
- static_cast<int_x2_type>(site2.y0(true)) *
- static_cast<int_x2_type>(site2.x1(true));
- c[2] = static_cast<int_x2_type>(site3.x0(true)) *
- static_cast<int_x2_type>(site3.y1(true)) -
- static_cast<int_x2_type>(site3.y0(true)) *
- static_cast<int_x2_type>(site3.x1(true));
-
- for (int i = 0; i < 3; ++i) {
- cB[i] = a[i] * a[i] + b[i] * b[i];
- }
+ private:
+ // Evaluates A[3] + A[0] * sqrt(B[0]) + A[1] * sqrt(B[1]) +
+ // A[2] * sqrt(B[3] * (sqrt(B[0] * B[1]) + B[2])).
+ template <typename _int, typename _fpt>
+ _fpt sqrt_expr_evaluator_pss4(_int *A, _int *B) {
+ _int cA[4], cB[4];
+ if (is_zero(A[3])) {
+ _fpt lh = sqrt_expr_.eval2(A, B);
+ cA[0] = 1;
+ cB[0] = B[0] * B[1];
+ cA[1] = B[2];
+ cB[1] = 1;
+ _fpt rh = sqrt_expr_.eval1(A+2, B+3) *
+ get_sqrt(sqrt_expr_.eval2(cA, cB));
+ if ((!is_neg(lh) && !is_neg(rh)) || (!is_pos(lh) && !is_pos(rh)))
+ return lh + rh;
+ cA[0] = A[0] * A[0] * B[0] + A[1] * A[1] * B[1] -
+ A[2] * A[2] * B[3] * B[2];
+ cB[0] = 1;
+ cA[1] = A[0] * A[1] * 2 - A[2] * A[2] * B[3];
+ cB[1] = B[0] * B[1];
+ _fpt numer = sqrt_expr_.eval2(cA, cB);
+ return numer / (lh - rh);
+ }
+ cA[0] = 1;
+ cB[0] = B[0] * B[1];
+ cA[1] = B[2];
+ cB[1] = 1;
+ _fpt rh = sqrt_expr_.eval1(A+2, B+3) * get_sqrt(sqrt_expr_.eval2(cA, cB));
+ cA[0] = A[0];
+ cB[0] = B[0];
+ cA[1] = A[1];
+ cB[1] = B[1];
+ cA[2] = A[3];
+ cB[2] = 1;
+ _fpt lh = sqrt_expr_.eval3(cA, cB);
+ if ((!is_neg(lh) && !is_neg(rh)) || (!is_pos(lh) && !is_pos(rh)))
+ return lh + rh;
+ cA[0] = A[3] * A[0] * 2;
+ cA[1] = A[3] * A[1] * 2;
+ cA[2] = A[0] * A[0] * B[0] + A[1] * A[1] * B[1] +
+ A[3] * A[3] - A[2] * A[2] * B[2] * B[3];
+ cA[3] = A[0] * A[1] * 2 - A[2] * A[2] * B[3];
+ cB[3] = B[0] * B[1];
+ _fpt numer = sqrt_expr_evaluator_pss3<_int, _fpt>(cA, cB);
+ return numer / (lh - rh);
+ }

- for (int i = 0; i < 3; ++i) {
- int j = (i+1) % 3;
- int k = (i+2) % 3;
- cA[i] = a[j] * b[k] - a[k] * b[j];
- }
- fpt_type denom = to_fpt(sqrt_expr_.eval3(cA, cB));
+ template <typename _int, typename _fpt>
+ // Evaluates A[0] * sqrt(B[0]) + A[1] * sqrt(B[1]) +
+ // A[2] + A[3] * sqrt(B[0] * B[1]).
+ // B[3] = B[0] * B[1].
+ _fpt sqrt_expr_evaluator_pss3(_int *A, _int *B) {
+ _int cA[2], cB[2];
+ _fpt lh = sqrt_expr_.eval2(A, B);
+ _fpt rh = sqrt_expr_.eval2(A+2, B+2);
+ if ((!is_neg(lh) && !is_neg(rh)) || (!is_pos(lh) && !is_pos(rh)))
+ return lh + rh;
+ cA[0] = A[0] * A[0] * B[0] + A[1] * A[1] * B[1] -
+ A[2] * A[2] - A[3] * A[3] * B[0] * B[1];
+ cB[0] = 1;
+ cA[1] = (A[0] * A[1] - A[2] * A[3]) * 2;
+ cB[1] = B[3];
+ _fpt numer = sqrt_expr_.eval2(cA, cB);
+ return numer / (lh - rh);
+ }

- if (recompute_c_y) {
- for (int i = 0; i < 3; ++i) {
- int j = (i+1) % 3;
- int k = (i+2) % 3;
- cA[i] = b[j] * c[k] - b[k] * c[j];
- }
- fpt_type c_y = to_fpt(sqrt_expr_.eval3(cA, cB));
- c_event.y(c_y / denom);
- }
+ robust_sqrt_expr_type sqrt_expr_;
+ to_fpt_converter to_fpt;
+ };
+
+ template <typename Site, typename Circle>
+ class lazy_circle_formation_functor {
+ public:
+ typedef robust_fpt<fpt_type> robust_fpt_type;
+ typedef robust_dif<robust_fpt_type> robust_dif_type;
+ typedef typename Site::point_type point_type;
+ typedef Site site_type;
+ typedef Circle circle_type;
+ typedef mp_circle_formation_functor<site_type, circle_type>
+ exact_circle_formation_functor_type;
+
+ void ppp(const site_type &site1,
+ const site_type &site2,
+ const site_type &site3,
+ circle_type &c_event) {
+ fpt_type dif_x1 = to_fpt(site1.x()) - to_fpt(site2.x());
+ fpt_type dif_x2 = to_fpt(site2.x()) - to_fpt(site3.x());
+ fpt_type dif_y1 = to_fpt(site1.y()) - to_fpt(site2.y());
+ fpt_type dif_y2 = to_fpt(site2.y()) - to_fpt(site3.y());
+ fpt_type orientation =
+ robust_cross_product(dif_x1, dif_y1, dif_x2, dif_y2);
+ robust_fpt_type inv_orientation(to_fpt(0.5) / orientation, 2.0);
+ fpt_type sum_x1 = to_fpt(site1.x()) + to_fpt(site2.x());
+ fpt_type sum_x2 = to_fpt(site2.x()) + to_fpt(site3.x());
+ fpt_type sum_y1 = to_fpt(site1.y()) + to_fpt(site2.y());
+ fpt_type sum_y2 = to_fpt(site2.y()) + to_fpt(site3.y());
+ fpt_type dif_x3 = to_fpt(site1.x()) - to_fpt(site3.x());
+ fpt_type dif_y3 = to_fpt(site1.y()) - to_fpt(site3.y());
+ robust_dif_type c_x, c_y;
+ c_x += robust_fpt_type(dif_x1 * sum_x1 * dif_y2, 2.0);
+ c_x += robust_fpt_type(dif_y1 * sum_y1 * dif_y2, 2.0);
+ c_x -= robust_fpt_type(dif_x2 * sum_x2 * dif_y1, 2.0);
+ c_x -= robust_fpt_type(dif_y2 * sum_y2 * dif_y1, 2.0);
+ c_y += robust_fpt_type(dif_x2 * sum_x2 * dif_x1, 2.0);
+ c_y += robust_fpt_type(dif_y2 * sum_y2 * dif_x1, 2.0);
+ c_y -= robust_fpt_type(dif_x1 * sum_x1 * dif_x2, 2.0);
+ c_y -= robust_fpt_type(dif_y1 * sum_y1 * dif_x2, 2.0);
+ robust_dif_type lower_x(c_x);
+ lower_x -= robust_fpt_type(get_sqrt(
+ (dif_x1 * dif_x1 + dif_y1 * dif_y1) *
+ (dif_x2 * dif_x2 + dif_y2 * dif_y2) *
+ (dif_x3 * dif_x3 + dif_y3 * dif_y3)), 5.0);
+ c_event = circle_type(
+ c_x.dif().fpv() * inv_orientation.fpv(),
+ c_y.dif().fpv() * inv_orientation.fpv(),
+ lower_x.dif().fpv() * inv_orientation.fpv());
+ bool recompute_c_x = c_x.dif().ulp() > ULPS;
+ bool recompute_c_y = c_y.dif().ulp() > ULPS;
+ bool recompute_lower_x = lower_x.dif().ulp() > ULPS;
+ if (recompute_c_x || recompute_c_y || recompute_lower_x) {
+ exact_circle_formation_functor_.ppp(
+ site1, site2, site3, c_event,
+ recompute_c_x, recompute_c_y, recompute_lower_x);
+ }
+ }

- if (recompute_c_x || recompute_lower_x) {
- cA[3] = 0;
- for (int i = 0; i < 3; ++i) {
- int j = (i+1) % 3;
- int k = (i+2) % 3;
- cA[i] = a[j] * c[k] - a[k] * c[j];
- if (recompute_lower_x) {
- cA[3] = cA[3] + cA[i] * b[i];
- }
- }
-
- if (recompute_c_x) {
- fpt_type c_x = to_fpt(sqrt_expr_.eval3(cA, cB));
- c_event.x(c_x / denom);
- }
-
- if (recompute_lower_x) {
- cB[3] = 1;
- fpt_type lower_x = to_fpt(sqrt_expr_.eval4(cA, cB));
- c_event.lower_x(lower_x / denom);
- }
- }
+ void pps(const site_type &site1,
+ const site_type &site2,
+ const site_type &site3,
+ int segment_index,
+ circle_type &c_event) {
+ fpt_type line_a = to_fpt(site3.point1(true).y()) -
+ to_fpt(site3.point0(true).y());
+ fpt_type line_b = to_fpt(site3.point0(true).x()) -
+ to_fpt(site3.point1(true).x());
+ fpt_type vec_x = to_fpt(site2.y()) - to_fpt(site1.y());
+ fpt_type vec_y = to_fpt(site1.x()) - to_fpt(site2.x());
+ robust_fpt_type teta(
+ robust_cross_product(line_a, line_b, -vec_y, vec_x), 1.0);
+ robust_fpt_type A(robust_cross_product(
+ line_a, line_b,
+ to_fpt(site3.point1().y()) - to_fpt(site1.y()),
+ to_fpt(site1.x()) - to_fpt(site3.point1().x())), 1.0);
+ robust_fpt_type B(robust_cross_product(
+ line_a, line_b,
+ to_fpt(site3.point1().y()) - to_fpt(site2.y()),
+ to_fpt(site2.x()) - to_fpt(site3.point1().x())), 1.0);
+ robust_fpt_type denom(
+ robust_cross_product(vec_x, vec_y, line_a, line_b), 1.0);
+ robust_fpt_type inv_segm_len(to_fpt(1.0) /
+ get_sqrt(line_a * line_a + line_b * line_b), 3.0);
+ robust_dif_type t;
+ if (ot::eval(denom) == ot::COLLINEAR) {
+ t += teta / (robust_fpt_type(8.0, false) * A);
+ t -= A / (robust_fpt_type(2.0, false) * teta);
+ } else {
+ robust_fpt_type det = ((teta * teta + denom * denom) * A * B).sqrt();
+ if (segment_index == 2) {
+ t -= det / (denom * denom);
+ } else {
+ t += det / (denom * denom);
         }
+ t += teta * (A + B) / (robust_fpt_type(2.0, false) * denom * denom);
+ }
+ robust_dif_type c_x, c_y;
+ c_x += robust_fpt_type(to_fpt(0.5) * (to_fpt(site1.x()) +
+ to_fpt(site2.x())), false);
+ c_x += robust_fpt_type(vec_x, false) * t;
+ c_y += robust_fpt_type(to_fpt(0.5) * (to_fpt(site1.y()) +
+ to_fpt(site2.y())), false);
+ c_y += robust_fpt_type(vec_y, false) * t;
+ robust_dif_type r, lower_x(c_x);
+ r -= robust_fpt_type(line_a, false) * robust_fpt_type(site3.x0(), false);
+ r -= robust_fpt_type(line_b, false) * robust_fpt_type(site3.y0(), false);
+ r += robust_fpt_type(line_a, false) * c_x;
+ r += robust_fpt_type(line_b, false) * c_y;
+ if (r.pos().fpv() < r.neg().fpv())
+ r = -r;
+ lower_x += r * inv_segm_len;
+ c_event = circle_type(
+ c_x.dif().fpv(), c_y.dif().fpv(), lower_x.dif().fpv());
+ bool recompute_c_x = c_x.dif().ulp() > ULPS;
+ bool recompute_c_y = c_y.dif().ulp() > ULPS;
+ bool recompute_lower_x = lower_x.dif().ulp() > ULPS;
+ if (recompute_c_x || recompute_c_y || recompute_lower_x) {
+ exact_circle_formation_functor_.pps(
+ site1, site2, site3, segment_index, c_event,
+ recompute_c_x, recompute_c_y, recompute_lower_x);
+ }
+ }

- private:
- // Evaluates A[3] + A[0] * sqrt(B[0]) + A[1] * sqrt(B[1]) +
- // A[2] * sqrt(B[3] * (sqrt(B[0] * B[1]) + B[2])).
- template <typename _int, typename _fpt>
- _fpt sqrt_expr_evaluator_pss4(_int *A, _int *B) {
- _int cA[4], cB[4];
- if (is_zero(A[3])) {
- _fpt lh = sqrt_expr_.eval2(A, B);
- cA[0] = 1;
- cB[0] = B[0] * B[1];
- cA[1] = B[2];
- cB[1] = 1;
- _fpt rh = sqrt_expr_.eval1(A+2, B+3) * get_sqrt(sqrt_expr_.eval2(cA, cB));
- if ((!is_neg(lh) && !is_neg(rh)) || (!is_pos(lh) && !is_pos(rh))) {
- return lh + rh;
- }
- cA[0] = A[0] * A[0] * B[0] + A[1] * A[1] * B[1] -
- A[2] * A[2] * B[3] * B[2];
- cB[0] = 1;
- cA[1] = A[0] * A[1] * 2 - A[2] * A[2] * B[3];
- cB[1] = B[0] * B[1];
- _fpt numer = sqrt_expr_.eval2(cA, cB);
- return numer / (lh - rh);
- }
- cA[0] = 1;
- cB[0] = B[0] * B[1];
- cA[1] = B[2];
- cB[1] = 1;
- _fpt rh = sqrt_expr_.eval1(A+2, B+3) * get_sqrt(sqrt_expr_.eval2(cA, cB));
- cA[0] = A[0];
- cB[0] = B[0];
- cA[1] = A[1];
- cB[1] = B[1];
- cA[2] = A[3];
- cB[2] = 1;
- _fpt lh = sqrt_expr_.eval3(cA, cB);
- if ((!is_neg(lh) && !is_neg(rh)) || (!is_pos(lh) && !is_pos(rh))) {
- return lh + rh;
- }
- cA[0] = A[3] * A[0] * 2;
- cA[1] = A[3] * A[1] * 2;
- cA[2] = A[0] * A[0] * B[0] + A[1] * A[1] * B[1] +
- A[3] * A[3] - A[2] * A[2] * B[2] * B[3];
- cA[3] = A[0] * A[1] * 2 - A[2] * A[2] * B[3];
- cB[3] = B[0] * B[1];
- _fpt numer = sqrt_expr_evaluator_pss3<_int, _fpt>(cA, cB);
- return numer / (lh - rh);
- }
-
- template <typename _int, typename _fpt>
- // Evaluates A[0] * sqrt(B[0]) + A[1] * sqrt(B[1]) +
- // A[2] + A[3] * sqrt(B[0] * B[1]).
- // B[3] = B[0] * B[1].
- _fpt sqrt_expr_evaluator_pss3(_int *A, _int *B) {
- _int cA[2], cB[2];
- _fpt lh = sqrt_expr_.eval2(A, B);
- _fpt rh = sqrt_expr_.eval2(A+2, B+2);
- if ((!is_neg(lh) && !is_neg(rh)) || (!is_pos(lh) && !is_pos(rh))) {
- return lh + rh;
- }
- cA[0] = A[0] * A[0] * B[0] + A[1] * A[1] * B[1] -
- A[2] * A[2] - A[3] * A[3] * B[0] * B[1];
- cB[0] = 1;
- cA[1] = (A[0] * A[1] - A[2] * A[3]) * 2;
- cB[1] = B[3];
- _fpt numer = sqrt_expr_.eval2(cA, cB);
- return numer / (lh - rh);
+ void pss(const site_type &site1,
+ const site_type &site2,
+ const site_type &site3,
+ int point_index,
+ circle_type &c_event) {
+ const point_type &segm_start1 = site2.point1(true);
+ const point_type &segm_end1 = site2.point0(true);
+ const point_type &segm_start2 = site3.point0(true);
+ const point_type &segm_end2 = site3.point1(true);
+ fpt_type a1 = to_fpt(segm_end1.x()) - to_fpt(segm_start1.x());
+ fpt_type b1 = to_fpt(segm_end1.y()) - to_fpt(segm_start1.y());
+ fpt_type a2 = to_fpt(segm_end2.x()) - to_fpt(segm_start2.x());
+ fpt_type b2 = to_fpt(segm_end2.y()) - to_fpt(segm_start2.y());
+ bool recompute_c_x, recompute_c_y, recompute_lower_x;
+ robust_fpt_type orientation(robust_cross_product(b1, a1, b2, a2), 1.0);
+ if (ot::eval(orientation) == ot::COLLINEAR) {
+ robust_fpt_type a(a1 * a1 + b1 * b1, 2.0);
+ robust_fpt_type c(robust_cross_product(
+ b1, a1,
+ to_fpt(segm_start2.y()) - to_fpt(segm_start1.y()),
+ to_fpt(segm_start2.x()) - to_fpt(segm_start1.x())), 1.0);
+ robust_fpt_type det(
+ robust_cross_product(
+ a1, b1,
+ to_fpt(site1.x()) - to_fpt(segm_start1.x()),
+ to_fpt(site1.y()) - to_fpt(segm_start1.y())) *
+ robust_cross_product(
+ b1, a1,
+ to_fpt(site1.y()) - to_fpt(segm_start2.y()),
+ to_fpt(site1.x()) - to_fpt(segm_start2.x())),
+ 3.0);
+ robust_dif_type t;
+ t -= robust_fpt_type(a1, false) * robust_fpt_type((
+ to_fpt(segm_start1.x()) + to_fpt(segm_start2.x())) * to_fpt(0.5) -
+ to_fpt(site1.x()), false);
+ t -= robust_fpt_type(b1, false) * robust_fpt_type((
+ to_fpt(segm_start1.y()) + to_fpt(segm_start2.y())) * to_fpt(0.5) -
+ to_fpt(site1.y()), false);
+ if (point_index == 2) {
+ t += det.sqrt();
+ } else {
+ t -= det.sqrt();
         }
-
- robust_sqrt_expr_type sqrt_expr_;
- to_fpt_converter to_fpt;
- };
-
- template <typename Site, typename Circle>
- class lazy_circle_formation_functor {
- public:
- typedef robust_fpt<fpt_type> robust_fpt_type;
- typedef robust_dif<robust_fpt_type> robust_dif_type;
- typedef typename Site::point_type point_type;
- typedef Site site_type;
- typedef Circle circle_type;
- typedef mp_circle_formation_functor<site_type, circle_type>
- exact_circle_formation_functor_type;
-
- void ppp(const site_type &site1,
- const site_type &site2,
- const site_type &site3,
- circle_type &c_event) {
- fpt_type dif_x1 = to_fpt(site1.x()) - to_fpt(site2.x());
- fpt_type dif_x2 = to_fpt(site2.x()) - to_fpt(site3.x());
- fpt_type dif_y1 = to_fpt(site1.y()) - to_fpt(site2.y());
- fpt_type dif_y2 = to_fpt(site2.y()) - to_fpt(site3.y());
- fpt_type orientation = robust_cross_product(dif_x1, dif_y1, dif_x2, dif_y2);
- robust_fpt_type inv_orientation(to_fpt(0.5) / orientation, 2.0);
- fpt_type sum_x1 = to_fpt(site1.x()) + to_fpt(site2.x());
- fpt_type sum_x2 = to_fpt(site2.x()) + to_fpt(site3.x());
- fpt_type sum_y1 = to_fpt(site1.y()) + to_fpt(site2.y());
- fpt_type sum_y2 = to_fpt(site2.y()) + to_fpt(site3.y());
- fpt_type dif_x3 = to_fpt(site1.x()) - to_fpt(site3.x());
- fpt_type dif_y3 = to_fpt(site1.y()) - to_fpt(site3.y());
- robust_dif_type c_x, c_y;
- c_x += robust_fpt_type(dif_x1 * sum_x1 * dif_y2, 2.0);
- c_x += robust_fpt_type(dif_y1 * sum_y1 * dif_y2, 2.0);
- c_x -= robust_fpt_type(dif_x2 * sum_x2 * dif_y1, 2.0);
- c_x -= robust_fpt_type(dif_y2 * sum_y2 * dif_y1, 2.0);
- c_y += robust_fpt_type(dif_x2 * sum_x2 * dif_x1, 2.0);
- c_y += robust_fpt_type(dif_y2 * sum_y2 * dif_x1, 2.0);
- c_y -= robust_fpt_type(dif_x1 * sum_x1 * dif_x2, 2.0);
- c_y -= robust_fpt_type(dif_y1 * sum_y1 * dif_x2, 2.0);
- robust_dif_type lower_x(c_x);
- lower_x -= robust_fpt_type(get_sqrt((dif_x1 * dif_x1 + dif_y1 * dif_y1) *
- (dif_x2 * dif_x2 + dif_y2 * dif_y2) *
- (dif_x3 * dif_x3 + dif_y3 * dif_y3)), 5.0);
- c_event = circle_type(c_x.dif().fpv() * inv_orientation.fpv(),
- c_y.dif().fpv() * inv_orientation.fpv(),
- lower_x.dif().fpv() * inv_orientation.fpv());
- bool recompute_c_x = c_x.dif().ulp() > ULPS;
- bool recompute_c_y = c_y.dif().ulp() > ULPS;
- bool recompute_lower_x = lower_x.dif().ulp() > ULPS;
- if (recompute_c_x || recompute_c_y || recompute_lower_x) {
- exact_circle_formation_functor_.ppp(
- site1, site2, site3, c_event, recompute_c_x, recompute_c_y, recompute_lower_x);
- }
+ t /= a;
+ robust_dif_type c_x, c_y;
+ c_x += robust_fpt_type(to_fpt(0.5) * (
+ to_fpt(segm_start1.x()) + to_fpt(segm_start2.x())), false);
+ c_x += robust_fpt_type(a1, false) * t;
+ c_y += robust_fpt_type(to_fpt(0.5) * (
+ to_fpt(segm_start1.y()) + to_fpt(segm_start2.y())), false);
+ c_y += robust_fpt_type(b1, false) * t;
+ robust_dif_type lower_x(c_x);
+ if (is_neg(c)) {
+ lower_x -= robust_fpt_type(0.5, false) * c / a.sqrt();
+ } else {
+ lower_x += robust_fpt_type(0.5, false) * c / a.sqrt();
         }
-
- void pps(const site_type &site1,
- const site_type &site2,
- const site_type &site3,
- int segment_index,
- circle_type &c_event) {
- fpt_type line_a = to_fpt(site3.point1(true).y()) - to_fpt(site3.point0(true).y());
- fpt_type line_b = to_fpt(site3.point0(true).x()) - to_fpt(site3.point1(true).x());
- fpt_type vec_x = to_fpt(site2.y()) - to_fpt(site1.y());
- fpt_type vec_y = to_fpt(site1.x()) - to_fpt(site2.x());
- robust_fpt_type teta(robust_cross_product(line_a, line_b, -vec_y, vec_x), 1.0);
- robust_fpt_type A(robust_cross_product(
- line_a, line_b,
- to_fpt(site3.point1().y()) - to_fpt(site1.y()),
- to_fpt(site1.x()) - to_fpt(site3.point1().x())), 1.0);
- robust_fpt_type B(robust_cross_product(
- line_a, line_b,
- to_fpt(site3.point1().y()) - to_fpt(site2.y()),
- to_fpt(site2.x()) - to_fpt(site3.point1().x())), 1.0);
- robust_fpt_type denom(robust_cross_product(vec_x, vec_y, line_a, line_b), 1.0);
- robust_fpt_type inv_segm_len(to_fpt(1.0) / get_sqrt(line_a * line_a + line_b * line_b), 3.0);
- robust_dif_type t;
- if (ot::eval(denom) == ot::COLLINEAR) {
- t += teta / (robust_fpt_type(8.0, false) * A);
- t -= A / (robust_fpt_type(2.0, false) * teta);
- } else {
- robust_fpt_type det = ((teta * teta + denom * denom) * A * B).sqrt();
- if (segment_index == 2) {
- t -= det / (denom * denom);
- } else {
- t += det / (denom * denom);
- }
- t += teta * (A + B) / (robust_fpt_type(2.0, false) * denom * denom);
- }
- robust_dif_type c_x, c_y;
- c_x += robust_fpt_type(to_fpt(0.5) * (to_fpt(site1.x()) + to_fpt(site2.x())), false);
- c_x += robust_fpt_type(vec_x, false) * t;
- c_y += robust_fpt_type(to_fpt(0.5) * (to_fpt(site1.y()) + to_fpt(site2.y())), false);
- c_y += robust_fpt_type(vec_y, false) * t;
- robust_dif_type r, lower_x(c_x);
- r -= robust_fpt_type(line_a, false) * robust_fpt_type(site3.x0(), false);
- r -= robust_fpt_type(line_b, false) * robust_fpt_type(site3.y0(), false);
- r += robust_fpt_type(line_a, false) * c_x;
- r += robust_fpt_type(line_b, false) * c_y;
- if (r.pos().fpv() < r.neg().fpv()) {
- r = -r;
- }
- lower_x += r * inv_segm_len;
- c_event = circle_type(c_x.dif().fpv(), c_y.dif().fpv(), lower_x.dif().fpv());
- bool recompute_c_x = c_x.dif().ulp() > ULPS;
- bool recompute_c_y = c_y.dif().ulp() > ULPS;
- bool recompute_lower_x = lower_x.dif().ulp() > ULPS;
- if (recompute_c_x || recompute_c_y || recompute_lower_x) {
- exact_circle_formation_functor_.pps(
- site1, site2, site3, segment_index, c_event,
- recompute_c_x, recompute_c_y, recompute_lower_x);
- }
+ recompute_c_x = c_x.dif().ulp() > ULPS;
+ recompute_c_y = c_y.dif().ulp() > ULPS;
+ recompute_lower_x = lower_x.dif().ulp() > ULPS;
+ c_event =
+ circle_type(c_x.dif().fpv(), c_y.dif().fpv(), lower_x.dif().fpv());
+ } else {
+ robust_fpt_type sqr_sum1(get_sqrt(a1 * a1 + b1 * b1), 2.0);
+ robust_fpt_type sqr_sum2(get_sqrt(a2 * a2 + b2 * b2), 2.0);
+ robust_fpt_type a(robust_cross_product(a1, b1, -b2, a2), 1.0);
+ if (!is_neg(a)) {
+ a += sqr_sum1 * sqr_sum2;
+ } else {
+ a = (orientation * orientation) / (sqr_sum1 * sqr_sum2 - a);
         }
-
- void pss(const site_type &site1,
- const site_type &site2,
- const site_type &site3,
- int point_index,
- circle_type &c_event) {
- const point_type &segm_start1 = site2.point1(true);
- const point_type &segm_end1 = site2.point0(true);
- const point_type &segm_start2 = site3.point0(true);
- const point_type &segm_end2 = site3.point1(true);
- fpt_type a1 = to_fpt(segm_end1.x()) - to_fpt(segm_start1.x());
- fpt_type b1 = to_fpt(segm_end1.y()) - to_fpt(segm_start1.y());
- fpt_type a2 = to_fpt(segm_end2.x()) - to_fpt(segm_start2.x());
- fpt_type b2 = to_fpt(segm_end2.y()) - to_fpt(segm_start2.y());
- bool recompute_c_x, recompute_c_y, recompute_lower_x;
- robust_fpt_type orientation(robust_cross_product(b1, a1, b2, a2), 1.0);
- if (ot::eval(orientation) == ot::COLLINEAR) {
- robust_fpt_type a(a1 * a1 + b1 * b1, 2.0);
- robust_fpt_type c(robust_cross_product(
- b1, a1,
- to_fpt(segm_start2.y()) - to_fpt(segm_start1.y()),
- to_fpt(segm_start2.x()) - to_fpt(segm_start1.x())), 1.0);
- robust_fpt_type det(
- robust_cross_product(
- a1, b1,
- to_fpt(site1.x()) - to_fpt(segm_start1.x()),
- to_fpt(site1.y()) - to_fpt(segm_start1.y())) *
- robust_cross_product(
- b1, a1,
- to_fpt(site1.y()) - to_fpt(segm_start2.y()),
- to_fpt(site1.x()) - to_fpt(segm_start2.x())),
- 3.0);
- robust_dif_type t;
- t -= robust_fpt_type(a1, false) * robust_fpt_type(
- (to_fpt(segm_start1.x()) + to_fpt(segm_start2.x())) * to_fpt(0.5) -
- to_fpt(site1.x()), false);
- t -= robust_fpt_type(b1, false) * robust_fpt_type((
- to_fpt(segm_start1.y()) + to_fpt(segm_start2.y())) * to_fpt(0.5) -
- to_fpt(site1.y()), false);
- if (point_index == 2) {
- t += det.sqrt();
- } else {
- t -= det.sqrt();
- }
- t /= a;
- robust_dif_type c_x, c_y;
- c_x += robust_fpt_type(to_fpt(0.5) * (
- to_fpt(segm_start1.x()) + to_fpt(segm_start2.x())), false);
- c_x += robust_fpt_type(a1, false) * t;
- c_y += robust_fpt_type(to_fpt(0.5) * (
- to_fpt(segm_start1.y()) + to_fpt(segm_start2.y())), false);
- c_y += robust_fpt_type(b1, false) * t;
- robust_dif_type lower_x(c_x);
- if (is_neg(c)) {
- lower_x -= robust_fpt_type(0.5, false) * c / a.sqrt();
- } else {
- lower_x += robust_fpt_type(0.5, false) * c / a.sqrt();
- }
- recompute_c_x = c_x.dif().ulp() > ULPS;
- recompute_c_y = c_y.dif().ulp() > ULPS;
- recompute_lower_x = lower_x.dif().ulp() > ULPS;
- c_event = circle_type(c_x.dif().fpv(), c_y.dif().fpv(), lower_x.dif().fpv());
- } else {
- robust_fpt_type sqr_sum1(get_sqrt(a1 * a1 + b1 * b1), 2.0);
- robust_fpt_type sqr_sum2(get_sqrt(a2 * a2 + b2 * b2), 2.0);
- robust_fpt_type a(robust_cross_product(a1, b1, -b2, a2), 1.0);
- if (!is_neg(a)) {
- a += sqr_sum1 * sqr_sum2;
- } else {
- a = (orientation * orientation) / (sqr_sum1 * sqr_sum2 - a);
- }
- robust_fpt_type or1(robust_cross_product(
- b1, a1,
- to_fpt(segm_end1.y()) - to_fpt(site1.y()),
- to_fpt(segm_end1.x()) - to_fpt(site1.x())), 1.0);
- robust_fpt_type or2(robust_cross_product(
- a2, b2,
- to_fpt(segm_end2.x()) - to_fpt(site1.x()),
- to_fpt(segm_end2.y()) - to_fpt(site1.y())), 1.0);
- robust_fpt_type det = robust_fpt_type(2.0, false) * a * or1 * or2;
- robust_fpt_type c1(robust_cross_product(
- b1, a1,
- to_fpt(segm_end1.y()), to_fpt(segm_end1.x())), 1.0);
- robust_fpt_type c2(robust_cross_product(
- a2, b2,
- to_fpt(segm_end2.x()), to_fpt(segm_end2.y())), 1.0);
- robust_fpt_type inv_orientation = robust_fpt_type(1.0, false) / orientation;
- robust_dif_type t, b, ix, iy;
- ix += robust_fpt_type(a2, false) * c1 * inv_orientation;
- ix += robust_fpt_type(a1, false) * c2 * inv_orientation;
- iy += robust_fpt_type(b1, false) * c2 * inv_orientation;
- iy += robust_fpt_type(b2, false) * c1 * inv_orientation;
-
- b += ix * (robust_fpt_type(a1, false) * sqr_sum2);
- b += ix * (robust_fpt_type(a2, false) * sqr_sum1);
- b += iy * (robust_fpt_type(b1, false) * sqr_sum2);
- b += iy * (robust_fpt_type(b2, false) * sqr_sum1);
- b -= sqr_sum1 * robust_fpt_type(robust_cross_product(
- a2, b2,
- to_fpt(-site1.y()), to_fpt(site1.x())), 1.0);
- b -= sqr_sum2 * robust_fpt_type(robust_cross_product(
- a1, b1,
- to_fpt(-site1.y()), to_fpt(site1.x())), 1.0);
- t -= b;
- if (point_index == 2) {
- t += det.sqrt();
- } else {
- t -= det.sqrt();
- }
- t /= (a * a);
- robust_dif_type c_x(ix), c_y(iy);
- c_x += t * (robust_fpt_type(a1, false) * sqr_sum2);
- c_x += t * (robust_fpt_type(a2, false) * sqr_sum1);
- c_y += t * (robust_fpt_type(b1, false) * sqr_sum2);
- c_y += t * (robust_fpt_type(b2, false) * sqr_sum1);
- if (t.pos().fpv() < t.neg().fpv()) {
- t = -t;
- }
- robust_dif_type lower_x(c_x);
- if (is_neg(orientation)) {
- lower_x -= t * orientation;
- } else {
- lower_x += t * orientation;
- }
- recompute_c_x = c_x.dif().ulp() > ULPS;
- recompute_c_y = c_y.dif().ulp() > ULPS;
- recompute_lower_x = lower_x.dif().ulp() > ULPS;
- c_event = circle_type(c_x.dif().fpv(), c_y.dif().fpv(), lower_x.dif().fpv());
- }
- if (recompute_c_x || recompute_c_y || recompute_lower_x) {
- exact_circle_formation_functor_.pss(
- site1, site2, site3, point_index, c_event,
- recompute_c_x, recompute_c_y, recompute_lower_x);
- }
+ robust_fpt_type or1(robust_cross_product(
+ b1, a1,
+ to_fpt(segm_end1.y()) - to_fpt(site1.y()),
+ to_fpt(segm_end1.x()) - to_fpt(site1.x())), 1.0);
+ robust_fpt_type or2(robust_cross_product(
+ a2, b2,
+ to_fpt(segm_end2.x()) - to_fpt(site1.x()),
+ to_fpt(segm_end2.y()) - to_fpt(site1.y())), 1.0);
+ robust_fpt_type det = robust_fpt_type(2.0, false) * a * or1 * or2;
+ robust_fpt_type c1(robust_cross_product(
+ b1, a1, to_fpt(segm_end1.y()), to_fpt(segm_end1.x())), 1.0);
+ robust_fpt_type c2(robust_cross_product(
+ a2, b2, to_fpt(segm_end2.x()), to_fpt(segm_end2.y())), 1.0);
+ robust_fpt_type inv_orientation =
+ robust_fpt_type(1.0, false) / orientation;
+ robust_dif_type t, b, ix, iy;
+ ix += robust_fpt_type(a2, false) * c1 * inv_orientation;
+ ix += robust_fpt_type(a1, false) * c2 * inv_orientation;
+ iy += robust_fpt_type(b1, false) * c2 * inv_orientation;
+ iy += robust_fpt_type(b2, false) * c1 * inv_orientation;
+
+ b += ix * (robust_fpt_type(a1, false) * sqr_sum2);
+ b += ix * (robust_fpt_type(a2, false) * sqr_sum1);
+ b += iy * (robust_fpt_type(b1, false) * sqr_sum2);
+ b += iy * (robust_fpt_type(b2, false) * sqr_sum1);
+ b -= sqr_sum1 * robust_fpt_type(robust_cross_product(
+ a2, b2, to_fpt(-site1.y()), to_fpt(site1.x())), 1.0);
+ b -= sqr_sum2 * robust_fpt_type(robust_cross_product(
+ a1, b1, to_fpt(-site1.y()), to_fpt(site1.x())), 1.0);
+ t -= b;
+ if (point_index == 2) {
+ t += det.sqrt();
+ } else {
+ t -= det.sqrt();
         }
-
- void sss(const site_type &site1,
- const site_type &site2,
- const site_type &site3,
- circle_type &c_event) {
- robust_fpt_type a1(to_fpt(site1.x1(true)) - to_fpt(site1.x0(true)), 0.0);
- robust_fpt_type b1(to_fpt(site1.y1(true)) - to_fpt(site1.y0(true)), 0.0);
- robust_fpt_type c1(robust_cross_product(
- site1.x0(true), site1.y0(true), site1.x1(true), site1.y1(true)), 1.0);
-
- robust_fpt_type a2(to_fpt(site2.x1(true)) - to_fpt(site2.x0(true)), 0.0);
- robust_fpt_type b2(to_fpt(site2.y1(true)) - to_fpt(site2.y0(true)), 0.0);
- robust_fpt_type c2(robust_cross_product(
- site2.x0(true), site2.y0(true), site2.x1(true), site2.y1(true)), 1.0);
-
- robust_fpt_type a3(to_fpt(site3.x1(true)) - to_fpt(site3.x0(true)), 0.0);
- robust_fpt_type b3(to_fpt(site3.y1(true)) - to_fpt(site3.y0(true)), 0.0);
- robust_fpt_type c3(robust_cross_product(
- site3.x0(true), site3.y0(true), site3.x1(true), site3.y1(true)), 1.0);
-
- robust_fpt_type len1 = (a1 * a1 + b1 * b1).sqrt();
- robust_fpt_type len2 = (a2 * a2 + b2 * b2).sqrt();
- robust_fpt_type len3 = (a3 * a3 + b3 * b3).sqrt();
- robust_fpt_type cross_12(robust_cross_product(
- a1.fpv(), b1.fpv(), a2.fpv(), b2.fpv()), 1.0);
- robust_fpt_type cross_23(robust_cross_product(
- a2.fpv(), b2.fpv(), a3.fpv(), b3.fpv()), 1.0);
- robust_fpt_type cross_31(robust_cross_product(
- a3.fpv(), b3.fpv(), a1.fpv(), b1.fpv()), 1.0);
- robust_dif_type denom, c_x, c_y, r;
-
- // denom = cross_12 * len3 + cross_23 * len1 + cross_31 * len2.
- denom += cross_12 * len3;
- denom += cross_23 * len1;
- denom += cross_31 * len2;
-
- // denom * r = (b2 * c_x - a2 * c_y - c2 * denom) / len2.
- r -= cross_12 * c3;
- r -= cross_23 * c1;
- r -= cross_31 * c2;
-
- c_x += a1 * c2 * len3;
- c_x -= a2 * c1 * len3;
- c_x += a2 * c3 * len1;
- c_x -= a3 * c2 * len1;
- c_x += a3 * c1 * len2;
- c_x -= a1 * c3 * len2;
- c_y += b1 * c2 * len3;
- c_y -= b2 * c1 * len3;
- c_y += b2 * c3 * len1;
- c_y -= b3 * c2 * len1;
- c_y += b3 * c1 * len2;
- c_y -= b1 * c3 * len2;
- robust_dif_type lower_x(c_x + r);
- bool recompute_c_x = c_x.dif().ulp() > ULPS;
- bool recompute_c_y = c_y.dif().ulp() > ULPS;
- bool recompute_lower_x = lower_x.dif().ulp() > ULPS;
- bool recompute_denom = denom.dif().ulp() > ULPS;
- c_event = circle_type(c_x.dif().fpv() / denom.dif().fpv(),
- c_y.dif().fpv() / denom.dif().fpv(),
- lower_x.dif().fpv() / denom.dif().fpv());
- if (recompute_c_x || recompute_c_y || recompute_lower_x || recompute_denom) {
- exact_circle_formation_functor_.sss(
- site1, site2, site3, c_event,
- recompute_c_x, recompute_c_y, recompute_lower_x);
- }
+ t /= (a * a);
+ robust_dif_type c_x(ix), c_y(iy);
+ c_x += t * (robust_fpt_type(a1, false) * sqr_sum2);
+ c_x += t * (robust_fpt_type(a2, false) * sqr_sum1);
+ c_y += t * (robust_fpt_type(b1, false) * sqr_sum2);
+ c_y += t * (robust_fpt_type(b2, false) * sqr_sum1);
+ if (t.pos().fpv() < t.neg().fpv()) {
+ t = -t;
+ }
+ robust_dif_type lower_x(c_x);
+ if (is_neg(orientation)) {
+ lower_x -= t * orientation;
+ } else {
+ lower_x += t * orientation;
         }
+ recompute_c_x = c_x.dif().ulp() > ULPS;
+ recompute_c_y = c_y.dif().ulp() > ULPS;
+ recompute_lower_x = lower_x.dif().ulp() > ULPS;
+ c_event = circle_type(
+ c_x.dif().fpv(), c_y.dif().fpv(), lower_x.dif().fpv());
+ }
+ if (recompute_c_x || recompute_c_y || recompute_lower_x) {
+ exact_circle_formation_functor_.pss(
+ site1, site2, site3, point_index, c_event,
+ recompute_c_x, recompute_c_y, recompute_lower_x);
+ }
+ }

- private:
- exact_circle_formation_functor_type exact_circle_formation_functor_;
- to_fpt_converter to_fpt;
- };
+ void sss(const site_type &site1,
+ const site_type &site2,
+ const site_type &site3,
+ circle_type &c_event) {
+ robust_fpt_type a1(to_fpt(site1.x1(true)) - to_fpt(site1.x0(true)), 0.0);
+ robust_fpt_type b1(to_fpt(site1.y1(true)) - to_fpt(site1.y0(true)), 0.0);
+ robust_fpt_type c1(robust_cross_product(
+ site1.x0(true), site1.y0(true),
+ site1.x1(true), site1.y1(true)), 1.0);
+
+ robust_fpt_type a2(to_fpt(site2.x1(true)) - to_fpt(site2.x0(true)), 0.0);
+ robust_fpt_type b2(to_fpt(site2.y1(true)) - to_fpt(site2.y0(true)), 0.0);
+ robust_fpt_type c2(robust_cross_product(
+ site2.x0(true), site2.y0(true),
+ site2.x1(true), site2.y1(true)), 1.0);
+
+ robust_fpt_type a3(to_fpt(site3.x1(true)) - to_fpt(site3.x0(true)), 0.0);
+ robust_fpt_type b3(to_fpt(site3.y1(true)) - to_fpt(site3.y0(true)), 0.0);
+ robust_fpt_type c3(robust_cross_product(
+ site3.x0(true), site3.y0(true),
+ site3.x1(true), site3.y1(true)), 1.0);
+
+ robust_fpt_type len1 = (a1 * a1 + b1 * b1).sqrt();
+ robust_fpt_type len2 = (a2 * a2 + b2 * b2).sqrt();
+ robust_fpt_type len3 = (a3 * a3 + b3 * b3).sqrt();
+ robust_fpt_type cross_12(robust_cross_product(
+ a1.fpv(), b1.fpv(), a2.fpv(), b2.fpv()), 1.0);
+ robust_fpt_type cross_23(robust_cross_product(
+ a2.fpv(), b2.fpv(), a3.fpv(), b3.fpv()), 1.0);
+ robust_fpt_type cross_31(robust_cross_product(
+ a3.fpv(), b3.fpv(), a1.fpv(), b1.fpv()), 1.0);
+ robust_dif_type denom, c_x, c_y, r;
+
+ // denom = cross_12 * len3 + cross_23 * len1 + cross_31 * len2.
+ denom += cross_12 * len3;
+ denom += cross_23 * len1;
+ denom += cross_31 * len2;
+
+ // denom * r = (b2 * c_x - a2 * c_y - c2 * denom) / len2.
+ r -= cross_12 * c3;
+ r -= cross_23 * c1;
+ r -= cross_31 * c2;
+
+ c_x += a1 * c2 * len3;
+ c_x -= a2 * c1 * len3;
+ c_x += a2 * c3 * len1;
+ c_x -= a3 * c2 * len1;
+ c_x += a3 * c1 * len2;
+ c_x -= a1 * c3 * len2;
+ c_y += b1 * c2 * len3;
+ c_y -= b2 * c1 * len3;
+ c_y += b2 * c3 * len1;
+ c_y -= b3 * c2 * len1;
+ c_y += b3 * c1 * len2;
+ c_y -= b1 * c3 * len2;
+ robust_dif_type lower_x(c_x + r);
+ bool recompute_c_x = c_x.dif().ulp() > ULPS;
+ bool recompute_c_y = c_y.dif().ulp() > ULPS;
+ bool recompute_lower_x = lower_x.dif().ulp() > ULPS;
+ bool recompute_denom = denom.dif().ulp() > ULPS;
+ c_event = circle_type(
+ c_x.dif().fpv() / denom.dif().fpv(),
+ c_y.dif().fpv() / denom.dif().fpv(),
+ lower_x.dif().fpv() / denom.dif().fpv());
+ if (recompute_c_x || recompute_c_y ||
+ recompute_lower_x || recompute_denom) {
+ exact_circle_formation_functor_.sss(
+ site1, site2, site3, c_event,
+ recompute_c_x, recompute_c_y, recompute_lower_x);
+ }
+ }

- template <typename Site,
- typename Circle,
- typename CEP = circle_existence_predicate<Site>,
- typename CFF = lazy_circle_formation_functor<Site, Circle> >
- class circle_formation_predicate {
- public:
- typedef Site site_type;
- typedef Circle circle_type;
- typedef CEP circle_existence_predicate_type;
- typedef CFF circle_formation_functor_type;
-
- // Create a circle event from the given three sites.
- // Returns true if the circle event exists, else false.
- // If exists circle event is saved into the c_event variable.
- bool operator()(const site_type &site1, const site_type &site2,
- const site_type &site3, circle_type &circle) {
- if (!site1.is_segment()) {
- if (!site2.is_segment()) {
- if (!site3.is_segment()) {
- // (point, point, point) sites.
- if (!circle_existence_predicate_.ppp(site1, site2, site3))
- return false;
- circle_formation_functor_.ppp(site1, site2, site3, circle);
- } else {
- // (point, point, segment) sites.
- if (!circle_existence_predicate_.pps(site1, site2, site3, 3))
- return false;
- circle_formation_functor_.pps(site1, site2, site3, 3, circle);
- }
- } else {
- if (!site3.is_segment()) {
- // (point, segment, point) sites.
- if (!circle_existence_predicate_.pps(site1, site3, site2, 2))
- return false;
- circle_formation_functor_.pps(site1, site3, site2, 2, circle);
- } else {
- // (point, segment, segment) sites.
- if (!circle_existence_predicate_.pss(site1, site2, site3, 1))
- return false;
- circle_formation_functor_.pss(site1, site2, site3, 1, circle);
- }
- }
- } else {
- if (!site2.is_segment()) {
- if (!site3.is_segment()) {
- // (segment, point, point) sites.
- if (!circle_existence_predicate_.pps(site2, site3, site1, 1))
- return false;
- circle_formation_functor_.pps(site2, site3, site1, 1, circle);
- } else {
- // (segment, point, segment) sites.
- if (!circle_existence_predicate_.pss(site2, site1, site3, 2))
- return false;
- circle_formation_functor_.pss(site2, site1, site3, 2, circle);
- }
- } else {
- if (!site3.is_segment()) {
- // (segment, segment, point) sites.
- if (!circle_existence_predicate_.pss(site3, site1, site2, 3))
- return false;
- circle_formation_functor_.pss(site3, site1, site2, 3, circle);
- } else {
- // (segment, segment, segment) sites.
- if (!circle_existence_predicate_.sss(site1, site2, site3))
- return false;
- circle_formation_functor_.sss(site1, site2, site3, circle);
- }
- }
- }
- return true;
+ private:
+ exact_circle_formation_functor_type exact_circle_formation_functor_;
+ to_fpt_converter to_fpt;
+ };
+
+ template <typename Site,
+ typename Circle,
+ typename CEP = circle_existence_predicate<Site>,
+ typename CFF = lazy_circle_formation_functor<Site, Circle> >
+ class circle_formation_predicate {
+ public:
+ typedef Site site_type;
+ typedef Circle circle_type;
+ typedef CEP circle_existence_predicate_type;
+ typedef CFF circle_formation_functor_type;
+
+ // Create a circle event from the given three sites.
+ // Returns true if the circle event exists, else false.
+ // If exists circle event is saved into the c_event variable.
+ bool operator()(const site_type &site1, const site_type &site2,
+ const site_type &site3, circle_type &circle) {
+ if (!site1.is_segment()) {
+ if (!site2.is_segment()) {
+ if (!site3.is_segment()) {
+ // (point, point, point) sites.
+ if (!circle_existence_predicate_.ppp(site1, site2, site3))
+ return false;
+ circle_formation_functor_.ppp(site1, site2, site3, circle);
+ } else {
+ // (point, point, segment) sites.
+ if (!circle_existence_predicate_.pps(site1, site2, site3, 3))
+ return false;
+ circle_formation_functor_.pps(site1, site2, site3, 3, circle);
+ }
+ } else {
+ if (!site3.is_segment()) {
+ // (point, segment, point) sites.
+ if (!circle_existence_predicate_.pps(site1, site3, site2, 2))
+ return false;
+ circle_formation_functor_.pps(site1, site3, site2, 2, circle);
+ } else {
+ // (point, segment, segment) sites.
+ if (!circle_existence_predicate_.pss(site1, site2, site3, 1))
+ return false;
+ circle_formation_functor_.pss(site1, site2, site3, 1, circle);
+ }
+ }
+ } else {
+ if (!site2.is_segment()) {
+ if (!site3.is_segment()) {
+ // (segment, point, point) sites.
+ if (!circle_existence_predicate_.pps(site2, site3, site1, 1))
+ return false;
+ circle_formation_functor_.pps(site2, site3, site1, 1, circle);
+ } else {
+ // (segment, point, segment) sites.
+ if (!circle_existence_predicate_.pss(site2, site1, site3, 2))
+ return false;
+ circle_formation_functor_.pss(site2, site1, site3, 2, circle);
+ }
+ } else {
+ if (!site3.is_segment()) {
+ // (segment, segment, point) sites.
+ if (!circle_existence_predicate_.pss(site3, site1, site2, 3))
+ return false;
+ circle_formation_functor_.pss(site3, site1, site2, 3, circle);
+ } else {
+ // (segment, segment, segment) sites.
+ if (!circle_existence_predicate_.sss(site1, site2, site3))
+ return false;
+ circle_formation_functor_.sss(site1, site2, site3, circle);
+ }
         }
+ }
+ return true;
+ }

- private:
- circle_existence_predicate_type circle_existence_predicate_;
- circle_formation_functor_type circle_formation_functor_;
- };
+ private:
+ circle_existence_predicate_type circle_existence_predicate_;
+ circle_formation_functor_type circle_formation_functor_;
+ };
};
} // detail
} // polygon

Modified: sandbox/gtl/boost/polygon/detail/voronoi_robust_fpt.hpp
==============================================================================
--- sandbox/gtl/boost/polygon/detail/voronoi_robust_fpt.hpp (original)
+++ sandbox/gtl/boost/polygon/detail/voronoi_robust_fpt.hpp 2012-03-25 15:08:13 EDT (Sun, 25 Mar 2012)
@@ -153,7 +153,7 @@
   robust_fpt& operator*=(const robust_fpt &that) {
     this->re_ += that.re_ + ROUNDING_ERROR;
     this->fpv_ *= that.fpv_;
- return *this;
+ return *this;
   }

   robust_fpt& operator/=(const robust_fpt &that) {

Next message: vicente.botet_at_[hidden]: "[Boost-commit] svn:boost r77545 - trunk/boost/chrono/detail/inlined/win"
Previous message: vicente.botet_at_[hidden]: "[Boost-commit] svn:boost r77543 - in trunk/libs/utility/doc: . html"

Date view	Thread view	Subject view	Author view

Boost-Commit list run by bdawes at acm.org, david.abrahams at rcn.com, gregod at cs.rpi.edu, cpdaniel at pacbell.net, john at johnmaddock.co.uk