Subject: [Boost-commit] svn:boost r74942 - in sandbox/gtl: boost/polygon boost/polygon/detail libs/polygon/test
From: sydorchuk.andriy_at_[hidden]
Date: 2011-10-13 19:51:47
Author: asydorchuk
Date: 2011-10-13 19:51:44 EDT (Thu, 13 Oct 2011)
New Revision: 74942
URL: http://svn.boost.org/trac/boost/changeset/74942
Log:
Integrated voronoi_fpt_kernel_tests.
Changed erc to robust_dif.
Added robust_dif tests.
Refactored sqrt_expr_evaluator(robust_sqrt_expr) - won't require static fields in the future.
Fixed voronoi_clipping_test warnings on Linux.
Added:
sandbox/gtl/libs/polygon/test/voronoi_fpt_kernel_test.cpp (contents, props changed)
Text files modified:
sandbox/gtl/boost/polygon/detail/voronoi_calc_kernel.hpp | 74 ++++---
sandbox/gtl/boost/polygon/detail/voronoi_fpt_kernel.hpp | 386 +++++++++++++++++++++------------------
sandbox/gtl/boost/polygon/voronoi_builder.hpp | 4
sandbox/gtl/libs/polygon/test/Jamfile.v2 | 1
sandbox/gtl/libs/polygon/test/voronoi_clipping_test.cpp | 44 ++--
5 files changed, 275 insertions(+), 234 deletions(-)
Modified: sandbox/gtl/boost/polygon/detail/voronoi_calc_kernel.hpp
==============================================================================
--- sandbox/gtl/boost/polygon/detail/voronoi_calc_kernel.hpp (original)
+++ sandbox/gtl/boost/polygon/detail/voronoi_calc_kernel.hpp 2011-10-13 19:51:44 EDT (Thu, 13 Oct 2011)
@@ -10,6 +10,10 @@
#ifndef BOOST_POLYGON_VORONOI_CALC_KERNEL
#define BOOST_POLYGON_VORONOI_CALC_KERNEL
+#pragma warning(disable:4800)
+#include <gmpxx.h>
+
+#include "mpt_wrapper.hpp"
#include "voronoi_fpt_kernel.hpp"
namespace boost {
@@ -61,7 +65,7 @@
class voronoi_calc_kernel<int> {
public:
typedef double fpt_type;
- typedef polygon_ulong_long_type ulong_long_type;
+ typedef unsigned long long ulong_long_type;
enum kOrientation {
RIGHT = -1,
@@ -543,6 +547,10 @@
typedef Site site_type;
typedef Circle circle_type;
+ typedef mpt_wrapper<mpz_class, 8> mpt_type;
+ typedef mpt_wrapper<mpf_class, 2> mpf_type;
+ typedef robust_sqrt_expr<mpt_type, mpf_type> robust_sqrt_expr_type;
+
bool ppp(const site_type &site1,
const site_type &site2,
const site_type &site3,
@@ -649,7 +657,7 @@
c_event.y(0.25 * cA[2].get_d() / denom.get_d());
}
if (recompute_lower_x) {
- c_event.lower_x(0.25 * sqrt_expr_evaluator<2>::eval<mpt_type, mpf_type>(cA, cB).get_d() /
+ c_event.lower_x(0.25 * sqrt_expr_.eval2(cA, cB).get_d() /
(denom.get_d() * std::sqrt(segm_len.get_d())));
}
return true;
@@ -664,8 +672,7 @@
cA[1] = (segment_index == 2) ? -vec_x : vec_x;
cB[1] = det;
if (recompute_c_x) {
- c_event.x(0.5 * sqrt_expr_evaluator<2>::eval<mpt_type, mpf_type>(cA, cB).get_d() /
- denom_sqr);
+ c_event.x(0.5 * sqrt_expr_.eval2(cA, cB).get_d() / denom_sqr);
}
}
@@ -675,7 +682,7 @@
cA[3] = (segment_index == 2) ? -vec_y : vec_y;
cB[3] = det;
if (recompute_c_y) {
- c_event.y(0.5 * sqrt_expr_evaluator<2>::eval<mpt_type, mpf_type>(&cA[2], &cB[2]).get_d() /
+ c_event.y(0.5 * sqrt_expr_.eval2(&cA[2], &cB[2]).get_d() /
denom_sqr);
}
}
@@ -687,7 +694,7 @@
cB[2] = 1;
cA[3] = (segment_index == 2) ? -teta : teta;
cB[3] = det;
- c_event.lower_x(0.5 * sqrt_expr_evaluator<4>::eval<mpt_type, mpf_type>(cA, cB).get_d() /
+ c_event.lower_x(0.5 * sqrt_expr_.eval4(cA, cB).get_d() /
(denom_sqr * std::sqrt(segm_len.get_d())));
}
@@ -733,7 +740,7 @@
cA[1] = a[0] * a[0] * (segm_start1.y() + segm_start2.y()) -
a[0] * b[0] * (segm_start1.x() + segm_start2.x() - 2.0 * site1.x()) +
b[0] * b[0] * (2.0 * site1.y());
- double c_y = sqrt_expr_evaluator<2>::eval<mpt_type, mpf_type>(cA, cB).get_d();
+ double c_y = sqrt_expr_.eval2(cA, cB).get_d();
c_event.y(c_y / denom);
}
@@ -744,14 +751,14 @@
a[0] * a[0] * (2.0 * site1.x());
if (recompute_c_x) {
- double c_x = sqrt_expr_evaluator<2>::eval<mpt_type, mpf_type>(cA, cB).get_d();
+ double c_x = sqrt_expr_.eval2(cA, cB).get_d();
c_event.x(c_x / denom);
}
if (recompute_lower_x) {
cA[2] = (c[0] < 0) ? -c[0] : c[0];
cB[2] = a[0] * a[0] + b[0] * b[0];
- double lower_x = sqrt_expr_evaluator<3>::eval<mpt_type, mpf_type>(cA, cB).get_d();
+ double lower_x = sqrt_expr_.eval3(cA, cB).get_d();
c_event.lower_x(lower_x / denom);
}
}
@@ -844,7 +851,7 @@
int k = (i+2) % 3;
cross[i] = a[j] * b[k] - a[k] * b[j];
}
- double denom = sqrt_expr_evaluator<3>::eval<mpt_type, mpf_type>(cross, sqr_len).get_d();
+ double denom = sqrt_expr_.eval3(cross, sqr_len).get_d();
if (recompute_c_y) {
for (int i = 0; i < 3; ++i) {
@@ -852,7 +859,7 @@
int k = (i+2) % 3;
cross[i] = b[j] * c[k] - b[k] * c[j];
}
- double c_y = sqrt_expr_evaluator<3>::eval<mpt_type, mpf_type>(cross, sqr_len).get_d();
+ double c_y = sqrt_expr_.eval3(cross, sqr_len).get_d();
c_event.y(c_y / denom);
}
@@ -868,13 +875,13 @@
}
if (recompute_c_x) {
- double c_x = sqrt_expr_evaluator<3>::eval<mpt_type, mpf_type>(cross, sqr_len).get_d();
+ double c_x = sqrt_expr_.eval3(cross, sqr_len).get_d();
c_event.x(c_x / denom);
}
if (recompute_lower_x) {
sqr_len[3] = 1;
- double lower_x = sqrt_expr_evaluator<4>::eval<mpt_type, mpf_type>(cross, sqr_len).get_d();
+ double lower_x = sqrt_expr_.eval4(cross, sqr_len).get_d();
c_event.lower_x(lower_x / denom);
}
}
@@ -901,13 +908,13 @@
static mpt cA[4], cB[4];
static mpf lh, rh, numer;
if (A[3] == 0) {
- lh = sqrt_expr_evaluator<2>::eval<mpt, mpf>(A, B);
+ lh = sqrt_expr_.eval2(A, B);
cA[0] = 1;
cB[0] = B[0] * B[1];
cA[1] = B[2];
cB[1] = 1;
rh = A[2].get_d() * std::sqrt(B[3].get_d() *
- sqrt_expr_evaluator<2>::eval<mpt, mpf>(cA, cB).get_d());
+ sqrt_expr_.eval2(cA, cB).get_d());
if (((lh >= 0) && (rh >= 0)) || ((lh <= 0) && (rh <= 0))) {
return lh + rh;
}
@@ -916,7 +923,7 @@
cB[0] = 1;
cA[1] = A[0] * A[1] * 2 - A[2] * A[2] * B[3];
cB[1] = B[0] * B[1];
- numer = sqrt_expr_evaluator<2>::eval<mpt, mpf>(cA, cB);
+ numer = sqrt_expr_.eval2(cA, cB);
return numer / (lh - rh);
}
cA[0] = A[3];
@@ -925,13 +932,13 @@
cB[1] = B[0];
cA[2] = A[1];
cB[2] = B[1];
- lh = sqrt_expr_evaluator<3>::eval<mpt, mpf>(cA, cB);
+ lh = sqrt_expr_.eval3(cA, cB);
cA[0] = 1;
cB[0] = B[0] * B[1];
cA[1] = B[2];
cB[1] = 1;
rh = A[2].get_d() * std::sqrt(B[3].get_d() *
- sqrt_expr_evaluator<2>::eval<mpt, mpf>(cA, cB).get_d());
+ sqrt_expr_.eval2(cA, cB).get_d());
if (((lh >= 0) && (rh >= 0)) || ((lh <= 0) && (rh <= 0))) {
return lh + rh;
}
@@ -944,9 +951,12 @@
cB[2] = B[1];
cA[3] = A[0] * A[1] * 2 - A[2] * A[2] * B[3];
cB[3] = B[0] * B[1];
- numer = sqrt_expr_evaluator<4>::eval<mpt, mpf>(cA, cB);
+ numer = sqrt_expr_.eval4(cA, cB);
return numer / (lh - rh);
}
+
+ private:
+ robust_sqrt_expr_type sqrt_expr_;
};
template <typename Site, typename Circle>
@@ -976,7 +986,7 @@
double sum_y2 = site2.y() + site3.y();
double dif_x3 = site1.x() - site3.x();
double dif_y3 = site1.y() - site3.y();
- epsilon_robust_comparator< robust_fpt<double> > c_x, c_y;
+ robust_dif< robust_fpt<double> > c_x, c_y;
c_x += robust_fpt<double>(dif_x1 * sum_x1 * dif_y2, 2.0);
c_x += robust_fpt<double>(dif_y1 * sum_y1 * dif_y2, 2.0);
c_x -= robust_fpt<double>(dif_x2 * sum_x2 * dif_y1, 2.0);
@@ -985,7 +995,7 @@
c_y += robust_fpt<double>(dif_y2 * sum_y2 * dif_x1, 2.0);
c_y -= robust_fpt<double>(dif_x1 * sum_x1 * dif_x2, 2.0);
c_y -= robust_fpt<double>(dif_y1 * sum_y1 * dif_x2, 2.0);
- epsilon_robust_comparator< robust_fpt<double> > lower_x(c_x);
+ robust_dif< robust_fpt<double> > lower_x(c_x);
lower_x -= robust_fpt<double>(std::sqrt(sqr_distance(dif_x1, dif_y1) *
sqr_distance(dif_x2, dif_y2) *
sqr_distance(dif_x3, dif_y3)), 5.0);
@@ -1022,7 +1032,7 @@
site2.x() - site3.point1().x()), 2.0);
robust_fpt<double> denom(robust_cross_product(vec_x, vec_y, line_a, line_b), 2.0);
robust_fpt<double> inv_segm_len(1.0 / std::sqrt(sqr_distance(line_a, line_b)), 3.0);
- epsilon_robust_comparator< robust_fpt<double> > t;
+ robust_dif< robust_fpt<double> > t;
if (get_orientation(denom) == COLLINEAR) {
t += teta / (robust_fpt<double>(8.0) * A);
t -= A / (robust_fpt<double>(2.0) * teta);
@@ -1035,12 +1045,12 @@
}
t += teta * (A + B) / (robust_fpt<double>(2.0) * denom * denom);
}
- epsilon_robust_comparator< robust_fpt<double> > c_x, c_y;
+ robust_dif< robust_fpt<double> > c_x, c_y;
c_x += robust_fpt<double>(0.5 * (site1.x() + site2.x()));
c_x += robust_fpt<double>(vec_x) * t;
c_y += robust_fpt<double>(0.5 * (site1.y() + site2.y()));
c_y += robust_fpt<double>(vec_y) * t;
- epsilon_robust_comparator< robust_fpt<double> > r, lower_x(c_x);
+ robust_dif< robust_fpt<double> > r, lower_x(c_x);
r -= robust_fpt<double>(line_a) * robust_fpt<double>(site3.x0());
r -= robust_fpt<double>(line_b) * robust_fpt<double>(site3.y0());
r += robust_fpt<double>(line_a) * c_x;
@@ -1084,7 +1094,7 @@
site1.y() - segm_start1.y()) *
robust_cross_product(b1, a1, site1.y() - segm_start2.y(),
site1.x() - segm_start2.x()), 5.0);
- epsilon_robust_comparator< robust_fpt<double> > t;
+ robust_dif< robust_fpt<double> > t;
t -= robust_fpt<double>(a1) * robust_fpt<double>((segm_start1.x() + segm_start2.x()) * 0.5 - site1.x());
t -= robust_fpt<double>(b1) * robust_fpt<double>((segm_start1.y() + segm_start2.y()) * 0.5 - site1.y());
if (point_index == 2) {
@@ -1093,12 +1103,12 @@
t -= det.sqrt();
}
t /= a;
- epsilon_robust_comparator< robust_fpt<double> > c_x, c_y;
+ robust_dif< robust_fpt<double> > c_x, c_y;
c_x += robust_fpt<double>(0.5 * (segm_start1.x() + segm_start2.x()));
c_x += robust_fpt<double>(a1) * t;
c_y += robust_fpt<double>(0.5 * (segm_start1.y() + segm_start2.y()));
c_y += robust_fpt<double>(b1) * t;
- epsilon_robust_comparator< robust_fpt<double> > lower_x(c_x);
+ robust_dif< robust_fpt<double> > lower_x(c_x);
lower_x += robust_fpt<double>(0.5) * c.fabs() / a.sqrt();
recompute_c_x = c_x.dif().ulp() > 128;
recompute_c_y = c_y.dif().ulp() > 128;
@@ -1121,7 +1131,7 @@
robust_fpt<double> c1(robust_cross_product(b1, a1, segm_end1.y(), segm_end1.x()), 2.0);
robust_fpt<double> c2(robust_cross_product(a2, b2, segm_end2.x(), segm_end2.y()), 2.0);
robust_fpt<double> inv_orientation = robust_fpt<double>(1.0) / orientation;
- epsilon_robust_comparator< robust_fpt<double> > t, b, ix, iy;
+ robust_dif< robust_fpt<double> > t, b, ix, iy;
ix += robust_fpt<double>(a2) * c1 * inv_orientation;
ix += robust_fpt<double>(a1) * c2 * inv_orientation;
iy += robust_fpt<double>(b1) * c2 * inv_orientation;
@@ -1140,13 +1150,13 @@
t -= det.sqrt();
}
t /= (a * a);
- epsilon_robust_comparator< robust_fpt<double> > c_x(ix), c_y(iy);
+ robust_dif< robust_fpt<double> > c_x(ix), c_y(iy);
c_x += t * (robust_fpt<double>(a1) * sqr_sum2);
c_x += t * (robust_fpt<double>(a2) * sqr_sum1);
c_y += t * (robust_fpt<double>(b1) * sqr_sum2);
c_y += t * (robust_fpt<double>(b2) * sqr_sum1);
t.abs();
- epsilon_robust_comparator< robust_fpt<double> > lower_x(c_x);
+ robust_dif< robust_fpt<double> > lower_x(c_x);
lower_x += t * orientation.fabs();
recompute_c_x = c_x.dif().ulp() > 128;
recompute_c_y = c_y.dif().ulp() > 128;
@@ -1185,7 +1195,7 @@
robust_fpt<double> cross_12(robust_cross_product(a1.fpv(), b1.fpv(), a2.fpv(), b2.fpv()), 2.0);
robust_fpt<double> cross_23(robust_cross_product(a2.fpv(), b2.fpv(), a3.fpv(), b3.fpv()), 2.0);
robust_fpt<double> cross_31(robust_cross_product(a3.fpv(), b3.fpv(), a1.fpv(), b1.fpv()), 2.0);
- epsilon_robust_comparator< robust_fpt<double> > denom, c_x, c_y, r;
+ robust_dif< robust_fpt<double> > denom, c_x, c_y, r;
// denom = cross_12 * len3 + cross_23 * len1 + cross_31 * len2.
denom += cross_12 * len3;
@@ -1209,7 +1219,7 @@
c_y -= b3 * c2 * len1;
c_y += b3 * c1 * len2;
c_y -= b1 * c3 * len2;
- epsilon_robust_comparator< robust_fpt<double> > lower_x(c_x + r);
+ robust_dif< robust_fpt<double> > lower_x(c_x + r);
bool recompute_c_x = c_x.dif().ulp() > 128;
bool recompute_c_y = c_y.dif().ulp() > 128;
bool recompute_lower_x = lower_x.dif().ulp() > 128;
Modified: sandbox/gtl/boost/polygon/detail/voronoi_fpt_kernel.hpp
==============================================================================
--- sandbox/gtl/boost/polygon/detail/voronoi_fpt_kernel.hpp (original)
+++ sandbox/gtl/boost/polygon/detail/voronoi_fpt_kernel.hpp 2011-10-13 19:51:44 EDT (Thu, 13 Oct 2011)
@@ -14,9 +14,6 @@
namespace polygon {
namespace detail {
- template <typename T>
- class epsilon_robust_comparator;
-
// Represents the result of the epsilon robust predicate.
// If the result is undefined some further processing is usually required.
enum kPredicateResult {
@@ -90,17 +87,22 @@
return value;
}
+ template <>
+ double get_d(const long double& value) {
+ return static_cast<double>(value);
+ }
+
template <typename FPT>
class robust_fpt {
public:
- typedef FPT floating_point_type;
- typedef double relative_error_type;
+ typedef FPT floating_point_type;
+ typedef double relative_error_type;
- // Rounding error is at most 1 EPS.
- static const relative_error_type ROUNDING_ERROR;
+ // Rounding error is at most 1 EPS.
+ static const relative_error_type ROUNDING_ERROR;
- // Constructors.
- robust_fpt() : fpv_(0.0), re_(0) {}
+ // Constructors.
+ robust_fpt() : fpv_(0.0), re_(0) {}
explicit robust_fpt(int fpv) : fpv_(fpv), re_(0) {}
explicit robust_fpt(floating_point_type fpv, bool rounded = true) : fpv_(fpv) {
re_ = rounded ? ROUNDING_ERROR : 0;
@@ -147,44 +149,44 @@
return this->fpv_ > static_cast<floating_point_type>(that);
}
- bool operator==(const robust_fpt &that) const {
- unsigned int ulp = static_cast<unsigned int>(ceil(this->re_ + that.re_));
- return almost_equal(this->fpv_, that.fpv_, ulp);
- }
-
- bool operator!=(const robust_fpt &that) const {
- return !(*this == that);
- }
-
- bool operator<(const robust_fpt &that) const {
- if (*this == that)
- return false;
- return this->fpv_ < that.fpv_;
- }
-
- bool operator>(const robust_fpt &that) const {
- return that < *this;
- }
-
- bool operator<=(const robust_fpt &that) const {
- return !(that < *this);
- }
-
- bool operator>=(const robust_fpt &that) const {
- return !(*this < that);
- }
-
- robust_fpt operator-() const {
- return robust_fpt(-fpv_, re_);
- }
-
- robust_fpt& operator=(const robust_fpt &that) {
- this->fpv_ = that.fpv_;
- this->re_ = that.re_;
- return *this;
- }
+ bool operator==(const robust_fpt &that) const {
+ unsigned int ulp = static_cast<unsigned int>(ceil(this->re_ + that.re_));
+ return almost_equal(this->fpv_, that.fpv_, ulp);
+ }
+
+ bool operator!=(const robust_fpt &that) const {
+ return !(*this == that);
+ }
+
+ bool operator<(const robust_fpt &that) const {
+ if (*this == that)
+ return false;
+ return this->fpv_ < that.fpv_;
+ }
+
+ bool operator>(const robust_fpt &that) const {
+ return that < *this;
+ }
+
+ bool operator<=(const robust_fpt &that) const {
+ return !(that < *this);
+ }
+
+ bool operator>=(const robust_fpt &that) const {
+ return !(*this < that);
+ }
+
+ robust_fpt operator-() const {
+ return robust_fpt(-fpv_, re_);
+ }
+
+ robust_fpt& operator=(const robust_fpt &that) {
+ this->fpv_ = that.fpv_;
+ this->re_ = that.re_;
+ return *this;
+ }
- robust_fpt& operator+=(const robust_fpt &that) {
+ robust_fpt& operator+=(const robust_fpt &that) {
floating_point_type fpv = this->fpv_ + that.fpv_;
if ((this->fpv_ >= 0 && that.fpv_ >= 0) ||
(this->fpv_ <= 0 && that.fpv_ <= 0))
@@ -194,10 +196,10 @@
this->re_ = std::fabs(get_d(temp)) + ROUNDING_ERROR;
}
this->fpv_ = fpv;
- return *this;
- }
+ return *this;
+ }
- robust_fpt& operator-=(const robust_fpt &that) {
+ robust_fpt& operator-=(const robust_fpt &that) {
floating_point_type fpv = this->fpv_ - that.fpv_;
if ((this->fpv_ >= 0 && that.fpv_ <= 0) ||
(this->fpv_ <= 0 && that.fpv_ >= 0))
@@ -207,20 +209,20 @@
this->re_ = std::fabs(get_d(temp)) + ROUNDING_ERROR;
}
this->fpv_ = fpv;
- return *this;
- }
+ return *this;
+ }
- robust_fpt& operator*=(const robust_fpt &that) {
- this->re_ += that.re_ + ROUNDING_ERROR;
- this->fpv_ *= that.fpv_;
+ robust_fpt& operator*=(const robust_fpt &that) {
+ this->re_ += that.re_ + ROUNDING_ERROR;
+ this->fpv_ *= that.fpv_;
return *this;
- }
+ }
- robust_fpt& operator/=(const robust_fpt &that) {
- this->re_ += that.re_ + ROUNDING_ERROR;
- this->fpv_ /= that.fpv_;
+ robust_fpt& operator/=(const robust_fpt &that) {
+ this->re_ += that.re_ + ROUNDING_ERROR;
+ this->fpv_ /= that.fpv_;
return *this;
- }
+ }
robust_fpt operator+(const robust_fpt &that) const {
floating_point_type fpv = this->fpv_ + that.fpv_;
@@ -296,17 +298,17 @@
// Solution: Use ERCs in case of defined comparison results and use
// high-precision libraries for undefined results.
template <typename T>
- class epsilon_robust_comparator {
+ class robust_dif {
public:
- epsilon_robust_comparator() :
+ robust_dif() :
positive_sum_(0),
negative_sum_(0) {}
- epsilon_robust_comparator(const T &value) :
+ robust_dif(const T &value) :
positive_sum_((value>0)?value:0),
negative_sum_((value<0)?-value:0) {}
- epsilon_robust_comparator(const T &pos, const T &neg) :
+ robust_dif(const T &pos, const T &neg) :
positive_sum_(pos),
negative_sum_(neg) {}
@@ -335,7 +337,11 @@
return false;
}
- epsilon_robust_comparator<T> &operator+=(const T &val) {
+ robust_dif<T> operator-() const {
+ return robust_dif(negative_sum_, positive_sum_);
+ }
+
+ robust_dif<T> &operator+=(const T &val) {
if (val >= 0)
positive_sum_ += val;
else
@@ -343,14 +349,14 @@
return *this;
}
- epsilon_robust_comparator<T> &operator+=(
- const epsilon_robust_comparator<T> &erc) {
- positive_sum_ += erc.positive_sum_;
- negative_sum_ += erc.negative_sum_;
+ robust_dif<T> &operator+=(
+ const robust_dif<T> &that) {
+ positive_sum_ += that.positive_sum_;
+ negative_sum_ += that.negative_sum_;
return *this;
}
- epsilon_robust_comparator<T> &operator-=(const T &val) {
+ robust_dif<T> &operator-=(const T &val) {
if (val >= 0)
negative_sum_ += val;
else
@@ -358,14 +364,14 @@
return *this;
}
- epsilon_robust_comparator<T> &operator-=(
- const epsilon_robust_comparator<T> &erc) {
- positive_sum_ += erc.negative_sum_;
- negative_sum_ += erc.positive_sum_;
+ robust_dif<T> &operator-=(
+ const robust_dif<T> &that) {
+ positive_sum_ += that.negative_sum_;
+ negative_sum_ += that.positive_sum_;
return *this;
}
- epsilon_robust_comparator<T> &operator*=(const T &val) {
+ robust_dif<T> &operator*=(const T &val) {
if (val >= 0) {
positive_sum_ *= val;
negative_sum_ *= val;
@@ -377,7 +383,18 @@
return *this;
}
- epsilon_robust_comparator<T> &operator/=(const T &val) {
+ robust_dif<T> &operator*=(
+ const robust_dif<T> &that) {
+ T positive_sum = this->positive_sum_ * that.positive_sum_ +
+ this->negative_sum_ * that.negative_sum_;
+ T negative_sum = this->positive_sum_ * that.negative_sum_ +
+ this->negative_sum_ * that.positive_sum_;
+ positive_sum_ = positive_sum;
+ negative_sum_ = negative_sum;
+ return *this;
+ }
+
+ robust_dif<T> &operator/=(const T &val) {
if (val >= 0) {
positive_sum_ /= val;
negative_sum_ /= val;
@@ -389,129 +406,145 @@
return *this;
}
- // Compare predicate with value using ulp precision.
- kPredicateResult compare(T value, int ulp = 0) const {
- T lhs = positive_sum_ - ((value < 0) ? value : 0);
- T rhs = negative_sum_ + ((value > 0) ? value : 0);
- if (almost_equal(lhs, rhs, ulp))
- return UNDEFINED;
- return (lhs < rhs) ? LESS : MORE;
- }
-
- // Compare two predicats using ulp precision.
- kPredicateResult compare(const epsilon_robust_comparator &rc,
- int ulp = 0) const {
- T lhs = positive_sum_ + rc.neg();
- T rhs = negative_sum_ + rc.pos();
- if (almost_equal(lhs, rhs, ulp))
- return UNDEFINED;
- return (lhs < rhs) ? LESS : MORE;
- }
-
private:
T positive_sum_;
T negative_sum_;
};
template<typename T>
- inline epsilon_robust_comparator<T> operator+(
- const epsilon_robust_comparator<T>& lhs,
- const epsilon_robust_comparator<T>& rhs) {
- return epsilon_robust_comparator<T>(lhs.pos() + rhs.pos(),
- lhs.neg() + rhs.neg());
+ robust_dif<T> operator+(
+ const robust_dif<T>& lhs,
+ const robust_dif<T>& rhs) {
+ return robust_dif<T>(lhs.pos() + rhs.pos(),
+ lhs.neg() + rhs.neg());
+ }
+
+ template<typename T>
+ robust_dif<T> operator+(
+ const robust_dif<T>& lhs, const T& rhs) {
+ if (rhs >= 0) {
+ return robust_dif<T>(lhs.pos() + rhs, lhs.neg());
+ } else {
+ return robust_dif<T>(lhs.pos(), lhs.neg() - rhs);
+ }
+ }
+
+ template<typename T>
+ robust_dif<T> operator+(
+ const T& lhs, const robust_dif<T>& rhs) {
+ if (lhs >= 0) {
+ return robust_dif<T>(lhs + rhs.pos(), rhs.neg());
+ } else {
+ return robust_dif<T>(rhs.pos(), rhs.neg() - lhs);
+ }
}
template<typename T>
- inline epsilon_robust_comparator<T> operator-(
- const epsilon_robust_comparator<T>& lhs,
- const epsilon_robust_comparator<T>& rhs) {
- return epsilon_robust_comparator<T>(lhs.pos() - rhs.neg(),
- lhs.neg() + rhs.pos());
+ robust_dif<T> operator-(
+ const robust_dif<T>& lhs,
+ const robust_dif<T>& rhs) {
+ return robust_dif<T>(lhs.pos() + rhs.neg(), lhs.neg() + rhs.pos());
+ }
+
+ template<typename T>
+ robust_dif<T> operator-(
+ const robust_dif<T>& lhs, const T& rhs) {
+ if (rhs >= 0) {
+ return robust_dif<T>(lhs.pos(), lhs.neg() + rhs);
+ } else {
+ return robust_dif<T>(lhs.pos() - rhs, lhs.neg());
+ }
}
template<typename T>
- inline epsilon_robust_comparator<T> operator*(
- const epsilon_robust_comparator<T>& lhs,
- const epsilon_robust_comparator<T>& rhs) {
- double res_pos = lhs.pos() * rhs.pos() + lhs.neg() * rhs.neg();
- double res_neg = lhs.pos() * rhs.neg() + lhs.neg() * rhs.pos();
- return epsilon_robust_comparator<T>(res_pos, res_neg);
+ robust_dif<T> operator-(
+ const T& lhs, const robust_dif<T>& rhs) {
+ if (lhs >= 0) {
+ return robust_dif<T>(lhs + rhs.neg(), rhs.pos());
+ } else {
+ return robust_dif<T>(rhs.neg(), rhs.pos() - lhs);
+ }
+ }
+
+ template<typename T>
+ robust_dif<T> operator*(
+ const robust_dif<T>& lhs,
+ const robust_dif<T>& rhs) {
+ T res_pos = lhs.pos() * rhs.pos() + lhs.neg() * rhs.neg();
+ T res_neg = lhs.pos() * rhs.neg() + lhs.neg() * rhs.pos();
+ return robust_dif<T>(res_pos, res_neg);
+ }
+
+ template<typename T>
+ robust_dif<T> operator*(
+ const robust_dif<T>& lhs, const T& val) {
+ if (val >= 0) {
+ return robust_dif<T>(lhs.pos() * val, lhs.neg() * val);
+ } else {
+ return robust_dif<T>(-lhs.neg() * val, -lhs.pos() * val);
+ }
}
template<typename T>
- inline epsilon_robust_comparator<T> operator*(
- const epsilon_robust_comparator<T>& lhs, const T& val) {
- if (val >= 0)
- return epsilon_robust_comparator<T>(lhs.pos() * val,
- lhs.neg() * val);
- return epsilon_robust_comparator<T>(-lhs.neg() * val,
- -lhs.pos() * val);
+ robust_dif<T> operator*(
+ const T& val, const robust_dif<T>& rhs) {
+ if (val >= 0) {
+ return robust_dif<T>(val * rhs.pos(), val * rhs.neg());
+ } else {
+ return robust_dif<T>(-val * rhs.neg(), -val * rhs.pos());
+ }
}
template<typename T>
- inline epsilon_robust_comparator<T> operator*(
- const T& val, const epsilon_robust_comparator<T>& rhs) {
- if (val >= 0)
- return epsilon_robust_comparator<T>(val * rhs.pos(),
- val * rhs.neg());
- return epsilon_robust_comparator<T>(-val * rhs.neg(),
- -val * rhs.pos());
+ robust_dif<T> operator/(
+ const robust_dif<T>& lhs, const T& val) {
+ if (val >= 0) {
+ return robust_dif<T>(lhs.pos() / val, lhs.neg() / val);
+ } else {
+ return robust_dif<T>(-lhs.neg() / val, -lhs.pos() / val);
+ }
}
///////////////////////////////////////////////////////////////////////////
- // VORONOI SQRT EXPR EVALUATOR ////////////////////////////////////////////
+ // VORONOI ROBUST SQRT EXPRESSION ////////////////////////////////////////
///////////////////////////////////////////////////////////////////////////
- template <int N>
- struct sqrt_expr_evaluator {
- template <typename mpt, typename mpf>
- static mpf eval(mpt *A, mpt *B);
- };
-
- // Evaluates expression:
- // A[0] * sqrt(B[0]).
- template <>
- struct sqrt_expr_evaluator<1> {
- template <typename mpt, typename mpf>
- static mpf eval(mpt *A, mpt *B) {
+ template <typename mpt, typename mpf>
+ class robust_sqrt_expr {
+ public:
+ // Evaluates expression:
+ // A[0] * sqrt(B[0]).
+ mpf eval1(mpt *A, mpt *B) {
#ifndef THREAD_SAFETY
static
#endif
- mpf a, b, ret_val;
- a = A[0];
- b = B[0];
- ret_val = a * sqrt(b);
- return ret_val;
+ mpf lhs, rhs, numer;
+ lhs = A[0];
+ rhs = B[0];
+ numer = lhs * sqrt(rhs);
+ return numer;
}
- };
- // Evaluates expression:
- // A[0] * sqrt(B[0]) + A[1] * sqrt(B[1]) with
- // 7 * EPS relative error in the worst case.
- template <>
- struct sqrt_expr_evaluator<2> {
- template <typename mpt, typename mpf>
- static mpf eval(mpt *A, mpt *B) {
+ // Evaluates expression:
+ // A[0] * sqrt(B[0]) + A[1] * sqrt(B[1]) with
+ // 7 * EPS relative error in the worst case.
+ mpf eval2(mpt *A, mpt *B) {
#ifndef THREAD_SAFETY
static
#endif
- mpf lhs, rhs, numerator;
- lhs = sqrt_expr_evaluator<1>::eval<mpt, mpf>(A, B);
- rhs = sqrt_expr_evaluator<1>::eval<mpt, mpf>(A + 1, B + 1);
+ mpf lhs, rhs, numer;
+ lhs = eval1(A, B);
+ rhs = eval1(A + 1, B + 1);
if ((lhs >= 0 && rhs >= 0) || (lhs <= 0 && rhs <= 0))
return lhs + rhs;
- numerator = A[0] * A[0] * B[0] - A[1] * A[1] * B[1];
- return numerator / (lhs - rhs);
+ numer = A[0] * A[0] * B[0] - A[1] * A[1] * B[1];
+ return numer / (lhs - rhs);
}
- };
- // Evaluates expression:
- // A[0] * sqrt(B[0]) + A[1] * sqrt(B[1]) + A[2] * sqrt(B[2])
- // with 16 * EPS relative error in the worst case.
- template <>
- struct sqrt_expr_evaluator<3> {
- template <typename mpt, typename mpf>
- static mpf eval(mpt *A, mpt *B) {
+ // Evaluates expression:
+ // A[0] * sqrt(B[0]) + A[1] * sqrt(B[1]) + A[2] * sqrt(B[2])
+ // with 16 * EPS relative error in the worst case.
+ mpf eval3(mpt *A, mpt *B) {
#ifndef THREAD_SAFETY
static
#endif
@@ -520,8 +553,8 @@
static
#endif
mpf lhs, rhs, numer;
- lhs = sqrt_expr_evaluator<2>::eval<mpt, mpf>(A, B);
- rhs = sqrt_expr_evaluator<1>::eval<mpt, mpf>(A + 2, B + 2);
+ lhs = eval2(A, B);
+ rhs = eval1(A + 2, B + 2);
if ((lhs >= 0 && rhs >= 0) || (lhs <= 0 && rhs <= 0))
return lhs + rhs;
cA[0] = A[0] * A[0] * B[0] + A[1] * A[1] * B[1];
@@ -529,18 +562,15 @@
cB[0] = 1;
cA[1] = A[0] * A[1] * 2;
cB[1] = B[0] * B[1];
- numer = sqrt_expr_evaluator<2>::eval<mpt, mpf>(cA, cB);
+ numer = eval2(cA, cB);
return numer / (lhs - rhs);
}
- };
- // Evaluates expression:
- // A[0] * sqrt(B[0]) + A[1] * sqrt(B[1]) + A[2] * sqrt(B[2]) + A[3] * sqrt(B[3])
- // with 25 * EPS relative error in the worst case.
- template <>
- struct sqrt_expr_evaluator<4> {
- template <typename mpt, typename mpf>
- static mpf eval(mpt *A, mpt *B) {
+
+ // Evaluates expression:
+ // A[0] * sqrt(B[0]) + A[1] * sqrt(B[1]) + A[2] * sqrt(B[2]) + A[3] * sqrt(B[3])
+ // with 25 * EPS relative error in the worst case.
+ mpf eval4(mpt *A, mpt *B) {
#ifndef THREAD_SAFETY
static
#endif
@@ -549,8 +579,8 @@
static
#endif
mpf lhs, rhs, numer;
- lhs = sqrt_expr_evaluator<2>::eval<mpt, mpf>(A, B);
- rhs = sqrt_expr_evaluator<2>::eval<mpt, mpf>(A + 2, B + 2);
+ lhs = eval2(A, B);
+ rhs = eval2(A + 2, B + 2);
if ((lhs >= 0 && rhs >= 0) || (lhs <= 0 && rhs <= 0))
return lhs + rhs;
cA[0] = A[0] * A[0] * B[0] + A[1] * A[1] * B[1];
@@ -560,9 +590,11 @@
cB[1] = B[0] * B[1];
cA[2] = A[2] * A[3] * -2;
cB[2] = B[2] * B[3];
- numer = sqrt_expr_evaluator<3>::eval<mpt, mpf>(cA, cB);
+ numer = eval3(cA, cB);
return numer / (lhs - rhs);
}
+
+ private:
};
} // detail
Modified: sandbox/gtl/boost/polygon/voronoi_builder.hpp
==============================================================================
--- sandbox/gtl/boost/polygon/voronoi_builder.hpp (original)
+++ sandbox/gtl/boost/polygon/voronoi_builder.hpp 2011-10-13 19:51:44 EDT (Thu, 13 Oct 2011)
@@ -18,12 +18,8 @@
#include <queue>
#include <vector>
-#pragma warning(disable:4800)
-#include <gmpxx.h>
-
#include "polygon.hpp"
-#include "detail/mpt_wrapper.hpp"
#include "detail/voronoi_calc_kernel.hpp"
#include "detail/voronoi_fpt_kernel.hpp"
#include "detail/voronoi_internal_structures.hpp"
Modified: sandbox/gtl/libs/polygon/test/Jamfile.v2
==============================================================================
--- sandbox/gtl/libs/polygon/test/Jamfile.v2 (original)
+++ sandbox/gtl/libs/polygon/test/Jamfile.v2 2011-10-13 19:51:44 EDT (Thu, 13 Oct 2011)
@@ -30,6 +30,7 @@
:
[ run voronoi_builder_test.cpp ]
[ run voronoi_clipping_test.cpp ]
+ [ run voronoi_fpt_kernel_test.cpp ]
[ run voronoi_internal_structures_test.cpp ]
;
Modified: sandbox/gtl/libs/polygon/test/voronoi_clipping_test.cpp
==============================================================================
--- sandbox/gtl/libs/polygon/test/voronoi_clipping_test.cpp (original)
+++ sandbox/gtl/libs/polygon/test/voronoi_clipping_test.cpp 2011-10-13 19:51:44 EDT (Thu, 13 Oct 2011)
@@ -17,6 +17,8 @@
typedef VH::brect_type rect_type;
typedef VH::point_type point_type;
+#define SZ(st) static_cast<int>(st.size())
+
// Test segment clipping.
BOOST_AUTO_TEST_CASE(segment_clipping_test1) {
rect_type rect(0, -1, 4, 2);
@@ -31,13 +33,13 @@
std::vector<point_type> intersections;
VH::find_intersections(origin, direction1_1, VH::SEGMENT, rect, intersections);
- BOOST_CHECK_EQUAL(intersections.size(), 2);
+ BOOST_CHECK_EQUAL(SZ(intersections), 2);
BOOST_CHECK_EQUAL(intersections[0].x() == 0 && intersections[0].y() == 2, true);
BOOST_CHECK_EQUAL(intersections[1].x() == 3 && intersections[1].y() == -1, true);
intersections.clear();
VH::find_intersections(origin, direction1_2, VH::SEGMENT, rect, intersections);
- BOOST_CHECK_EQUAL(intersections.size(), 1);
+ BOOST_CHECK_EQUAL(SZ(intersections), 1);
BOOST_CHECK_EQUAL(intersections[0].x() == 0 && intersections[0].y() == 2, true);
intersections.clear();
@@ -46,23 +48,23 @@
intersections.clear();
VH::find_intersections(origin, direction2, VH::SEGMENT, rect, intersections);
- BOOST_CHECK_EQUAL(intersections.size(), 2);
+ BOOST_CHECK_EQUAL(SZ(intersections), 2);
BOOST_CHECK_EQUAL(intersections[0].x() == 0 && intersections[0].y() == 1, true);
BOOST_CHECK_EQUAL(intersections[1].x() == 1 && intersections[1].y() == -1, true);
intersections.clear();
VH::find_intersections(origin, direction3, VH::SEGMENT, rect, intersections);
- BOOST_CHECK_EQUAL(intersections.size(), 1);
+ BOOST_CHECK_EQUAL(SZ(intersections), 1);
BOOST_CHECK_EQUAL(intersections[0].x() == 1 && intersections[0].y() == 2, true);
intersections.clear();
VH::find_intersections(origin, direction4, VH::SEGMENT, rect, intersections);
- BOOST_CHECK_EQUAL(intersections.size(), 1);
+ BOOST_CHECK_EQUAL(SZ(intersections), 1);
BOOST_CHECK_EQUAL(intersections[0].x() == 0 && intersections[0].y() == -1, true);
intersections.clear();
VH::find_intersections(origin, direction5, VH::SEGMENT, rect, intersections);
- BOOST_CHECK_EQUAL(intersections.size(), 1);
+ BOOST_CHECK_EQUAL(SZ(intersections), 1);
BOOST_CHECK_EQUAL(intersections[0].x() == 4 && intersections[0].y() == 2, true);
intersections.clear();
}
@@ -78,9 +80,9 @@
point_type direction(i, j);
VH::find_intersections(origin, direction, VH::SEGMENT, rect, intersections);
if (abs(i) >= 2 || abs(j) >= 2)
- BOOST_CHECK_EQUAL(intersections.size(), 1);
+ BOOST_CHECK_EQUAL(SZ(intersections), 1);
else
- BOOST_CHECK_EQUAL(intersections.size(), 0);
+ BOOST_CHECK_EQUAL(SZ(intersections), 0);
}
}
@@ -96,7 +98,7 @@
double y = 1.0 * j / 26;
point_type direction(x, y);
VH::find_intersections(origin, direction, VH::SEGMENT, rect, intersections);
- BOOST_CHECK_EQUAL(intersections.size(), 0);
+ BOOST_CHECK_EQUAL(SZ(intersections), 0);
}
}
@@ -112,30 +114,30 @@
std::vector<point_type> intersections;
VH::find_intersections(origin, direction1, VH::RAY, rect, intersections);
- BOOST_CHECK_EQUAL(intersections.size(), 2);
+ BOOST_CHECK_EQUAL(SZ(intersections), 2);
BOOST_CHECK_EQUAL(intersections[0].x() == 0 && intersections[0].y() == 2, true);
BOOST_CHECK_EQUAL(intersections[1].x() == 3 && intersections[1].y() == -1, true);
intersections.clear();
VH::find_intersections(origin, direction2, VH::RAY, rect, intersections);
- BOOST_CHECK_EQUAL(intersections.size(), 2);
+ BOOST_CHECK_EQUAL(SZ(intersections), 2);
BOOST_CHECK_EQUAL(intersections[0].x() == 0 && intersections[0].y() == 1, true);
BOOST_CHECK_EQUAL(intersections[1].x() == 1 && intersections[1].y() == -1, true);
intersections.clear();
VH::find_intersections(origin, direction3, VH::RAY, rect, intersections);
- BOOST_CHECK_EQUAL(intersections.size(), 2);
+ BOOST_CHECK_EQUAL(SZ(intersections), 2);
BOOST_CHECK_EQUAL(intersections[0].x() == 1 && intersections[0].y() == 2, true);
BOOST_CHECK_EQUAL(intersections[1].x() == 4 && intersections[1].y() == 0.5, true);
intersections.clear();
VH::find_intersections(origin, direction4, VH::RAY, rect, intersections);
- BOOST_CHECK_EQUAL(intersections.size(), 1);
+ BOOST_CHECK_EQUAL(SZ(intersections), 1);
BOOST_CHECK_EQUAL(intersections[0].x() == 0 && intersections[0].y() == -1, true);
intersections.clear();
VH::find_intersections(origin, direction5, VH::RAY, rect, intersections);
- BOOST_CHECK_EQUAL(intersections.size(), 1);
+ BOOST_CHECK_EQUAL(SZ(intersections), 1);
BOOST_CHECK_EQUAL(intersections[0].x() == 4 && intersections[0].y() == 2, true);
intersections.clear();
}
@@ -153,7 +155,7 @@
double y = 1.0 * j / 26;
point_type direction(x, y);
VH::find_intersections(origin, direction, VH::RAY, rect, intersections);
- BOOST_CHECK_EQUAL(intersections.size(), 1);
+ BOOST_CHECK_EQUAL(SZ(intersections), 1);
}
}
@@ -169,30 +171,30 @@
std::vector<point_type> intersections;
VH::find_intersections(origin, direction1, VH::LINE, rect, intersections);
- BOOST_CHECK_EQUAL(intersections.size(), 2);
+ BOOST_CHECK_EQUAL(SZ(intersections), 2);
BOOST_CHECK_EQUAL(intersections[0].x() == 3 && intersections[0].y() == -1, true);
BOOST_CHECK_EQUAL(intersections[1].x() == 0 && intersections[1].y() == 2, true);
intersections.clear();
VH::find_intersections(origin, direction2, VH::LINE, rect, intersections);
- BOOST_CHECK_EQUAL(intersections.size(), 2);
+ BOOST_CHECK_EQUAL(SZ(intersections), 2);
BOOST_CHECK_EQUAL(intersections[0].x() == 1 && intersections[0].y() == -1, true);
BOOST_CHECK_EQUAL(intersections[1].x() == 0 && intersections[1].y() == 1, true);
intersections.clear();
VH::find_intersections(origin, direction3, VH::LINE, rect, intersections);
- BOOST_CHECK_EQUAL(intersections.size(), 2);
+ BOOST_CHECK_EQUAL(SZ(intersections), 2);
BOOST_CHECK_EQUAL(intersections[0].x() == 4 && intersections[0].y() == 0.5, true);
BOOST_CHECK_EQUAL(intersections[1].x() == 1 && intersections[1].y() == 2, true);
intersections.clear();
VH::find_intersections(origin, direction4, VH::LINE, rect, intersections);
- BOOST_CHECK_EQUAL(intersections.size(), 1);
+ BOOST_CHECK_EQUAL(SZ(intersections), 1);
BOOST_CHECK_EQUAL(intersections[0].x() == 0 && intersections[0].y() == -1, true);
intersections.clear();
VH::find_intersections(origin, direction5, VH::LINE, rect, intersections);
- BOOST_CHECK_EQUAL(intersections.size(), 1);
+ BOOST_CHECK_EQUAL(SZ(intersections), 1);
BOOST_CHECK_EQUAL(intersections[0].x() == 4 && intersections[0].y() == 2, true);
intersections.clear();
}
@@ -210,6 +212,6 @@
double y = 1.0 * j / 26;
point_type direction(x, y);
VH::find_intersections(origin, direction, VH::LINE, rect, intersections);
- BOOST_CHECK_EQUAL(intersections.size(), 2);
+ BOOST_CHECK_EQUAL(SZ(intersections), 2);
}
}
Added: sandbox/gtl/libs/polygon/test/voronoi_fpt_kernel_test.cpp
==============================================================================
--- (empty file)
+++ sandbox/gtl/libs/polygon/test/voronoi_fpt_kernel_test.cpp 2011-10-13 19:51:44 EDT (Thu, 13 Oct 2011)
@@ -0,0 +1,286 @@
+// Boost.Polygon library voronoi_fpt_kernel_test.cpp file
+
+// Copyright Andrii Sydorchuk 2010-2011.
+// Distributed under the Boost Software License, Version 1.0.
+// (See accompanying file LICENSE_1_0.txt or copy at
+// http://www.boost.org/LICENSE_1_0.txt)
+
+// See http://www.boost.org for updates, documentation, and revision history.
+
+#include <cmath>
+#include <ctime>
+#include <vector>
+
+#define BOOST_TEST_MODULE voronoi_fpt_kernel_test
+#include <boost/mpl/list.hpp>
+#include <boost/random/mersenne_twister.hpp>
+#include <boost/test/test_case_template.hpp>
+
+#include "boost/polygon/detail/voronoi_fpt_kernel.hpp"
+#include "boost/polygon/detail/voronoi_calc_kernel.hpp"
+using namespace boost::polygon::detail;
+
+typedef boost::mpl::list<double, mpf_class, mpt_wrapper<mpf_class, 8> > test_types;
+
+BOOST_AUTO_TEST_CASE_TEMPLATE(robust_fpt_constructors_test1, T, test_types) {
+ typedef robust_fpt<T> rfpt_type;
+ rfpt_type a = rfpt_type();
+ BOOST_CHECK_EQUAL(a.fpv() == static_cast<T>(0), true);
+ BOOST_CHECK_EQUAL(a.re() == 0.0, true);
+ BOOST_CHECK_EQUAL(a.ulp() == 0, true);
+}
+
+BOOST_AUTO_TEST_CASE_TEMPLATE(robust_fpt_constructors_test2, T, test_types) {
+ typedef robust_fpt<T> rfpt_type;
+ rfpt_type a(static_cast<T>(10));
+ BOOST_CHECK_EQUAL(a.fpv() == static_cast<T>(10), true);
+ BOOST_CHECK_EQUAL(a.re() == 1.0, true);
+ BOOST_CHECK_EQUAL(a.ulp() == 1.0, true);
+}
+
+BOOST_AUTO_TEST_CASE_TEMPLATE(robust_fpt_constructors_test3, T, test_types) {
+ typedef robust_fpt<T> rfpt_type;
+ rfpt_type a(static_cast<T>(10), false);
+ BOOST_CHECK_EQUAL(a.fpv() == static_cast<T>(10), true);
+ BOOST_CHECK_EQUAL(a.re() == 0.0, true);
+ BOOST_CHECK_EQUAL(a.ulp() == 0.0, true);
+}
+
+BOOST_AUTO_TEST_CASE_TEMPLATE(robust_fpt_constructors_test4, T, test_types) {
+ typedef robust_fpt<T> rfpt_type;
+ rfpt_type a(static_cast<T>(10), 3.0);
+ BOOST_CHECK_EQUAL(a.fpv() == static_cast<T>(10), true);
+ BOOST_CHECK_EQUAL(a.re() == 3.0, true);
+ BOOST_CHECK_EQUAL(a.ulp() == 3.0, true);
+
+ rfpt_type b(static_cast<T>(10), 2.75);
+ BOOST_CHECK_EQUAL(b.fpv() == static_cast<T>(10), true);
+ BOOST_CHECK_EQUAL(b.re() == 2.75, true);
+ BOOST_CHECK_EQUAL(b.ulp() == 2.75, true);
+}
+
+BOOST_AUTO_TEST_CASE_TEMPLATE(robust_fpt_sum_test1, T, test_types) {
+ typedef robust_fpt<T> rfpt_type;
+ rfpt_type a(static_cast<T>(2), 5.0);
+ rfpt_type b(static_cast<T>(3), 4.0);
+ rfpt_type c = a + b;
+ BOOST_CHECK_EQUAL(c.fpv() == static_cast<T>(5), true);
+ BOOST_CHECK_EQUAL(c.re() == 6.0, true);
+ BOOST_CHECK_EQUAL(c.ulp() == 6.0, true);
+
+ c += b;
+ BOOST_CHECK_EQUAL(c.fpv() == static_cast<T>(8), true);
+ BOOST_CHECK_EQUAL(c.re() == 7.0, true);
+ BOOST_CHECK_EQUAL(c.ulp() == 7.0, true);
+}
+
+BOOST_AUTO_TEST_CASE_TEMPLATE(robust_fpt_sum_test2, T, test_types) {
+ typedef robust_fpt<T> rfpt_type;
+ rfpt_type a(static_cast<T>(3), 2.0);
+ rfpt_type b(static_cast<T>(-2), 3.0);
+ rfpt_type c = a + b;
+ BOOST_CHECK_EQUAL(c.fpv() == static_cast<T>(1), true);
+ BOOST_CHECK_EQUAL(c.re() == 13.0, true);
+ BOOST_CHECK_EQUAL(c.ulp() == 13.0, true);
+
+ c += b;
+ BOOST_CHECK_EQUAL(c.fpv() == static_cast<T>(-1), true);
+ BOOST_CHECK_EQUAL(c.re() == 20.0, true);
+ BOOST_CHECK_EQUAL(c.ulp() == 20.0, true);
+}
+
+BOOST_AUTO_TEST_CASE_TEMPLATE(robust_fpt_dif_test1, T, test_types) {
+ typedef robust_fpt<T> rfpt_type;
+ rfpt_type a(static_cast<T>(2), 5.0);
+ rfpt_type b(static_cast<T>(-3), 4.0);
+ rfpt_type c = a - b;
+ BOOST_CHECK_EQUAL(c.fpv() == static_cast<T>(5), true);
+ BOOST_CHECK_EQUAL(c.re() == 6.0, true);
+ BOOST_CHECK_EQUAL(c.ulp() == 6.0, true);
+
+ c -= b;
+ BOOST_CHECK_EQUAL(c.fpv() == static_cast<T>(8), true);
+ BOOST_CHECK_EQUAL(c.re() == 7.0, true);
+ BOOST_CHECK_EQUAL(c.ulp() == 7.0, true);
+}
+
+BOOST_AUTO_TEST_CASE_TEMPLATE(robust_fpt_dif_test2, T, test_types) {
+ typedef robust_fpt<T> rfpt_type;
+ rfpt_type a(static_cast<T>(3), 2.0);
+ rfpt_type b(static_cast<T>(2), 3.0);
+ rfpt_type c = a - b;
+ BOOST_CHECK_EQUAL(c.fpv() == static_cast<T>(1), true);
+ BOOST_CHECK_EQUAL(c.re() == 13.0, true);
+ BOOST_CHECK_EQUAL(c.ulp() == 13.0, true);
+
+ c -= b;
+ BOOST_CHECK_EQUAL(c.fpv() == static_cast<T>(-1), true);
+ BOOST_CHECK_EQUAL(c.re() == 20.0, true);
+ BOOST_CHECK_EQUAL(c.ulp() == 20.0, true);
+}
+
+BOOST_AUTO_TEST_CASE_TEMPLATE(robust_fpt_mult_test3, T, test_types) {
+ typedef robust_fpt<T> rfpt_type;
+ rfpt_type a(static_cast<T>(2), 3.0);
+ rfpt_type b(static_cast<T>(4));
+ rfpt_type c = a * b;
+ BOOST_CHECK_EQUAL(c.fpv() == static_cast<T>(8), true);
+ BOOST_CHECK_EQUAL(c.re() == 5.0, true);
+ BOOST_CHECK_EQUAL(c.ulp() == 5.0, true);
+
+ c *= b;
+ BOOST_CHECK_EQUAL(c.fpv() == static_cast<T>(32), true);
+ BOOST_CHECK_EQUAL(c.re() == 7.0, true);
+ BOOST_CHECK_EQUAL(c.ulp() == 7.0, true);
+}
+
+BOOST_AUTO_TEST_CASE_TEMPLATE(robust_fpt_div_test1, T, test_types) {
+ typedef robust_fpt<T> rfpt_type;
+ rfpt_type a(static_cast<T>(2), 3.0);
+ rfpt_type b(static_cast<T>(4));
+ rfpt_type c = a / b;
+ BOOST_CHECK_EQUAL(c.fpv() == static_cast<T>(0.5), true);
+ BOOST_CHECK_EQUAL(c.re() == 5.0, true);
+ BOOST_CHECK_EQUAL(c.ulp() == 5.0, true);
+
+ c /= b;
+ BOOST_CHECK_EQUAL(c.fpv() == static_cast<T>(0.125), true);
+ BOOST_CHECK_EQUAL(c.re() == 7.0, true);
+ BOOST_CHECK_EQUAL(c.ulp() == 7.0, true);
+}
+
+BOOST_AUTO_TEST_CASE(robust_dif_constructors_test) {
+ robust_dif<int> rd1;
+ BOOST_CHECK_EQUAL(rd1.pos(), 0);
+ BOOST_CHECK_EQUAL(rd1.neg(), 0);
+ BOOST_CHECK_EQUAL(rd1.dif(), 0);
+
+ robust_dif<int> rd2(1);
+ BOOST_CHECK_EQUAL(rd2.pos(), 1);
+ BOOST_CHECK_EQUAL(rd2.neg(), 0);
+ BOOST_CHECK_EQUAL(rd2.dif(), 1);
+
+ robust_dif<int> rd3(-1);
+ BOOST_CHECK_EQUAL(rd3.pos(), 0);
+ BOOST_CHECK_EQUAL(rd3.neg(), 1);
+ BOOST_CHECK_EQUAL(rd3.dif(), -1);
+
+ robust_dif<int> rd4(1, 2);
+ BOOST_CHECK_EQUAL(rd4.pos(), 1);
+ BOOST_CHECK_EQUAL(rd4.neg(), 2);
+ BOOST_CHECK_EQUAL(rd4.dif(), -1);
+}
+
+BOOST_AUTO_TEST_CASE(robust_dif_operators_test1) {
+ robust_dif<int> a(5, 2), b(1, 10);
+ int dif_a = a.dif();
+ int dif_b = b.dif();
+ robust_dif<int> sum = a + b;
+ robust_dif<int> dif = a - b;
+ robust_dif<int> mult = a * b;
+ robust_dif<int> umin = -a;
+ BOOST_CHECK_EQUAL(sum.dif(), dif_a + dif_b);
+ BOOST_CHECK_EQUAL(dif.dif(), dif_a - dif_b);
+ BOOST_CHECK_EQUAL(mult.dif(), dif_a * dif_b);
+ BOOST_CHECK_EQUAL(umin.dif(), -dif_a);
+}
+
+BOOST_AUTO_TEST_CASE(robust_dif_operators_test2) {
+ robust_dif<int> a(5, 2);
+ for (int b = -3; b <= 3; b += 6) {
+ int dif_a = a.dif();
+ int dif_b = b;
+ robust_dif<int> sum = a + b;
+ robust_dif<int> dif = a - b;
+ robust_dif<int> mult = a * b;
+ robust_dif<int> div = a / b;
+ BOOST_CHECK_EQUAL(sum.dif(), dif_a + dif_b);
+ BOOST_CHECK_EQUAL(dif.dif(), dif_a - dif_b);
+ BOOST_CHECK_EQUAL(mult.dif(), dif_a * dif_b);
+ BOOST_CHECK_EQUAL(div.dif(), dif_a / dif_b);
+ }
+}
+
+BOOST_AUTO_TEST_CASE(robust_dif_operators_test3) {
+ robust_dif<int> b(5, 2);
+ for (int a = -3; a <= 3; a += 6) {
+ int dif_a = a;
+ int dif_b = b.dif();
+ robust_dif<int> sum = a + b;
+ robust_dif<int> dif = a - b;
+ robust_dif<int> mult = a * b;
+ BOOST_CHECK_EQUAL(sum.dif(), dif_a + dif_b);
+ BOOST_CHECK_EQUAL(dif.dif(), dif_a - dif_b);
+ BOOST_CHECK_EQUAL(mult.dif(), dif_a * dif_b);
+ }
+}
+
+BOOST_AUTO_TEST_CASE(robust_dif_operators_test4) {
+ std::vector< robust_dif<int> > a4(4, robust_dif<int>(5, 2));
+ std::vector< robust_dif<int> > b4(4, robust_dif<int>(1, 2));
+ std::vector< robust_dif<int> > c4 = a4;
+ c4[0] += b4[0];
+ c4[1] -= b4[1];
+ c4[2] *= b4[2];
+ BOOST_CHECK_EQUAL(c4[0].dif(), a4[0].dif() + b4[0].dif());
+ BOOST_CHECK_EQUAL(c4[1].dif(), a4[1].dif() - b4[1].dif());
+ BOOST_CHECK_EQUAL(c4[2].dif(), a4[2].dif() * b4[2].dif());
+ a4[0] += b4[0].dif();
+ a4[1] -= b4[1].dif();
+ a4[2] *= b4[2].dif();
+ a4[3] /= b4[3].dif();
+ BOOST_CHECK_EQUAL(c4[0].dif(), a4[0].dif());
+ BOOST_CHECK_EQUAL(c4[1].dif(), a4[1].dif());
+ BOOST_CHECK_EQUAL(c4[2].dif(), a4[2].dif());
+ BOOST_CHECK_EQUAL(c4[3].dif() / b4[3].dif(), a4[3].dif());
+}
+
+template <typename mpz, typename mpf>
+class sqrt_expr_tester {
+public:
+ bool run() {
+ static boost::mt19937 gen(static_cast<unsigned int>(time(NULL)));
+ bool ret_val = true;
+ for (int i = 0; i < MX_SQRTS; ++i) {
+ a[i] = gen() % 1000000;
+ double random_int = gen() % 1000000;
+ b[i] = random_int * random_int;
+ }
+ int mask = (1 << MX_SQRTS);
+ for (int i = 0; i < mask; i++) {
+ double expected_val = 0.0;
+ for (int j = 0; j < MX_SQRTS; j++) {
+ if (i & (1 << j)) {
+ A[j] = a[j];
+ B[j] = b[j];
+ expected_val += a[j] * sqrt(b[j]);
+ } else {
+ A[j] = -a[j];
+ B[j] = b[j];
+ expected_val -= a[j] * sqrt(b[j]);
+ }
+ }
+ mpf received_val = (sqrt_expr_.eval4(A, B));
+ ret_val &= almost_equal(expected_val, received_val.get_d(), 1);
+ }
+ return ret_val;
+ }
+private:
+ static const int MX_SQRTS = 4;
+ robust_sqrt_expr<mpz, mpf> sqrt_expr_;
+ mpz A[MX_SQRTS];
+ mpz B[MX_SQRTS];
+ double a[MX_SQRTS];
+ double b[MX_SQRTS];
+};
+
+BOOST_AUTO_TEST_CASE(mpz_sqrt_evaluator_test) {
+ typedef mpt_wrapper<mpz_class, 8> mpz_wrapper_type;
+ typedef mpt_wrapper<mpf_class, 2> mpf_wrapper_type;
+ sqrt_expr_tester<mpz_class, mpf_class> tester1;
+ sqrt_expr_tester<mpz_wrapper_type, mpf_wrapper_type> tester2;
+ for (int i = 0; i < 2000; ++i) {
+ BOOST_CHECK(tester1.run());
+ BOOST_CHECK(tester2.run());
+ }
+}