|
|
Boost-Commit : |
Subject: [Boost-commit] svn:boost r76135 - in sandbox/gtl: boost/polygon boost/polygon/detail libs/polygon libs/polygon/test
From: sydorchuk.andriy_at_[hidden]
Date: 2011-12-24 12:20:32
Author: asydorchuk
Date: 2011-12-24 12:20:30 EST (Sat, 24 Dec 2011)
New Revision: 76135
URL: http://svn.boost.org/trac/boost/changeset/76135
Log:
Adding extended exponent floating point type implementation.
Adding extended int (big integer) implementation.
Adding new tests for extended exponent fpt and extended int.
Removing dependencies on gmp (not using anymore).
Removing mpt_wrapper.
Removed:
sandbox/gtl/boost/polygon/detail/mpt_wrapper.hpp
Text files modified:
sandbox/gtl/boost/polygon/detail/voronoi_calc_utils.hpp | 458 ++++++++++----------
sandbox/gtl/boost/polygon/detail/voronoi_robust_fpt.hpp | 888 +++++++++++++++++++++++++++++++++++----
sandbox/gtl/boost/polygon/directed_line_segment_set_data.hpp | 2
sandbox/gtl/boost/polygon/voronoi_diagram.hpp | 2
sandbox/gtl/libs/polygon/Jamroot | 31 -
sandbox/gtl/libs/polygon/test/voronoi_calc_utils_test.cpp | 20
sandbox/gtl/libs/polygon/test/voronoi_robust_fpt_test.cpp | 462 +++++++++++++++++---
sandbox/gtl/libs/polygon/test/voronoi_test_helper.hpp | 8
8 files changed, 1424 insertions(+), 447 deletions(-)
Deleted: sandbox/gtl/boost/polygon/detail/mpt_wrapper.hpp
==============================================================================
--- sandbox/gtl/boost/polygon/detail/mpt_wrapper.hpp 2011-12-24 12:20:30 EST (Sat, 24 Dec 2011)
+++ (empty file)
@@ -1,221 +0,0 @@
-// Boost polygon/detail/mpt_wrapper.hpp header file
-
-// Copyright Andrii Sydorchuk 2010.
-// 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.
-
-#ifndef BOOST_POLYGON_MPT_WRAPPER
-#define BOOST_POLYGON_MPT_WRAPPER
-
-#include <cmath>
-#include <string>
-
-namespace boost {
-namespace polygon {
-namespace detail {
-
- // Allows to compute expression without allocation of additional memory.
- // Expression evaluation should contain less then N temporary variables
- // (e.g. (a1 * b1) + (a2 * b2) - requires two temporary variables to hold
- // a1 * b1 and a2 * b2). This class is not thread safe, use mpz_class or
- // mpq_class instead (however there'll be at least 2 times slowdown).
- template <typename mpt, int N>
- class mpt_wrapper {
- public:
- mpt_wrapper() {}
-
- explicit mpt_wrapper(int input) : m_(input) {}
-
- explicit mpt_wrapper(double input) : m_(input) {}
-
- mpt_wrapper(const mpt& input) : m_(input) {}
-
- mpt_wrapper(const mpt_wrapper& w) : m_(w.m_) {}
-
- template <typename mpt2, int N2>
- mpt_wrapper(const mpt_wrapper<mpt2, N2>& w) : m_(w.m_) {}
-
- mpt_wrapper& operator=(int that) {
- m_ = that;
- return *this;
- }
-
- mpt_wrapper& operator=(double that) {
- m_ = that;
- return *this;
- }
-
- mpt_wrapper& operator=(const mpt_wrapper& that) {
- m_ = that.m_;
- return *this;
- }
-
- template <typename mpt2, int N2>
- mpt_wrapper& operator=(const mpt_wrapper<mpt2, N2>& that) {
- m_ = that.get_mpt();
- return *this;
- }
-
- double get_d() const {
- return m_.get_d();
- }
-
- std::string get_str() const {
- return m_.get_str();
- }
-
- const mpt& get_mpt() const {
- return m_;
- }
-
- bool operator==(const mpt_wrapper& that) const {
- return m_ == that.m_;
- }
-
- bool operator!=(const mpt_wrapper& that) const {
- return m_ != that.m_;
- }
-
- bool operator<(const mpt_wrapper& that) const {
- return m_ < that.m_;
- }
-
- bool operator<=(const mpt_wrapper& that) const {
- return m_ <= that.m_;
- }
-
- bool operator>(const mpt_wrapper& that) const {
- return m_ > that.m_;
- }
-
- bool operator>=(const mpt_wrapper& that) const {
- return m_ >= that.m_;
- }
-
- bool operator==(int that) const {
- if (that == 0)
- return m_.__get_mp()->_mp_size == 0;
- temp_[cur_] = that;
- return m_ == temp_[cur_].m_;
- }
-
- bool operator!=(int that) const {
- return !(*this == that);
- }
-
- bool operator<(int that) const {
- if (that == 0)
- return m_.__get_mp()->_mp_size < 0;
- temp_[cur_] = that;
- return m_ < temp_[cur_].m_;
- }
-
- bool operator<=(int that) const {
- if (that == 0)
- return m_.__get_mp()->_mp_size <= 0;
- temp_[cur_] = that;
- return m_ <= temp_[cur_].m_;
- }
-
- bool operator>(int that) const {
- if (that == 0)
- return m_.__get_mp()->_mp_size > 0;
- temp_[cur_] = that;
- return m_ > temp_[cur_].m_;
- }
-
- bool operator>=(int that) const {
- if (that == 0)
- return m_.__get_mp()->_mp_size >= 0;
- temp_[cur_] = that;
- return m_ >= temp_[cur_].m_;
- }
-
- mpt_wrapper& operator-() const {
- temp_[cur_].m_ = -this->m_;
- return temp_[next_cur()];
- }
-
- mpt_wrapper& operator+(const mpt_wrapper& that) const {
- temp_[cur_].m_ = this->m_ + that.m_;
- return temp_[next_cur()];
- }
-
- mpt_wrapper& operator-(const mpt_wrapper& that) const {
- temp_[cur_].m_ = this->m_ - that.m_;
- return temp_[next_cur()];
- }
-
- mpt_wrapper& operator*(const mpt_wrapper& that) const {
- temp_[cur_].m_ = this->m_ * that.m_;
- return temp_[next_cur()];
- }
-
- mpt_wrapper& operator/(const mpt_wrapper& that) const {
- temp_[cur_].m_ = this->m_ / that.m_;
- return temp_[next_cur()];
- }
-
- mpt_wrapper& operator*(double that) const {
- temp_[cur_].m_ = that;
- temp_[cur_].m_ *= this->m_;
- return temp_[next_cur()];
- }
-
- mpt_wrapper& operator+=(const mpt_wrapper& that) {
- this->m_ += that.m_;
- return *this;
- }
-
- mpt_wrapper& operator-=(const mpt_wrapper& that) {
- this->m_ -= that.m_;
- return *this;
- }
-
- mpt_wrapper& operator*=(const mpt_wrapper& that) {
- this->m_ *= that.m_;
- return *this;
- }
-
- mpt_wrapper& operator/=(const mpt_wrapper& that) {
- this->m_ /= that.m_;
- return *this;
- }
-
- mpt_wrapper& get_sqrt() {
- temp_[cur_].m_ = sqrt(m_);
- return temp_[next_cur()];
- }
-
- private:
- static int next_cur() {
- int ret_val = cur_++;
- if (cur_ == N)
- cur_ = 0;
- return ret_val;
- }
-
- mpt m_;
- static int cur_;
- static mpt_wrapper temp_[N];
- };
-
- template <typename mpt, int N>
- int mpt_wrapper<mpt, N>::cur_ = 0;
-
- template <typename mpt, int N>
- mpt_wrapper<mpt, N> mpt_wrapper<mpt, N>::temp_[N];
-
- template<int N>
- mpt_wrapper<mpf_class, N>& sqrt(mpt_wrapper<mpf_class, N>& value) {
- return value.get_sqrt();
- }
-
-} // detail
-} // polygon
-} // boost
-
-#endif
Modified: sandbox/gtl/boost/polygon/detail/voronoi_calc_utils.hpp
==============================================================================
--- sandbox/gtl/boost/polygon/detail/voronoi_calc_utils.hpp (original)
+++ sandbox/gtl/boost/polygon/detail/voronoi_calc_utils.hpp 2011-12-24 12:20:30 EST (Sat, 24 Dec 2011)
@@ -10,10 +10,6 @@
#ifndef BOOST_POLYGON_VORONOI_CALC_UTILS
#define BOOST_POLYGON_VORONOI_CALC_UTILS
-#pragma warning(disable:4800)
-#include <gmpxx.h>
-
-#include "mpt_wrapper.hpp"
#include "voronoi_robust_fpt.hpp"
namespace boost {
@@ -23,23 +19,31 @@
template <typename T>
class voronoi_calc_utils;
-// Predicate utils. Operates with the types that could
-// be converted to the int32 without precision loss.
+// Predicate utils. Operates with the coordinate types that could
+// be converted to the 32-bit signed integer without precision loss.
template <>
-class voronoi_calc_utils<int> {
+class voronoi_calc_utils<int32> {
public:
- typedef double fpt_type;
- typedef unsigned long long ulong_long_type;
+ typedef fpt64 fpt_type;
- static const unsigned int ULPS = 64;
- static const unsigned int ULPSx2 = (ULPS << 1);
+ static const uint32 ULPS;
+ static const uint32 ULPSx2;
static const fpt_type fULPS;
static const fpt_type fULPSx2;
+ // 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
+ };
+
+ // Represents orientation test result.
enum kOrientation {
RIGHT = -1,
COLLINEAR = 0,
- LEFT = 1,
+ LEFT = 1
};
// Value is a determinant of two vectors (e.g. x1 * y2 - x2 * y1).
@@ -51,19 +55,19 @@
}
// Compute robust cross_product: a1 * b2 - b1 * a2.
- // The result is correct with epsilon relative error equal to 2EPS.
+ // 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_) {
- ulong_long_type a1, b1, a2, b2;
+ uint64 a1, b1, a2, b2;
bool a1_plus, a2_plus, b1_plus, b2_plus;
a1_plus = convert_to_65_bit(a1_, a1);
b1_plus = convert_to_65_bit(b1_, b1);
a2_plus = convert_to_65_bit(a2_, a2);
b2_plus = convert_to_65_bit(b2_, b2);
- ulong_long_type expr_l = a1 * b2;
+ uint64 expr_l = a1 * b2;
bool expr_l_plus = (a1_plus == b2_plus) ? true : false;
- ulong_long_type expr_r = b1 * a2;
+ uint64 expr_r = b1 * a2;
bool expr_r_plus = (a2_plus == b1_plus) ? true : false;
if (expr_l == 0) expr_l_plus = true;
@@ -234,6 +238,7 @@
}
}
}
+
private:
typedef typename Site::point_type point_type;
@@ -424,6 +429,7 @@
}
}
}
+
private:
// Get the newer site.
const site_type &get_comparison_site(const node_type &node) const {
@@ -449,7 +455,7 @@
}
return std::make_pair(node.right_site().y(), -1);
}
- private:
+
distance_predicate_type predicate_;
};
@@ -523,15 +529,12 @@
};
template <typename Site, typename Circle>
- class gmp_circle_formation_functor {
+ class mp_circle_formation_functor {
public:
typedef typename Site::point_type point_type;
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;
+ typedef robust_sqrt_expr<eint4096, efpt64> robust_sqrt_expr_type;
void ppp(const site_type &site1,
const site_type &site2,
@@ -540,62 +543,60 @@
bool recompute_c_x = true,
bool recompute_c_y = true,
bool recompute_lower_x = true) {
- static mpt_type mpz_dif_x[3], mpz_dif_y[3], mpz_sum_x[2], mpz_sum_y[2],
- mpz_numerator[3], mpz_c_x, mpz_c_y, mpz_sqr_r, denom,
- cA[2], cB[2];
- mpz_dif_x[0] = static_cast<fpt_type>(site1.x()) -
- static_cast<fpt_type>(site2.x());
- mpz_dif_x[1] = static_cast<fpt_type>(site2.x()) -
- static_cast<fpt_type>(site3.x());
- mpz_dif_x[2] = static_cast<fpt_type>(site1.x()) -
- static_cast<fpt_type>(site3.x());
- mpz_dif_y[0] = static_cast<fpt_type>(site1.y()) -
- static_cast<fpt_type>(site2.y());
- mpz_dif_y[1] = static_cast<fpt_type>(site2.y()) -
- static_cast<fpt_type>(site3.y());
- mpz_dif_y[2] = static_cast<fpt_type>(site1.y()) -
- static_cast<fpt_type>(site3.y());
- mpz_sum_x[0] = static_cast<fpt_type>(site1.x()) +
- static_cast<fpt_type>(site2.x());
- mpz_sum_x[1] = static_cast<fpt_type>(site2.x()) +
- static_cast<fpt_type>(site3.x());
- mpz_sum_y[0] = static_cast<fpt_type>(site1.y()) +
- static_cast<fpt_type>(site2.y());
- mpz_sum_y[1] = static_cast<fpt_type>(site2.y()) +
- static_cast<fpt_type>(site3.y());
-
- denom = (mpz_dif_x[0] * mpz_dif_y[1] - mpz_dif_x[1] * mpz_dif_y[0]) * 2.0;
- mpz_numerator[1] = mpz_dif_x[0] * mpz_sum_x[0] + mpz_dif_y[0] * mpz_sum_y[0];
- mpz_numerator[2] = mpz_dif_x[1] * mpz_sum_x[1] + mpz_dif_y[1] * mpz_sum_y[1];
+ typedef eint256 eint;
+ eint dif_x[3], dif_y[3], sum_x[2], sum_y[2];
+ dif_x[0] = static_cast<int64>(site1.x()) -
+ static_cast<int64>(site2.x());
+ dif_x[1] = static_cast<int64>(site2.x()) -
+ static_cast<int64>(site3.x());
+ dif_x[2] = static_cast<int64>(site1.x()) -
+ static_cast<int64>(site3.x());
+ dif_y[0] = static_cast<int64>(site1.y()) -
+ static_cast<int64>(site2.y());
+ dif_y[1] = static_cast<int64>(site2.y()) -
+ static_cast<int64>(site3.y());
+ dif_y[2] = static_cast<int64>(site1.y()) -
+ static_cast<int64>(site3.y());
+ sum_x[0] = static_cast<int64>(site1.x()) +
+ static_cast<int64>(site2.x());
+ sum_x[1] = static_cast<int64>(site2.x()) +
+ static_cast<int64>(site3.x());
+ sum_y[0] = static_cast<int64>(site1.y()) +
+ static_cast<int64>(site2.y());
+ sum_y[1] = static_cast<int64>(site2.y()) +
+ static_cast<int64>(site3.y());
+ fpt_type denom = get_d(dif_x[0] * dif_y[1] - dif_x[1] * dif_y[0]) * 2.0;
+ eint numer1 = dif_x[0] * sum_x[0] + dif_y[0] * sum_y[0];
+ eint numer2 = dif_x[1] * sum_x[1] + dif_y[1] * sum_y[1];
if (recompute_c_x || recompute_lower_x) {
- mpz_c_x = mpz_numerator[1] * mpz_dif_y[1] - mpz_numerator[2] * mpz_dif_y[0];
- circle.x(mpz_c_x.get_d() / denom.get_d());
+ eint c_x = numer1 * dif_y[1] - numer2 * dif_y[0];
+ circle.x(get_d(c_x) / denom);
if (recompute_lower_x) {
// Evaluate radius of the circle.
- mpz_sqr_r = 1.0;
- for (int i = 0; i < 3; ++i)
- mpz_sqr_r *= mpz_dif_x[i] * mpz_dif_x[i] + mpz_dif_y[i] * mpz_dif_y[i];
- fpt_type r = std::sqrt(mpz_sqr_r.get_d());
+ eint 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(get_d(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 (circle.x() >= static_cast<fpt_type>(0.0)) {
- circle.lower_x(circle.x() + r / fabs(denom.get_d()));
+ circle.lower_x(circle.x() + r / fabs(denom));
} else {
- mpz_numerator[0] = mpz_c_x * mpz_c_x - mpz_sqr_r;
- fpt_type lower_x = mpz_numerator[0].get_d() /
- (denom.get_d() * (mpz_c_x.get_d() + r));
+ eint numer = c_x * c_x - sqr_r;
+ fpt_type lower_x = get_d(numer) /
+ (denom * (get_d(c_x) + r));
circle.lower_x(lower_x);
}
}
}
if (recompute_c_y) {
- mpz_c_y = mpz_numerator[2] * mpz_dif_x[0] - mpz_numerator[1] * mpz_dif_x[1];
- circle.y(mpz_c_y.get_d() / denom.get_d());
+ eint c_y = numer2 * dif_x[0] - numer1 * dif_x[1];
+ circle.y(get_d(c_y) / denom);
}
}
@@ -608,59 +609,60 @@
bool recompute_c_x = true,
bool recompute_c_y = true,
bool recompute_lower_x = true) {
- static mpt_type line_a, line_b, segm_len, vec_x, vec_y, sum_x, sum_y, teta,
- denom, A, B, sum_AB, dif[2], numer, cA[4], cB[4], det;
- line_a = static_cast<fpt_type>(site3.point1(true).y()) -
- static_cast<fpt_type>(site3.point0(true).y());
- line_b = static_cast<fpt_type>(site3.point0(true).x()) -
- static_cast<fpt_type>(site3.point1(true).x());
- segm_len = line_a * line_a + line_b * line_b;
- vec_x = static_cast<fpt_type>(site2.y()) -
- static_cast<fpt_type>(site1.y());
- vec_y = static_cast<fpt_type>(site1.x()) -
- static_cast<fpt_type>(site2.x());
- sum_x = static_cast<fpt_type>(site1.x()) +
- static_cast<fpt_type>(site2.x());
- sum_y = static_cast<fpt_type>(site1.y()) +
- static_cast<fpt_type>(site2.y());
- teta = line_a * vec_x + line_b * vec_y;
- denom = vec_x * line_b - vec_y * line_a;
-
- dif[0] = static_cast<fpt_type>(site3.point1().y()) -
- static_cast<fpt_type>(site1.y());
- dif[1] = static_cast<fpt_type>(site1.x()) -
- static_cast<fpt_type>(site3.point1().x());
- A = line_a * dif[1] - line_b * dif[0];
- dif[0] = static_cast<fpt_type>(site3.point1().y()) -
- static_cast<fpt_type>(site2.y());
- dif[1] = static_cast<fpt_type>(site2.x()) -
- static_cast<fpt_type>(site3.point1().x());
- B = line_a * dif[1] - line_b * dif[0];
- sum_AB = A + B;
-
- if (denom == 0) {
- numer = teta * teta - sum_AB * sum_AB;
- denom = teta * sum_AB;
+ typedef eint4096 eint;
+ eint cA[4], cB[4];
+ eint line_a = static_cast<int64>(site3.point1(true).y()) -
+ static_cast<int64>(site3.point0(true).y());
+ eint line_b = static_cast<int64>(site3.point0(true).x()) -
+ static_cast<int64>(site3.point1(true).x());
+ eint segm_len = line_a * line_a + line_b * line_b;
+ eint vec_x = static_cast<int64>(site2.y()) -
+ static_cast<int64>(site1.y());
+ eint vec_y = static_cast<int64>(site1.x()) -
+ static_cast<int64>(site2.x());
+ eint sum_x = static_cast<int64>(site1.x()) +
+ static_cast<int64>(site2.x());
+ eint sum_y = static_cast<int64>(site1.y()) +
+ static_cast<int64>(site2.y());
+ eint teta = line_a * vec_x + line_b * vec_y;
+ eint denom = vec_x * line_b - vec_y * line_a;
+
+ eint dif0 = static_cast<int64>(site3.point1().y()) -
+ static_cast<int64>(site1.y());
+ eint dif1 = static_cast<int64>(site1.x()) -
+ static_cast<int64>(site3.point1().x());
+ eint A = line_a * dif1 - line_b * dif0;
+ dif0 = static_cast<int64>(site3.point1().y()) -
+ static_cast<int64>(site2.y());
+ dif1 = static_cast<int64>(site2.x()) -
+ static_cast<int64>(site3.point1().x());
+ eint B = line_a * dif1 - line_b * dif0;
+ eint sum_AB = A + B;
+
+ if (is_zero(denom)) {
+ eint numer = teta * teta - sum_AB * sum_AB;
+ eint 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 = 1.0 / get_d(denom);
if (recompute_c_x) {
- c_event.x(0.25 * cA[0].get_d() / denom.get_d());
+ c_event.x(0.25 * get_d(cA[0]) * inv_denom);
}
if (recompute_c_y) {
- c_event.y(0.25 * cA[2].get_d() / denom.get_d());
+ c_event.y(0.25 * get_d(cA[2]) * inv_denom);
}
if (recompute_lower_x) {
- c_event.lower_x(0.25 * sqrt_expr_.eval2(cA, cB).get_d() /
- (denom.get_d() * std::sqrt(segm_len.get_d())));
+ c_event.lower_x(0.25 * get_d(sqrt_expr_.eval2(cA, cB)) * inv_denom /
+ get_sqrt(get_d(segm_len)));
}
return;
}
- det = (teta * teta + denom * denom) * A * B * 4;
- fpt_type denom_sqr = denom.get_d() * denom.get_d();
+ eint det = (teta * teta + denom * denom) * A * B * 4;
+ fpt_type inv_denom_sqr = 1.0 / sqr_value(get_d(denom));
if (recompute_c_x || recompute_lower_x) {
cA[0] = sum_x * denom * denom + teta * sum_AB * vec_x;
@@ -668,7 +670,7 @@
cA[1] = (segment_index == 2) ? -vec_x : vec_x;
cB[1] = det;
if (recompute_c_x) {
- c_event.x(0.5 * sqrt_expr_.eval2(cA, cB).get_d() / denom_sqr);
+ c_event.x(0.5 * get_d(sqrt_expr_.eval2(cA, cB)) * inv_denom_sqr);
}
}
@@ -678,20 +680,20 @@
cA[3] = (segment_index == 2) ? -vec_y : vec_y;
cB[3] = det;
if (recompute_c_y) {
- c_event.y(0.5 * sqrt_expr_.eval2(&cA[2], &cB[2]).get_d() /
- denom_sqr);
+ c_event.y(0.5 * get_d(sqrt_expr_.eval2(&cA[2], &cB[2])) *
+ inv_denom_sqr);
}
}
if (recompute_lower_x) {
- cB[0] *= segm_len;
- cB[1] *= segm_len;
+ 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(0.5 * sqrt_expr_.eval4(cA, cB).get_d() /
- (denom_sqr * std::sqrt(segm_len.get_d())));
+ c_event.lower_x(0.5 * get_d(sqrt_expr_.eval4(cA, cB)) * inv_denom_sqr /
+ get_sqrt(get_d(segm_len)));
}
}
@@ -704,88 +706,90 @@
bool recompute_c_x = true,
bool recompute_c_y = true,
bool recompute_lower_x = true) {
- static mpt_type a[2], b[2], c[2], cA[4], cB[4], orientation, dx, dy, ix, iy;
+ typedef eint4096 eint;
+ eint 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<fpt_type>(segm_end1.x()) -
- static_cast<fpt_type>(segm_start1.x());
- b[0] = static_cast<fpt_type>(segm_end1.y()) -
- static_cast<fpt_type>(segm_start1.y());
- a[1] = static_cast<fpt_type>(segm_end2.x()) -
- static_cast<fpt_type>(segm_start2.x());
- b[1] = static_cast<fpt_type>(segm_end2.y()) -
- static_cast<fpt_type>(segm_start2.y());
- orientation = a[1] * b[0] - a[0] * b[1];
- if (orientation == 0) {
- fpt_type denom = (a[0] * a[0] + b[0] * b[0]).get_d() * 2.0;
- c[0] = b[0] * (static_cast<fpt_type>(segm_start2.x()) -
- static_cast<fpt_type>(segm_start1.x())) -
- a[0] * (static_cast<fpt_type>(segm_start2.y()) -
- static_cast<fpt_type>(segm_start1.y()));
- dx = a[0] * (static_cast<fpt_type>(site1.y()) -
- static_cast<fpt_type>(segm_start1.y())) -
- b[0] * (static_cast<fpt_type>(site1.x()) -
- static_cast<fpt_type>(segm_start1.x()));
- dy = b[0] * (static_cast<fpt_type>(site1.x()) -
- static_cast<fpt_type>(segm_start2.x())) -
- a[0] * (static_cast<fpt_type>(site1.y()) -
- static_cast<fpt_type>(segm_start2.y()));
+ a[0] = static_cast<int64>(segm_end1.x()) -
+ static_cast<int64>(segm_start1.x());
+ b[0] = static_cast<int64>(segm_end1.y()) -
+ static_cast<int64>(segm_start1.y());
+ a[1] = static_cast<int64>(segm_end2.x()) -
+ static_cast<int64>(segm_start2.x());
+ b[1] = static_cast<int64>(segm_end2.y()) -
+ static_cast<int64>(segm_start2.y());
+ eint orientation = a[1] * b[0] - a[0] * b[1];
+ if (is_zero(orientation)) {
+ fpt_type denom = get_d(a[0] * a[0] + b[0] * b[0]) * 2;
+ c[0] = b[0] * (static_cast<int64>(segm_start2.x()) -
+ static_cast<int64>(segm_start1.x())) -
+ a[0] * (static_cast<int64>(segm_start2.y()) -
+ static_cast<int64>(segm_start1.y()));
+ eint dx = a[0] * (static_cast<int64>(site1.y()) -
+ static_cast<int64>(segm_start1.y())) -
+ b[0] * (static_cast<int64>(site1.x()) -
+ static_cast<int64>(segm_start1.x()));
+ eint dy = b[0] * (static_cast<int64>(site1.x()) -
+ static_cast<int64>(segm_start2.x())) -
+ a[0] * (static_cast<int64>(site1.y()) -
+ static_cast<int64>(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<fpt_type>(segm_start1.y()) +
- static_cast<fpt_type>(segm_start2.y())) -
- a[0] * b[0] * (static_cast<fpt_type>(segm_start1.x()) +
- static_cast<fpt_type>(segm_start2.x()) -
- 2.0 * static_cast<fpt_type>(site1.x())) +
- b[0] * b[0] * (2.0 * static_cast<fpt_type>(site1.y()));
- fpt_type c_y = sqrt_expr_.eval2(cA, cB).get_d();
+ cA[1] = a[0] * a[0] * (static_cast<int64>(segm_start1.y()) +
+ static_cast<int64>(segm_start2.y())) -
+ a[0] * b[0] * (static_cast<int64>(segm_start1.x()) +
+ static_cast<int64>(segm_start2.x()) -
+ static_cast<int64>(site1.x()) * 2) +
+ b[0] * b[0] * (static_cast<int64>(site1.y()) * 2);
+ fpt_type c_y = get_d(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<fpt_type>(segm_start1.x()) +
- static_cast<fpt_type>(segm_start2.x())) -
- a[0] * b[0] * (static_cast<fpt_type>(segm_start1.y()) +
- static_cast<fpt_type>(segm_start2.y()) -
- 2.0 * static_cast<fpt_type>(site1.y())) +
- a[0] * a[0] * (2.0 * static_cast<fpt_type>(site1.x()));
+ cA[1] = b[0] * b[0] * (static_cast<int64>(segm_start1.x()) +
+ static_cast<int64>(segm_start2.x())) -
+ a[0] * b[0] * (static_cast<int64>(segm_start1.y()) +
+ static_cast<int64>(segm_start2.y()) -
+ static_cast<int64>(site1.y()) * 2) +
+ a[0] * a[0] * (static_cast<int64>(site1.x()) * 2);
if (recompute_c_x) {
- fpt_type c_x = sqrt_expr_.eval2(cA, cB).get_d();
+ fpt_type c_x = get_d(sqrt_expr_.eval2(cA, cB));
c_event.x(c_x / denom);
}
if (recompute_lower_x) {
- cA[2] = (c[0] < 0) ? -c[0] : c[0];
+ cA[2] = is_neg(c[0]) ? -c[0] : c[0];
cB[2] = a[0] * a[0] + b[0] * b[0];
- fpt_type lower_x = sqrt_expr_.eval3(cA, cB).get_d();
+ fpt_type lower_x = get_d(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();
- ix = a[0] * c[1] + a[1] * c[0];
- iy = b[0] * c[1] + b[1] * c[0];
- dx = ix - orientation * site1.x();
- dy = iy - orientation * site1.y();
- if ((dx == 0) && (dy == 0)) {
- fpt_type denom = orientation.get_d();
- fpt_type c_x = ix.get_d() / denom;
- fpt_type c_y = iy.get_d() / denom;
+ c[0] = b[0] * static_cast<int32>(segm_end1.x()) -
+ a[0] * static_cast<int32>(segm_end1.y());
+ c[1] = a[1] * static_cast<int32>(segm_end2.y()) -
+ b[1] * static_cast<int32>(segm_end2.x());
+ eint ix = a[0] * c[1] + a[1] * c[0];
+ eint iy = b[0] * c[1] + b[1] * c[0];
+ eint dx = ix - orientation * static_cast<int32>(site1.x());
+ eint dy = iy - orientation * static_cast<int32>(site1.y());
+ if (is_zero(dx) && is_zero(dy)) {
+ fpt_type denom = get_d(orientation);
+ fpt_type c_x = get_d(ix) / denom;
+ fpt_type c_y = get_d(iy) / denom;
c_event = circle_type(c_x, c_y, c_x);
return;
}
- fpt_type sign = ((point_index == 2) ? 1 : -1) * ((orientation < 0) ? 1 : -1);
+ eint 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;
@@ -794,14 +798,14 @@
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 = sqrt_expr_evaluator_pss<mpt_type, mpf_type>(cA, cB).get_d();
- fpt_type denom = temp * orientation.get_d();
+ fpt_type temp = get_d(sqrt_expr_evaluator_pss<eint, efpt64>(cA, cB));
+ fpt_type denom = temp * get_d(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 = sqrt_expr_evaluator_pss<mpt_type, mpf_type>(cA, cB).get_d();
+ fpt_type cy = get_d(sqrt_expr_evaluator_pss<eint, efpt64>(cA, cB));
c_event.y(cy / denom);
}
@@ -811,13 +815,13 @@
cA[2] = ix * sign;
if (recompute_c_x) {
- fpt_type cx = sqrt_expr_evaluator_pss<mpt_type, mpf_type>(cA, cB).get_d();
+ fpt_type cx = get_d(sqrt_expr_evaluator_pss<eint, efpt64>(cA, cB));
c_event.x(cx / denom);
}
if (recompute_lower_x) {
cA[3] = orientation * (dx * dx + dy * dy) * (temp < 0 ? -1 : 1);
- fpt_type lower_x = sqrt_expr_evaluator_pss<mpt_type, mpf_type>(cA, cB).get_d();
+ fpt_type lower_x = get_d(sqrt_expr_evaluator_pss<eint, efpt64>(cA, cB));
c_event.lower_x(lower_x / denom);
}
}
@@ -831,26 +835,30 @@
bool recompute_c_x = true,
bool recompute_c_y = true,
bool recompute_lower_x = true) {
- static mpt_type a[3], b[3], c[3], cA[4], cB[4];
+ typedef eint4096 eint;
+ eint 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<fpt_type>(site1.x1(true)) -
- static_cast<fpt_type>(site1.x0(true));
- a[1] = static_cast<fpt_type>(site2.x1(true)) -
- static_cast<fpt_type>(site2.x0(true));
- a[2] = static_cast<fpt_type>(site3.x1(true)) -
- static_cast<fpt_type>(site3.x0(true));
-
- b[0] = static_cast<fpt_type>(site1.y1(true)) -
- static_cast<fpt_type>(site1.y0(true));
- b[1] = static_cast<fpt_type>(site2.y1(true)) -
- static_cast<fpt_type>(site2.y0(true));
- b[2] = static_cast<fpt_type>(site3.y1(true)) -
- static_cast<fpt_type>(site3.y0(true));
-
- c[0] = mpt_cross(site1.x0(true), site1.y0(true), site1.x1(true), site1.y1(true));
- c[1] = mpt_cross(site2.x0(true), site2.y0(true), site2.x1(true), site2.y1(true));
- c[2] = mpt_cross(site3.x0(true), site3.y0(true), site3.x1(true), site3.y1(true));
+ a[0] = static_cast<int64>(site1.x1(true)) -
+ static_cast<int64>(site1.x0(true));
+ a[1] = static_cast<int64>(site2.x1(true)) -
+ static_cast<int64>(site2.x0(true));
+ a[2] = static_cast<int64>(site3.x1(true)) -
+ static_cast<int64>(site3.x0(true));
+
+ b[0] = static_cast<int64>(site1.y1(true)) -
+ static_cast<int64>(site1.y0(true));
+ b[1] = static_cast<int64>(site2.y1(true)) -
+ static_cast<int64>(site2.y0(true));
+ b[2] = static_cast<int64>(site3.y1(true)) -
+ static_cast<int64>(site3.y0(true));
+
+ c[0] = static_cast<int64>(site1.x0(true)) * static_cast<int64>(site1.y1(true)) -
+ static_cast<int64>(site1.y0(true)) * static_cast<int64>(site1.x1(true));
+ c[1] = static_cast<int64>(site2.x0(true)) * static_cast<int64>(site2.y1(true)) -
+ static_cast<int64>(site2.y0(true)) * static_cast<int64>(site2.x1(true));
+ c[2] = static_cast<int64>(site3.x0(true)) * static_cast<int64>(site3.y1(true)) -
+ static_cast<int64>(site3.y0(true)) * static_cast<int64>(site3.x1(true));
for (int i = 0; i < 3; ++i) {
cB[i] = a[i] * a[i] + b[i] * b[i];
@@ -861,7 +869,7 @@
int k = (i+2) % 3;
cA[i] = a[j] * b[k] - a[k] * b[j];
}
- fpt_type denom = sqrt_expr_.eval3(cA, cB).get_d();
+ fpt_type denom = get_d(sqrt_expr_.eval3(cA, cB));
if (recompute_c_y) {
for (int i = 0; i < 3; ++i) {
@@ -869,7 +877,7 @@
int k = (i+2) % 3;
cA[i] = b[j] * c[k] - b[k] * c[j];
}
- fpt_type c_y = sqrt_expr_.eval3(cA, cB).get_d();
+ fpt_type c_y = get_d(sqrt_expr_.eval3(cA, cB));
c_event.y(c_y / denom);
}
@@ -880,57 +888,45 @@
int k = (i+2) % 3;
cA[i] = a[j] * c[k] - a[k] * c[j];
if (recompute_lower_x) {
- cA[3] += cA[i] * b[i];
+ cA[3] = cA[3] + cA[i] * b[i];
}
}
if (recompute_c_x) {
- fpt_type c_x = sqrt_expr_.eval3(cA, cB).get_d();
+ fpt_type c_x = get_d(sqrt_expr_.eval3(cA, cB));
c_event.x(c_x / denom);
}
if (recompute_lower_x) {
cB[3] = 1;
- fpt_type lower_x = sqrt_expr_.eval4(cA, cB).get_d();
+ fpt_type lower_x = get_d(sqrt_expr_.eval4(cA, cB));
c_event.lower_x(lower_x / denom);
}
}
}
- private:
- template <typename T>
- mpt_wrapper<mpz_class, 8> &mpt_cross(T a0, T b0, T a1, T b1) {
- static mpt_type temp[2];
- temp[0] = a0;
- temp[1] = b0;
- temp[0] = temp[0] * b1;
- temp[1] = temp[1] * a1;
- temp[0] -= temp[1];
- return temp[0];
- }
+ 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 mpt, typename mpf>
- mpf sqrt_expr_evaluator_pss(mpt *A, mpt *B) {
- static mpt cA[4], cB[4];
- static mpf lh, rh, numer;
- if (A[3] == 0) {
- lh = sqrt_expr_.eval2(A, B);
+ template <typename _int, typename _fpt>
+ _fpt sqrt_expr_evaluator_pss(_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;
- rh = A[2].get_d() * std::sqrt(B[3].get_d() *
- sqrt_expr_.eval2(cA, cB).get_d());
- if (((lh >= 0) && (rh >= 0)) || ((lh <= 0) && (rh <= 0))) {
+ _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];
- cA[0] -= A[2] * A[2] * B[3] * B[2];
+ 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];
- numer = sqrt_expr_.eval2(cA, cB);
+ _fpt numer = sqrt_expr_.eval2(cA, cB);
return numer / (lh - rh);
}
cA[0] = A[3];
@@ -939,18 +935,17 @@
cB[1] = B[0];
cA[2] = A[1];
cB[2] = B[1];
- lh = sqrt_expr_.eval3(cA, cB);
+ _fpt 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_.eval2(cA, cB).get_d());
- if (((lh >= 0) && (rh >= 0)) || ((lh <= 0) && (rh <= 0))) {
+ _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];
- cA[0] += A[3] * A[3] - A[2] * A[2] * B[2] * B[3];
+ cA[0] = A[0] * A[0] * B[0] + A[1] * A[1] * B[1] +
+ A[3] * A[3] - A[2] * A[2] * B[2] * B[3];
cB[0] = 1;
cA[1] = A[3] * A[0] * 2;
cB[1] = B[0];
@@ -958,10 +953,10 @@
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_.eval4(cA, cB);
+ _fpt numer = sqrt_expr_.eval4(cA, cB);
return numer / (lh - rh);
}
- private:
+
robust_sqrt_expr_type sqrt_expr_;
};
@@ -973,7 +968,7 @@
typedef typename Site::point_type point_type;
typedef Site site_type;
typedef Circle circle_type;
- typedef gmp_circle_formation_functor<site_type, circle_type>
+ typedef mp_circle_formation_functor<site_type, circle_type>
exact_circle_formation_functor_type;
void ppp(const site_type &site1,
@@ -1295,11 +1290,7 @@
recompute_c_x, recompute_c_y, recompute_lower_x);
}
}
- private:
- template <typename T>
- T sqr_distance(T dif_x, T dif_y) {
- return dif_x * dif_x + dif_y * dif_y;
- }
+
private:
exact_circle_formation_functor_type exact_circle_formation_functor_;
};
@@ -1379,21 +1370,34 @@
circle_existence_predicate_type circle_existence_predicate_;
circle_formation_functor_type circle_formation_functor_;
};
+
private:
// Convert value to 64-bit unsigned integer.
// Return true if the value is positive, else false.
template <typename T>
- static bool convert_to_65_bit(T value, ulong_long_type &res) {
+ static bool convert_to_65_bit(T value, uint64 &res) {
if (value >= static_cast<T>(0)) {
- res = static_cast<ulong_long_type>(value);
+ res = static_cast<uint64>(value);
return true;
} else {
- res = static_cast<ulong_long_type>(-value);
+ res = static_cast<uint64>(-value);
return false;
}
}
+
+ template <typename T>
+ static T sqr_value(T value) {
+ return value * value;
+ }
+
+ template <typename T>
+ static T sqr_distance(T dif_x, T dif_y) {
+ return dif_x * dif_x + dif_y * dif_y;
+ }
};
+const uint32 voronoi_calc_utils<int>::ULPS = 64;
+const uint32 voronoi_calc_utils<int>::ULPSx2 = 128;
const voronoi_calc_utils<int>::fpt_type voronoi_calc_utils<int>::fULPS =
voronoi_calc_utils<int>::ULPS;
const voronoi_calc_utils<int>::fpt_type voronoi_calc_utils<int>::fULPSx2 =
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 2011-12-24 12:20:30 EST (Sat, 24 Dec 2011)
@@ -10,6 +10,10 @@
#ifndef BOOST_POLYGON_VORONOI_ROBUST_FPT
#define BOOST_POLYGON_VORONOI_ROBUST_FPT
+#include <cmath>
+
+#include <boost/cstdint.hpp>
+
// Geometry predicates with floating-point variables usually require
// high-precision predicates to retrieve the correct result.
// Epsilon robust predicates give the result within some epsilon relative
@@ -52,28 +56,27 @@
namespace boost {
namespace polygon {
namespace detail {
- // 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 boost::int32_t int32;
+ typedef boost::int64_t int64;
+ typedef boost::uint32_t uint32;
+ typedef boost::uint64_t uint64;
+ typedef float fpt32;
+ typedef double fpt64;
// If two floating-point numbers in the same format are ordered (x < y),
// then they are ordered the same way when their bits are reinterpreted as
// sign-magnitude integers. Values are considered to be almost equal if
// their integer reinterpretations differ in not more than maxUlps units.
- template <typename T>
- bool almost_equal(T a, T b, unsigned int ulps);
+ template <typename _fpt>
+ bool almost_equal(_fpt a, _fpt b, uint32 ulps);
template<>
- bool almost_equal<float>(float a, float b, unsigned int maxUlps) {
- unsigned int ll_a, ll_b;
+ bool almost_equal<fpt32>(fpt32 a, fpt32 b, uint32 maxUlps) {
+ uint32 ll_a, ll_b;
// Reinterpret double bits as 32-bit signed integer.
- memcpy(&ll_a, &a, sizeof(float));
- memcpy(&ll_b, &b, sizeof(float));
+ memcpy(&ll_a, &a, sizeof(fpt32));
+ memcpy(&ll_b, &b, sizeof(fpt32));
if (ll_a < 0x80000000)
ll_a = 0x80000000 - ll_a;
@@ -86,12 +89,12 @@
}
template<>
- bool almost_equal<double>(double a, double b, unsigned int maxUlps) {
- unsigned long long ll_a, ll_b;
+ bool almost_equal<fpt64>(fpt64 a, fpt64 b, uint32 maxUlps) {
+ uint64 ll_a, ll_b;
// Reinterpret double bits as 64-bit signed integer.
- memcpy(&ll_a, &a, sizeof(double));
- memcpy(&ll_b, &b, sizeof(double));
+ memcpy(&ll_a, &a, sizeof(fpt64));
+ memcpy(&ll_b, &b, sizeof(fpt64));
// Positive 0.0 is integer zero. Negative 0.0 is 0x8000000000000000.
// Map negative zero to an integer zero representation - making it
@@ -111,36 +114,41 @@
}
template <typename T>
- double get_d(const T& value) {
- return value.get_d();
+ fpt64 get_d(const T& that) {
+ return static_cast<fpt64>(that);
}
- template <>
- double get_d(const float& value) {
- return value;
+ template <typename T>
+ T get_sqrt(const T& that) {
+ return (std::sqrt)(that);
}
- template <>
- double get_d(const double& value) {
- return value;
+ template <typename T>
+ bool is_pos(const T& that) {
+ return that > 0;
}
- template <>
- double get_d(const long double& value) {
- return static_cast<double>(value);
+ template <typename T>
+ bool is_neg(const T& that) {
+ return that < 0;
}
- template <typename FPT>
+ template <typename T>
+ bool is_zero(const T& that) {
+ return that == 0;
+ }
+
+ template <typename _fpt>
class robust_fpt {
public:
- typedef FPT floating_point_type;
- typedef FPT relative_error_type;
+ typedef _fpt floating_point_type;
+ typedef _fpt relative_error_type;
// Rounding error is at most 1 EPS.
static const relative_error_type ROUNDING_ERROR;
robust_fpt() : fpv_(0.0), re_(0.0) {}
- explicit robust_fpt(int fpv) : fpv_(fpv), re_(0.0) {}
+ explicit robust_fpt(int32 fpv) : fpv_(fpv), re_(0.0) {}
explicit robust_fpt(floating_point_type fpv,
bool rounded = true) : fpv_(fpv) {
re_ = rounded ? ROUNDING_ERROR : 0;
@@ -157,7 +165,7 @@
floating_point_type value = static_cast<floating_point_type>(that);
return almost_equal(this->fpv_,
value,
- static_cast<unsigned int>(this->ulp()));
+ static_cast<uint32>(this->ulp()));
}
template <typename T>
@@ -190,7 +198,7 @@
}
bool operator==(const robust_fpt &that) const {
- unsigned int ulp = static_cast<unsigned int>(this->re_ + that.re_);
+ uint32 ulp = static_cast<uint32>(this->re_ + that.re_);
return almost_equal(this->fpv_, that.fpv_, ulp);
}
@@ -307,7 +315,7 @@
}
robust_fpt sqrt() const {
- return robust_fpt(std::sqrt(fpv_), re_ * 0.5 + ROUNDING_ERROR);
+ return robust_fpt(get_sqrt(fpv_), re_ * 0.5 + ROUNDING_ERROR);
}
robust_fpt fabs() const {
@@ -323,6 +331,11 @@
const typename robust_fpt<T>::relative_error_type
robust_fpt<T>::ROUNDING_ERROR = 1;
+ template <typename T>
+ robust_fpt<T> get_sqrt(const robust_fpt<T>& that) {
+ return that.sqrt();
+ }
+
// robust_dif consists of two not negative values: value1 and value2.
// The resulting expression is equal to the value1 - value2.
// Substraction of a positive value is equivalent to the addition to value2
@@ -523,97 +536,770 @@
}
}
+ // Manages exponent of the floating-point value.
+ template <typename _fpt>
+ struct fpt_exponent_accessor;
+
+ template <>
+ class fpt_exponent_accessor<fpt64> {
+ public:
+ static const int64 kExponentMask;
+ static const int64 kSignedMantissaMask;
+ static const int64 kMinExponent;
+ static const int64 kMaxExponent;
+ static const int64 kMaxSignificantExpDif;
+
+ static int64 set_exponent(fpt64& value, int64 exponent) {
+ int64 bits;
+ memcpy(&bits, &value, sizeof(fpt64));
+ int64 exp = ((bits & kExponentMask) >> 52) - 1023;
+ if (exp == exponent) {
+ return exp;
+ }
+ bits = (bits & kSignedMantissaMask) | ((exponent + 1023) << 52);
+ memcpy(&value, &bits, sizeof(fpt64));
+ return exp;
+ }
+ };
+
+ const int64 fpt_exponent_accessor<fpt64>::kExponentMask = 0x7ff0000000000000LL;
+ const int64 fpt_exponent_accessor<fpt64>::kSignedMantissaMask = 0x800fffffffffffffLL;
+ const int64 fpt_exponent_accessor<fpt64>::kMinExponent = -1023LL;
+ const int64 fpt_exponent_accessor<fpt64>::kMaxExponent = 1024LL;
+ const int64 fpt_exponent_accessor<fpt64>::kMaxSignificantExpDif = 54;
+
+ // Allows to extend floating-point type exponent boundaries to the 64 bit
+ // integer range. This class does not handle division by zero, subnormal
+ // numbers or NaNs.
+ template <typename _fpt>
+ class extended_exponent_fpt {
+ public:
+ typedef _fpt fpt_type;
+ typedef int64 exp_type;
+ typedef fpt_exponent_accessor<fpt_type> fea;
+
+ explicit extended_exponent_fpt(fpt_type value) {
+ if (value == 0.0) {
+ exponent_ = 0;
+ value_ = 0.0;
+ } else {
+ exponent_ = fea::set_exponent(value, 0);
+ value_ = value;
+ }
+ }
+
+ extended_exponent_fpt(fpt_type value, exp_type exponent) {
+ if (value == 0.0) {
+ exponent_ = 0;
+ value_ = 0.0;
+ } else {
+ exponent_ = fea::set_exponent(value, 0) + exponent;
+ value_ = value;
+ }
+ }
+
+ bool is_pos() const {
+ return value_ > 0;
+ }
+
+ bool is_neg() const {
+ return value_ < 0;
+ }
+
+ bool is_zero() const {
+ return value_ == 0;
+ }
+
+ extended_exponent_fpt abs() const {
+ if (value_ >= 0) {
+ return *this;
+ } else {
+ return extended_exponent_fpt(-value_, exponent_);
+ }
+ }
+
+ extended_exponent_fpt operator-() const {
+ return extended_exponent_fpt(-value_, exponent_);
+ }
+
+ extended_exponent_fpt operator+(const extended_exponent_fpt& that) const {
+ if (this->value_ == 0.0 ||
+ that.exponent_ > this->exponent_ + fea::kMaxSignificantExpDif) {
+ return that;
+ }
+ if (that.value_ == 0.0 ||
+ this->exponent_ > that.exponent_ + fea::kMaxSignificantExpDif) {
+ return *this;
+ }
+ if (this->exponent_ >= that.exponent_) {
+ exp_type exp_dif = this->exponent_ - that.exponent_;
+ fpt_type value = this->value_;
+ fea::set_exponent(value, exp_dif);
+ value += that.value_;
+ return extended_exponent_fpt(value, that.exponent_);
+ } else {
+ exp_type exp_dif = that.exponent_ - this->exponent_;
+ fpt_type value = that.value_;
+ fea::set_exponent(value, exp_dif);
+ value += this->value_;
+ return extended_exponent_fpt(value, this->exponent_);
+ }
+ }
+
+ extended_exponent_fpt operator-(const extended_exponent_fpt& that) const {
+ if (this->value_ == 0.0 ||
+ that.exponent_ > this->exponent_ + fea::kMaxSignificantExpDif) {
+ return extended_exponent_fpt(-that.value_, that.exponent_);
+ }
+ if (that.value_ == 0.0 ||
+ this->exponent_ > that.exponent_ + fea::kMaxSignificantExpDif) {
+ return *this;
+ }
+ if (this->exponent_ >= that.exponent_) {
+ exp_type exp_dif = this->exponent_ - that.exponent_;
+ fpt_type value = this->value_;
+ fea::set_exponent(value, exp_dif);
+ value -= that.value_;
+ return extended_exponent_fpt(value, that.exponent_);
+ } else {
+ exp_type exp_dif = that.exponent_ - this->exponent_;
+ fpt_type value = -that.value_;
+ fea::set_exponent(value, exp_dif);
+ value += this->value_;
+ return extended_exponent_fpt(value, this->exponent_);
+ }
+ }
+
+ extended_exponent_fpt operator*(const extended_exponent_fpt& that) const {
+ fpt_type value = this->value_ * that.value_;
+ exp_type exponent = this->exponent_ + that.exponent_;
+ return extended_exponent_fpt(value, exponent);
+ }
+
+ extended_exponent_fpt operator/(const extended_exponent_fpt& that) const {
+ fpt_type value = this->value_ / that.value_;
+ exp_type exponent = this->exponent_ - that.exponent_;
+ return extended_exponent_fpt(value, exponent);
+ }
+
+ extended_exponent_fpt& operator+=(const extended_exponent_fpt& that) {
+ return *this = *this + that;
+ }
+
+ extended_exponent_fpt& operator-=(const extended_exponent_fpt& that) {
+ return *this = *this - that;
+ }
+
+ extended_exponent_fpt& operator*=(const extended_exponent_fpt& that) {
+ return *this = *this * that;
+ }
+
+ extended_exponent_fpt& operator/=(const extended_exponent_fpt& that) {
+ return *this = *this / that;
+ }
+
+ extended_exponent_fpt sqrt() const {
+ fpt_type val = value_;
+ exp_type exp = exponent_;
+ if (exp & 1) {
+ val *= 2.0;
+ --exp;
+ }
+ return extended_exponent_fpt(get_sqrt(val), exp >> 1);
+ }
+
+ fpt_type d() const {
+ fpt_type ret_val = value_;
+ exp_type exp = exponent_;
+ if (ret_val == 0.0) {
+ return ret_val;
+ }
+ if (exp >= fea::kMaxExponent) {
+ ret_val = 1.0;
+ exp = fea::kMaxExponent;
+ } else if (exp <= fea::kMinExponent) {
+ ret_val = 1.0;
+ exp = fea::kMinExponent;
+ }
+ fea::set_exponent(ret_val, exp);
+ return get_d(ret_val);
+ }
+
+ private:
+ fpt_type value_;
+ exp_type exponent_;
+ };
+ typedef extended_exponent_fpt<double> efpt64;
+
+ template <typename _fpt>
+ extended_exponent_fpt<_fpt> get_sqrt(const extended_exponent_fpt<_fpt>& that) {
+ return that.sqrt();
+ }
+
+ template <typename _fpt>
+ fpt64 get_d(const extended_exponent_fpt<_fpt>& that) {
+ return that.d();
+ }
+
+ template <typename _fpt>
+ bool is_pos(const extended_exponent_fpt<_fpt>& that) {
+ return that.is_pos();
+ }
+
+ template <typename _fpt>
+ bool is_neg(const extended_exponent_fpt<_fpt>& that) {
+ return that.is_neg();
+ }
+
+ template <typename _fpt>
+ bool is_zero(const extended_exponent_fpt<_fpt>& that) {
+ return that.is_zero();
+ }
+
+ template<size_t N>
+ class extended_int {
+ public:
+ static const uint64 kUInt64LowMask;
+ static const uint64 kUInt64HighMask;
+
+ extended_int() {}
+
+ extended_int(int32 that) {
+ if (that > 0) {
+ this->chunks_[0] = that;
+ this->count_ = 1;
+ } else if (that < 0) {
+ this->chunks_[0] = -that;
+ this->count_ = -1;
+ } else {
+ this->chunks_[0] = this->count_ = 0;
+ }
+ }
+
+ extended_int(int64 that) {
+ if (that > 0) {
+ this->chunks_[0] = static_cast<uint32>(that & kUInt64LowMask);
+ uint32 high = (that & kUInt64HighMask) >> 32;
+ if (high) {
+ this->chunks_[1] = high;
+ this->count_ = 2;
+ } else {
+ this->count_ = 1;
+ }
+ } else if (that < 0) {
+ that = -that;
+ this->chunks_[0] = static_cast<uint32>(that & kUInt64LowMask);
+ uint32 high = (that & kUInt64HighMask) >> 32;
+ if (high) {
+ this->chunks_[1] = high;
+ this->count_ = -2;
+ } else {
+ this->count_ = -1;
+ }
+ } else {
+ this->chunks_[0] = this->count_ = 0;
+ }
+ }
+
+ extended_int(const std::vector<uint32>& chunks, bool plus = true) {
+ this->count_ = static_cast<int32>(chunks.size());
+ for (size_t i = 0; i < chunks.size(); ++i) {
+ this->chunks_[i] = chunks[chunks.size() - i - 1];
+ }
+ if (!plus) {
+ this->count_ = -this->count_;
+ }
+ }
+
+ template <size_t M>
+ extended_int(const extended_int<M>& that) {
+ this->count_ = that.count();
+ memcpy(this->chunks_, that.chunks(), that.size() * sizeof(uint32));
+ }
+
+ extended_int<N>& operator=(int32 that) {
+ if (that > 0) {
+ this->chunks_[0] = that;
+ this->count_ = 1;
+ } else if (that < 0) {
+ this->chunks_[0] = -that;
+ this->count_ = -1;
+ } else {
+ this->chunks_[0] = this->count_ = 0;
+ }
+ return *this;
+ }
+
+ extended_int<N>& operator=(int64 that) {
+ if (that > 0) {
+ this->chunks_[0] = static_cast<uint32>(that & kUInt64LowMask);
+ uint32 high = (that & kUInt64HighMask) >> 32;
+ if (high) {
+ this->chunks_[1] = high;
+ this->count_ = 2;
+ } else {
+ this->count_ = 1;
+ }
+ } else if (that < 0) {
+ that = -that;
+ this->chunks_[0] = static_cast<uint32>(that & kUInt64LowMask);
+ uint32 high = (that & kUInt64HighMask) >> 32;
+ if (high) {
+ this->chunks_[1] = high;
+ this->count_ = -2;
+ } else {
+ this->count_ = -1;
+ }
+ } else {
+ this->chunks_[0] = this->count_ = 0;
+ }
+ return *this;
+ }
+
+ template <size_t M>
+ extended_int<N>& operator=(const extended_int<M>& that) {
+ this->count_ = that.count();
+ size_t sz = (std::min)(N, that.size());
+ memcpy(this->chunks_, that.chunks(), sz * sizeof(uint32));
+ return *this;
+ }
+
+ bool is_pos() const {
+ return this->count_ > 0;
+ }
+
+ bool is_neg() const {
+ return this->count_ < 0;
+ }
+
+ bool is_zero() const {
+ return this->count_ == 0;
+ }
+
+ template <size_t M>
+ bool operator==(const extended_int<M>& that) const {
+ if (this->count_ != that.count()) {
+ return false;
+ }
+ for (size_t i = 0; i < this->size(); ++i) {
+ if (this->chunks_[i] != that.chunks()[i]) {
+ return false;
+ }
+ }
+ return true;
+ }
+
+ template <size_t M>
+ bool operator!=(const extended_int<M>& that) const {
+ return !(*this == that);
+ }
+
+ template <size_t M>
+ bool operator<(const extended_int<M>& that) const {
+ if (this->count_ != that.count()) {
+ return this->count_ < that.count();
+ }
+ size_t i = this->size();
+ if (!i) {
+ return false;
+ }
+ do {
+ --i;
+ if (this->chunks_[i] != that.chunks()[i]) {
+ return (this->chunks_[i] < that.chunks()[i]) ^ (this->count_ < 0);
+ }
+ } while (i);
+ return false;
+ }
+
+ template <size_t M>
+ bool operator>(const extended_int<M>& that) const {
+ return that < *this;
+ }
+
+ template <size_t M>
+ bool operator<=(const extended_int<M>& that) const {
+ return !(that < *this);
+ }
+
+ template <size_t M>
+ bool operator>=(const extended_int<M>& that) const {
+ return !(*this < that);
+ }
+
+ extended_int<N> operator-() const {
+ extended_int<N> ret_val = *this;
+ ret_val.neg();
+ return ret_val;
+ }
+
+ void neg() {
+ this->count_ = -this->count_;
+ }
+
+ template <size_t M>
+ extended_int<(N>M?N:M)> operator+(const extended_int<M>& that) const {
+ extended_int<(N>M?N:M)> ret_val;
+ ret_val.add(*this, that);
+ return ret_val;
+ }
+
+ template <size_t N1, size_t N2>
+ void add(const extended_int<N1>& e1, const extended_int<N2>& e2) {
+ if (!e1.count()) {
+ *this = e2;
+ return;
+ }
+ if (!e2.count()) {
+ *this = e1;
+ return;
+ }
+ if ((e1.count() > 0) ^ (e2.count() > 0)) {
+ dif(e1.chunks(), e1.size(), e2.chunks(), e2.size());
+ } else {
+ add(e1.chunks(), e1.size(), e2.chunks(), e2.size());
+ }
+ if (e1.count() < 0) {
+ this->count_ = -this->count_;
+ }
+ }
+
+ template <size_t M>
+ extended_int<(N>M?N:M)> operator-(const extended_int<M>& that) const {
+ extended_int<(N>M?N:M)> ret_val;
+ ret_val.dif(*this, that);
+ return ret_val;
+ }
+
+ template <size_t N1, size_t N2>
+ void dif(const extended_int<N1>& e1, const extended_int<N2> &e2) {
+ if (!e1.count()) {
+ *this = e2;
+ this->count_ = -this->count_;
+ return;
+ }
+ if (!e2.count()) {
+ *this = e1;
+ return;
+ }
+ if ((e1.count() > 0) ^ (e2.count() > 0)) {
+ add(e1.chunks(), e1.size(), e2.chunks(), e2.size());
+ } else {
+ dif(e1.chunks(), e1.size(), e2.chunks(), e2.size());
+ }
+ if (e1.count() < 0) {
+ this->count_ = -this->count_;
+ }
+ }
+
+ extended_int<N> operator*(int32 that) const {
+ extended_int<N> temp(that);
+ return (*this) * temp;
+ }
+
+ extended_int<N> operator*(int64 that) const {
+ extended_int<N> temp(that);
+ return (*this) * temp;
+ }
+
+ template <size_t M>
+ extended_int<(N>M?N:M)> operator*(const extended_int<M>& that) const {
+ extended_int<(N>M?N:M)> ret_val;
+ ret_val.mul(*this, that);
+ return ret_val;
+ }
+
+ template <size_t N1, size_t N2>
+ void mul(const extended_int<N1>& e1, const extended_int<N2>& e2) {
+ if (!e1.count() || !e2.count()) {
+ this->chunks_[0] = this->count_ = 0;
+ return;
+ }
+ mul(e1.chunks(), e1.size(), e2.chunks(), e2.size());
+ if ((e1.count() > 0) ^ (e2.count() > 0)) {
+ this->count_ = -this->count_;
+ }
+ }
+
+ const uint32* chunks() const {
+ return chunks_;
+ }
+
+ int32 count() const {
+ return count_;
+ }
+
+ size_t size() const {
+ return (std::abs)(count_);
+ }
+
+ std::pair<fpt64, int64> p() const {
+ std::pair<fpt64, int64> ret_val(0, 0);
+ size_t sz = this->size();
+ if (!sz) {
+ return ret_val;
+ } else {
+ if (sz == 1) {
+ ret_val.first = static_cast<fpt64>(this->chunks_[0]);
+ ret_val.second = 0;
+ } else if (sz == 2) {
+ ret_val.first = static_cast<fpt64>(this->chunks_[1]) *
+ static_cast<fpt64>(0x100000000LL) +
+ static_cast<fpt64>(this->chunks_[0]);
+ ret_val.second = 0;
+ } else {
+ for (size_t i = 1; i <= 3; ++i) {
+ ret_val.first *= static_cast<fpt64>(0x100000000LL);
+ ret_val.first += static_cast<fpt64>(this->chunks_[sz - i]);
+ }
+ ret_val.second = (sz - 3) << 5;
+ }
+ }
+ if (this->count_ < 0) {
+ ret_val.first = -ret_val.first;
+ }
+ return ret_val;
+ }
+
+ fpt64 d() const {
+ std::pair<fpt64, int64> p = this->p();
+ extended_exponent_fpt<fpt64> efpt(p.first, p.second);
+ return efpt.d();
+ }
+
+ private:
+ void add(const uint32* c1, size_t sz1, const uint32* c2, size_t sz2) {
+ if (sz1 < sz2) {
+ add(c2, sz2, c1, sz1);
+ return;
+ }
+ this->count_ = sz1;
+ uint64 temp = 0;
+ for (size_t i = 0; i < sz2; ++i) {
+ temp += static_cast<uint64>(c1[i]) +
+ static_cast<uint64>(c2[i]);
+ this->chunks_[i] = static_cast<uint32>(temp & kUInt64LowMask);
+ temp = (temp & kUInt64HighMask) >> 32;
+ }
+ for (size_t i = sz2; i < sz1; ++i) {
+ temp += static_cast<uint64>(c1[i]);
+ this->chunks_[i] = static_cast<uint32>(temp & kUInt64LowMask);
+ temp = (temp & kUInt64HighMask) >> 32;
+ }
+ if (temp && (this->count_ != N)) {
+ this->chunks_[this->count_] = static_cast<uint32>(temp & kUInt64LowMask);
+ ++this->count_;
+ }
+ }
+
+ void dif(const uint32* c1, size_t sz1, const uint32* c2, size_t sz2, bool rec = false) {
+ if (sz1 < sz2) {
+ dif(c2, sz2, c1, sz1, true);
+ this->count_ = -this->count_;
+ return;
+ } else if ((sz1 == sz2) && !rec) {
+ do {
+ --sz1;
+ if (c1[sz1] < c2[sz1]) {
+ ++sz1;
+ dif(c2, sz1, c1, sz1, true);
+ this->count_ = -this->count_;
+ return;
+ } else if (c1[sz1] > c2[sz1]) {
+ ++sz1;
+ break;
+ } else if (!sz1) {
+ this->chunks_[0] = this->count_ = 0;
+ return;
+ }
+ } while (sz1);
+ sz2 = sz1;
+ }
+ this->count_ = sz1-1;
+ bool flag = false;
+ for (size_t i = 0; i < sz2; ++i) {
+ this->chunks_[i] = c1[i] - c2[i] - (flag?1:0);
+ flag = (c1[i] < c2[i]) || ((c1[i] == c2[i]) && flag);
+ }
+ for (size_t i = sz2; i < sz1; ++i) {
+ this->chunks_[i] = c1[i] - (flag?1:0);
+ flag = !c1[i] && flag;
+ }
+ if (this->chunks_[this->count_]) {
+ ++this->count_;
+ }
+ }
+
+ void mul(const uint32* c1, size_t sz1, const uint32* c2, size_t sz2) {
+ uint64 cur = 0, nxt, tmp;
+ this->count_ = static_cast<int32>((std::min)(N, sz1 + sz2 - 1));
+ for (size_t shift = 0; shift < static_cast<size_t>(this->count_); ++shift) {
+ nxt = 0;
+ for (size_t first = 0; first <= shift; ++first) {
+ if (first >= sz1) {
+ break;
+ }
+ size_t second = shift - first;
+ if (second >= sz2) {
+ continue;
+ }
+ tmp = static_cast<uint64>(c1[first]) *
+ static_cast<uint64>(c2[second]);
+ cur += tmp & kUInt64LowMask;
+ nxt += (tmp & kUInt64HighMask) >> 32;
+ }
+ this->chunks_[shift] = static_cast<uint32>(cur & kUInt64LowMask);
+ cur = nxt + ((cur & kUInt64HighMask) >> 32);
+ }
+ if (cur && (this->count_ != N)) {
+ this->chunks_[this->count_] = static_cast<uint32>(cur & kUInt64LowMask);
+ ++this->count_;
+ }
+ }
+
+ uint32 chunks_[N];
+ int32 count_;
+ };
+
+ template <size_t N>
+ const uint64 extended_int<N>::kUInt64LowMask = 0x00000000ffffffffULL;
+ template <size_t N>
+ const uint64 extended_int<N>::kUInt64HighMask = 0xffffffff00000000ULL;
+
+ typedef extended_int<1> eint32;
+ typedef extended_int<2> eint64;
+ typedef extended_int<4> eint128;
+ typedef extended_int<8> eint256;
+ typedef extended_int<16> eint512;
+ typedef extended_int<32> eint1024;
+ typedef extended_int<64> eint2048;
+ typedef extended_int<128> eint4096;
+ typedef extended_int<256> eint8192;
+ typedef extended_int<512> eint16384;
+ typedef extended_int<1024> eint32768;
+
+ template <size_t N>
+ fpt64 get_d(const extended_int<N>& that) {
+ return that.d();
+ }
+
+ template <size_t N>
+ bool is_pos(const extended_int<N>& that) {
+ return that.count() > 0;
+ }
+
+ template <size_t N>
+ bool is_neg(const extended_int<N>& that) {
+ return that.count() < 0;
+ }
+
+ template <size_t N>
+ bool is_zero(const extended_int<N>& that) {
+ return !that.count();
+ }
+
+
+ template <typename typeA, typename typeB>
+ struct type_converter {
+ static typeB convert(const typeA& that) {
+ return static_cast<typeB>(that);
+ }
+ };
+
+ template <typename typeA>
+ struct type_converter<typeA, fpt64> {
+ static fpt64 convert(const typeA& that) {
+ return get_d(that);
+ }
+ };
+
+ template <size_t N>
+ struct type_converter< extended_int<N>, extended_exponent_fpt<fpt64> > {
+ static extended_exponent_fpt<fpt64> convert(const extended_int<N>& that) {
+ std::pair<fpt64, int64> p = that.p();
+ return extended_exponent_fpt<fpt64>(p.first, p.second);
+ }
+ };
+
// Used to compute expressions that operate with sqrts with predefined
// relative error. Evaluates expressions of the next type:
// sum(i = 1 .. n)(A[i] * sqrt(B[i])), 1 <= n <= 4.
- template <typename mpt, typename mpf>
+ template <typename _int, typename _fpt>
class robust_sqrt_expr {
public:
+ typedef type_converter<_int, _fpt> converter;
+ static const uint32 EVAL1_MAX_RELATIVE_ERROR;
+ static const uint32 EVAL2_MAX_RELATIVE_ERROR;
+ static const uint32 EVAL3_MAX_RELATIVE_ERROR;
+ static const uint32 EVAL4_MAX_RELATIVE_ERROR;
+
// Evaluates expression (re = 4 EPS):
// A[0] * sqrt(B[0]).
- mpf& eval1(mpt *A, mpt *B) {
- a[0] = A[0];
- b[0] = B[0];
- b[0] = sqrt(b[0]);
- return b[0] *= a[0];
+ _fpt eval1(_int *A, _int *B) {
+ _fpt a = converter::convert(A[0]);
+ _fpt b = converter::convert(B[0]);
+ return a * get_sqrt(b);
}
// Evaluates expression (re = 7 EPS):
// A[0] * sqrt(B[0]) + A[1] * sqrt(B[1]).
- mpf& eval2(mpt *A, mpt *B) {
- a[1] = eval1(A, B);
- b[1] = eval1(A + 1, B + 1);
- if ((a[1] >= 0 && b[1] >= 0) || (a[1] <= 0 && b[1] <= 0))
- return b[1] += a[1];
- for (int i = 0; i < 2; ++i) {
- temp[i] = A[i] * A[i];
- temp[i] *= B[i];
- }
- a[1] -= b[1];
- b[1] = temp[0] - temp[1];
- return b[1] /= a[1];
+ _fpt eval2(_int *A, _int *B) {
+ _fpt a = eval1(A, B);
+ _fpt b = eval1(A + 1, B + 1);
+ if ((!is_neg(a) && !is_neg(b)) ||
+ (!is_pos(a) && !is_pos(b)))
+ return a + b;
+ return converter::convert(A[0] * A[0] * B[0] - A[1] * A[1] * B[1]) / (a - b);
}
// Evaluates expression (re = 16 EPS):
// A[0] * sqrt(B[0]) + A[1] * sqrt(B[1]) + A[2] * sqrt(B[2]).
- mpf& eval3(mpt *A, mpt *B) {
- a[2] = eval2(A, B);
- b[2] = eval1(A + 2, B + 2);
- if ((a[2] >= 0 && b[2] >= 0) || (a[2] <= 0 && b[2] <= 0))
- return b[2] += a[2];
- for (int i = 0; i < 3; ++i) {
- temp[i] = A[i] * A[i];
- temp[i] *= B[i];
- }
- cA[0] = temp[0] + temp[1];
- cA[0] -= temp[2];
- cB[0] = 1;
- cA[1] = A[0] * A[1];
- cA[1] += cA[1];
- cB[1] = B[0] * B[1];
- a[2] -= b[2];
- b[2] = eval2(cA, cB);
- return b[2] /= a[2];
+ _fpt eval3(_int *A, _int *B) {
+ _fpt a = eval2(A, B);
+ _fpt b = eval1(A + 2, B + 2);
+ if ((!is_neg(a) && !is_neg(b)) ||
+ (!is_pos(a) && !is_pos(b)))
+ return a + b;
+ tA[3] = A[0] * A[0] * B[0] + A[1] * A[1] * B[1] - A[2] * A[2] * B[2];
+ tB[3] = 1;
+ tA[4] = A[0] * A[1] * 2;
+ tB[4] = B[0] * B[1];
+ return eval2(tA + 3, tB + 3) / (a - b);
}
// Evaluates expression (re = 25 EPS):
// A[0] * sqrt(B[0]) + A[1] * sqrt(B[1]) +
// A[2] * sqrt(B[2]) + A[3] * sqrt(B[3]).
- mpf& eval4(mpt *A, mpt *B) {
- a[3] = eval2(A, B);
- b[3] = eval2(A + 2, B + 2);
- if ((a[3] >= 0 && b[3] >= 0) || (a[3] <= 0 && b[3] <= 0))
- return b[3] += a[3];
- for (int i = 0; i < 4; ++i) {
- temp[i] = A[i] * A[i];
- temp[i] *= B[i];
- }
- dA[0] = temp[0] + temp[1];
- dA[0] -= temp[2];
- dA[0] -= temp[3];
- dB[0] = 1;
- dA[1] = A[0] * A[1];
- dA[1] += dA[1];
- dB[1] = B[0] * B[1];
- dA[2] = A[2] * A[3];
- dA[2] = -dA[2];
- dA[2] += dA[2];
- dB[2] = B[2] * B[3];
- a[3] -= b[3];
- b[3] = eval3(dA, dB);
- return b[3] /= a[3];
+ _fpt eval4(_int *A, _int *B) {
+ _fpt a = eval2(A, B);
+ _fpt b = eval2(A + 2, B + 2);
+ if ((!is_neg(a) && !is_neg(b)) ||
+ (!is_pos(a) && !is_pos(b)))
+ return a + b;
+ tA[0] = A[0] * A[0] * B[0] + A[1] * A[1] * B[1] -
+ A[2] * A[2] * B[2] - A[3] * A[3] * B[3];
+ tB[0] = 1;
+ tA[1] = A[0] * A[1] * 2;
+ tB[1] = B[0] * B[1];
+ tA[2] = A[2] * A[3] * -2;
+ tB[2] = B[2] * B[3];
+ return eval3(tA, tB) / (a - b);
}
private:
- mpf a[4];
- mpf b[4];
- mpt cA[2];
- mpt cB[2];
- mpt dA[3];
- mpt dB[3];
- mpt temp[4];
+ _int tA[5];
+ _int tB[5];
};
+
+ template <typename _int, typename _fpt>
+ const uint32 robust_sqrt_expr<_int, _fpt>::EVAL1_MAX_RELATIVE_ERROR = 4;
+ template <typename _int, typename _fpt>
+ const uint32 robust_sqrt_expr<_int, _fpt>::EVAL2_MAX_RELATIVE_ERROR = 7;
+ template <typename _int, typename _fpt>
+ const uint32 robust_sqrt_expr<_int, _fpt>::EVAL3_MAX_RELATIVE_ERROR = 16;
+ template <typename _int, typename _fpt>
+ const uint32 robust_sqrt_expr<_int, _fpt>::EVAL4_MAX_RELATIVE_ERROR = 25;
} // detail
} // polygon
} // boost
Modified: sandbox/gtl/boost/polygon/directed_line_segment_set_data.hpp
==============================================================================
--- sandbox/gtl/boost/polygon/directed_line_segment_set_data.hpp (original)
+++ sandbox/gtl/boost/polygon/directed_line_segment_set_data.hpp 2011-12-24 12:20:30 EST (Sat, 24 Dec 2011)
@@ -257,7 +257,7 @@
typename std::vector<Point>::iterator new_end = std::unique(tmp_points.begin(), tmp_points.end());
output_points.insert(output_points.end(), tmp_points.begin(), new_end);
return std::distance(tmp_points.begin(), new_end);
- };
+ }
private:
Modified: sandbox/gtl/boost/polygon/voronoi_diagram.hpp
==============================================================================
--- sandbox/gtl/boost/polygon/voronoi_diagram.hpp (original)
+++ sandbox/gtl/boost/polygon/voronoi_diagram.hpp 2011-12-24 12:20:30 EST (Sat, 24 Dec 2011)
@@ -145,7 +145,7 @@
enum kEdgeType {
SEGMENT = 0,
RAY = 1,
- LINE = 2,
+ LINE = 2
};
// Get a view rectangle based on the voronoi bounding rectangle.
Modified: sandbox/gtl/libs/polygon/Jamroot
==============================================================================
--- sandbox/gtl/libs/polygon/Jamroot (original)
+++ sandbox/gtl/libs/polygon/Jamroot 2011-12-24 12:20:30 EST (Sat, 24 Dec 2011)
@@ -8,39 +8,16 @@
import os ;
path-constant BOOST_ROOT : [ os.environ BOOST_ROOT ] ;
-path-constant GMP_PATH : "c:/Users/Slevin/SWE/Sweepline/gmp" ;
-#path-constant GMP_PATH : "/usr/local" ;
-
-lib gmp
- : # sources
- : # requirements
- <name>gmp <search>$(GMP_PATH)/lib
- : #default-build
- : # usage-requirements
- ;
-
-lib gmpxx
- : # sources
- : # requirements
- <name>gmpxx <search>$(GMP_PATH)/lib
- : # default-build
- : # usage-requirements
- ;
-
-alias gmp_libs : gmp gmpxx : <link>static ;
-
project
: requirements
-# <warnings>all
-# <toolset>intel:<warnings>on
-# <toolset>gcc:<cxxflags>"-pedantic -Wall -Wstrict-aliasing -fstrict-aliasing -Wno-long-long"
-# <toolset>msvc:<cxxflags>/W4
+ <warnings>all
+ <toolset>intel:<warnings>on
+ <toolset>gcc:<cxxflags>"-pedantic -Wall -Wstrict-aliasing -fstrict-aliasing -Wno-long-long"
+ <toolset>msvc:<cxxflags>/W4
<include>../..
<include>.
<include>$(BOOST_ROOT)
- <include>$(GMP_PATH)/include
- <library>gmp_libs
:
build-dir bin.v2
;
Modified: sandbox/gtl/libs/polygon/test/voronoi_calc_utils_test.cpp
==============================================================================
--- sandbox/gtl/libs/polygon/test/voronoi_calc_utils_test.cpp (original)
+++ sandbox/gtl/libs/polygon/test/voronoi_calc_utils_test.cpp 2011-12-24 12:20:30 EST (Sat, 24 Dec 2011)
@@ -31,9 +31,9 @@
node_comparison_type node_comparison;
typedef VCU::circle_existence_predicate<site_type> CEP_type;
-typedef VCU::gmp_circle_formation_functor<site_type, circle_type> GMP_CFF_type;
+typedef VCU::mp_circle_formation_functor<site_type, circle_type> MP_CFF_type;
typedef VCU::lazy_circle_formation_functor<site_type, circle_type> lazy_CFF_type;
-VCU::circle_formation_predicate<site_type, circle_type, CEP_type, GMP_CFF_type> gmp_predicate;
+VCU::circle_formation_predicate<site_type, circle_type, CEP_type, MP_CFF_type> mp_predicate;
VCU::circle_formation_predicate<site_type, circle_type, CEP_type, lazy_CFF_type> lazy_predicate;
#define CHECK_ORIENTATION(P1, P2, P3, R1, R2) \
@@ -68,7 +68,7 @@
#define CHECK_CIRCLE_FORMATION_PREDICATE(s1, s2, s3, c_x, c_y, l_x) \
{ circle_type c1, c2; \
- BOOST_CHECK_EQUAL(gmp_predicate(s1, s2, s3, c1), true); \
+ BOOST_CHECK_EQUAL(mp_predicate(s1, s2, s3, c1), true); \
BOOST_CHECK_EQUAL(lazy_predicate(s1, s2, s3, c2), true); \
CHECK_CIRCLE(c1, c_x, c_y, l_x); \
CHECK_CIRCLE(c2, c_x, c_y, l_x); }
@@ -397,7 +397,7 @@
CHECK_CIRCLE_EXISTENCE(site4, site5, site3, false);
}
-BOOST_AUTO_TEST_CASE(circle_formation_predicate_test5) {
+BOOST_AUTO_TEST_CASE(circle_formation_predicate_test4) {
site_type site1(-4, 0, -4, 20);
site_type site2(-2, 10);
site_type site3(4, 10);
@@ -405,7 +405,7 @@
CHECK_CIRCLE_FORMATION_PREDICATE(site3, site2, site1, 1.0, 14.0, 6.0);
}
-BOOST_AUTO_TEST_CASE(circle_formation_predicate_test6) {
+BOOST_AUTO_TEST_CASE(circle_formation_predicate_test5) {
site_type site1(1, 0, 7, 0);
site1.inverse();
site_type site2(-2, 4, 10, 4);
@@ -415,7 +415,7 @@
CHECK_CIRCLE_FORMATION_PREDICATE(site4, site2, site1, 1.0, 2.0, 3.0);
}
-BOOST_AUTO_TEST_CASE(circle_formation_predicate_test7) {
+BOOST_AUTO_TEST_CASE(circle_formation_predicate_test6) {
site_type site1(-1, 2, 8, -10);
site1.inverse();
site_type site2(-1, 0, 8, 12);
@@ -423,7 +423,7 @@
CHECK_CIRCLE_FORMATION_PREDICATE(site3, site2, site1, 6.0, 1.0, 11.0);
}
-BOOST_AUTO_TEST_CASE(circle_formation_predicate_test8) {
+BOOST_AUTO_TEST_CASE(circle_formation_predicate_test7) {
site_type site1(1, 0, 6, 0);
site1.inverse();
site_type site2(-6, 4, 0, 12);
@@ -431,7 +431,7 @@
CHECK_CIRCLE_FORMATION_PREDICATE(site3, site2, site1, 1.0, 5.0, 6.0);
}
-BOOST_AUTO_TEST_CASE(circle_formation_predicate_test9) {
+BOOST_AUTO_TEST_CASE(circle_formation_predicate_test8) {
site_type site1(1, 0, 5, 0);
site1.inverse();
site_type site2(0, 12, 8, 6);
@@ -439,7 +439,7 @@
CHECK_CIRCLE_FORMATION_PREDICATE(site3, site2, site1, 1.0, 5.0, 6.0);
}
-BOOST_AUTO_TEST_CASE(circle_formation_predicate_test10) {
+BOOST_AUTO_TEST_CASE(circle_formation_predicate_test9) {
site_type site1(0, 0, 4, 0);
site_type site2(0, 0, 0, 4);
site_type site3(0, 4, 4, 4);
@@ -447,7 +447,7 @@
CHECK_CIRCLE_FORMATION_PREDICATE(site1, site2, site3, 2.0, 2.0, 4.0);
}
-BOOST_AUTO_TEST_CASE(circle_formation_predicate_test11) {
+BOOST_AUTO_TEST_CASE(circle_formation_predicate_test10) {
site_type site1(1, 0, 41, 30);
site_type site2(-39, 30, 1, 60);
site_type site3(1, 60, 41, 30);
Modified: sandbox/gtl/libs/polygon/test/voronoi_robust_fpt_test.cpp
==============================================================================
--- sandbox/gtl/libs/polygon/test/voronoi_robust_fpt_test.cpp (original)
+++ sandbox/gtl/libs/polygon/test/voronoi_robust_fpt_test.cpp 2011-12-24 12:20:30 EST (Sat, 24 Dec 2011)
@@ -17,126 +17,125 @@
#include <boost/test/test_case_template.hpp>
#include "boost/polygon/detail/voronoi_robust_fpt.hpp"
-#include "boost/polygon/detail/voronoi_calc_utils.hpp"
using namespace boost::polygon::detail;
typedef robust_fpt<double> rfpt_type;
BOOST_AUTO_TEST_CASE(robust_fpt_constructors_test1) {
rfpt_type a = rfpt_type();
- BOOST_CHECK_EQUAL(a.fpv() == 0.0, true);
- BOOST_CHECK_EQUAL(a.re() == 0.0, true);
- BOOST_CHECK_EQUAL(a.ulp() == 0, true);
+ BOOST_CHECK_EQUAL(a.fpv(),0.0);
+ BOOST_CHECK_EQUAL(a.re(), 0.0);
+ BOOST_CHECK_EQUAL(a.ulp(), 0);
}
BOOST_AUTO_TEST_CASE(robust_fpt_constructors_test2) {
rfpt_type a(10.0);
- BOOST_CHECK_EQUAL(a.fpv() == 10.0, true);
- BOOST_CHECK_EQUAL(a.re() == 1.0, true);
- BOOST_CHECK_EQUAL(a.ulp() == 1.0, true);
+ BOOST_CHECK_EQUAL(a.fpv(), 10.0);
+ BOOST_CHECK_EQUAL(a.re(), 1.0);
+ BOOST_CHECK_EQUAL(a.ulp(), 1.0);
}
BOOST_AUTO_TEST_CASE(robust_fpt_constructors_test3) {
rfpt_type a(10.0, false);
- BOOST_CHECK_EQUAL(a.fpv() == 10.0, true);
- BOOST_CHECK_EQUAL(a.re() == 0.0, true);
- BOOST_CHECK_EQUAL(a.ulp() == 0.0, true);
+ BOOST_CHECK_EQUAL(a.fpv(), 10.0);
+ BOOST_CHECK_EQUAL(a.re(), 0.0);
+ BOOST_CHECK_EQUAL(a.ulp(), 0.0);
}
BOOST_AUTO_TEST_CASE(robust_fpt_constructors_test4) {
rfpt_type a(10.0, 3.0);
- BOOST_CHECK_EQUAL(a.fpv() == 10.0, true);
- BOOST_CHECK_EQUAL(a.re() == 3.0, true);
- BOOST_CHECK_EQUAL(a.ulp() == 3.0, true);
+ BOOST_CHECK_EQUAL(a.fpv(), 10.0);
+ BOOST_CHECK_EQUAL(a.re(), 3.0);
+ BOOST_CHECK_EQUAL(a.ulp(), 3.0);
rfpt_type b(10.0, 2.75);
- BOOST_CHECK_EQUAL(b.fpv() == 10.0, true);
- BOOST_CHECK_EQUAL(b.re() == 2.75, true);
- BOOST_CHECK_EQUAL(b.ulp() == 2.75, true);
+ BOOST_CHECK_EQUAL(b.fpv(), 10.0);
+ BOOST_CHECK_EQUAL(b.re(), 2.75);
+ BOOST_CHECK_EQUAL(b.ulp(), 2.75);
}
BOOST_AUTO_TEST_CASE(robust_fpt_sum_test1) {
rfpt_type a(2.0, 5.0);
rfpt_type b(3.0, 4.0);
rfpt_type c = a + b;
- BOOST_CHECK_EQUAL(c.fpv() == 5.0, true);
- BOOST_CHECK_EQUAL(c.re() == 6.0, true);
- BOOST_CHECK_EQUAL(c.ulp() == 6.0, true);
+ BOOST_CHECK_EQUAL(c.fpv(), 5.0);
+ BOOST_CHECK_EQUAL(c.re(), 6.0);
+ BOOST_CHECK_EQUAL(c.ulp(), 6.0);
c += b;
- BOOST_CHECK_EQUAL(c.fpv() == 8.0, true);
- BOOST_CHECK_EQUAL(c.re() == 7.0, true);
- BOOST_CHECK_EQUAL(c.ulp() == 7.0, true);
+ BOOST_CHECK_EQUAL(c.fpv(), 8.0);
+ BOOST_CHECK_EQUAL(c.re(), 7.0);
+ BOOST_CHECK_EQUAL(c.ulp(), 7.0);
}
BOOST_AUTO_TEST_CASE(robust_fpt_sum_test2) {
rfpt_type a(3.0, 2.0);
rfpt_type b(-2.0, 3.0);
rfpt_type c = a + b;
- BOOST_CHECK_EQUAL(c.fpv() == 1.0, true);
- BOOST_CHECK_EQUAL(c.re() == 13.0, true);
- BOOST_CHECK_EQUAL(c.ulp() == 13.0, true);
+ BOOST_CHECK_EQUAL(c.fpv(), 1.0);
+ BOOST_CHECK_EQUAL(c.re(), 13.0);
+ BOOST_CHECK_EQUAL(c.ulp(), 13.0);
c += b;
- BOOST_CHECK_EQUAL(c.fpv() == -1.0, true);
- BOOST_CHECK_EQUAL(c.re() == 20.0, true);
- BOOST_CHECK_EQUAL(c.ulp() == 20.0, true);
+ BOOST_CHECK_EQUAL(c.fpv(), -1.0);
+ BOOST_CHECK_EQUAL(c.re(), 20.0);
+ BOOST_CHECK_EQUAL(c.ulp(), 20.0);
}
BOOST_AUTO_TEST_CASE(robust_fpt_dif_test1) {
rfpt_type a(2.0, 5.0);
rfpt_type b(-3.0, 4.0);
rfpt_type c = a - b;
- BOOST_CHECK_EQUAL(c.fpv() == 5.0, true);
- BOOST_CHECK_EQUAL(c.re() == 6.0, true);
- BOOST_CHECK_EQUAL(c.ulp() == 6.0, true);
+ BOOST_CHECK_EQUAL(c.fpv(), 5.0);
+ BOOST_CHECK_EQUAL(c.re(), 6.0);
+ BOOST_CHECK_EQUAL(c.ulp(), 6.0);
c -= b;
- BOOST_CHECK_EQUAL(c.fpv() == 8.0, true);
- BOOST_CHECK_EQUAL(c.re() == 7.0, true);
- BOOST_CHECK_EQUAL(c.ulp() == 7.0, true);
+ BOOST_CHECK_EQUAL(c.fpv(), 8.0);
+ BOOST_CHECK_EQUAL(c.re(), 7.0);
+ BOOST_CHECK_EQUAL(c.ulp(), 7.0);
}
BOOST_AUTO_TEST_CASE(robust_fpt_dif_test2) {
rfpt_type a(3.0, 2.0);
rfpt_type b(2.0, 3.0);
rfpt_type c = a - b;
- BOOST_CHECK_EQUAL(c.fpv() == 1.0, true);
- BOOST_CHECK_EQUAL(c.re() == 13.0, true);
- BOOST_CHECK_EQUAL(c.ulp() == 13.0, true);
+ BOOST_CHECK_EQUAL(c.fpv(), 1.0);
+ BOOST_CHECK_EQUAL(c.re(), 13.0);
+ BOOST_CHECK_EQUAL(c.ulp(), 13.0);
c -= b;
- BOOST_CHECK_EQUAL(c.fpv() == -1.0, true);
- BOOST_CHECK_EQUAL(c.re() == 20.0, true);
- BOOST_CHECK_EQUAL(c.ulp() == 20.0, true);
+ BOOST_CHECK_EQUAL(c.fpv(), -1.0);
+ BOOST_CHECK_EQUAL(c.re(), 20.0);
+ BOOST_CHECK_EQUAL(c.ulp(), 20.0);
}
BOOST_AUTO_TEST_CASE(robust_fpt_mult_test3) {
rfpt_type a(2.0, 3.0);
rfpt_type b(4.0);
rfpt_type c = a * b;
- BOOST_CHECK_EQUAL(c.fpv() == 8.0, true);
- BOOST_CHECK_EQUAL(c.re() == 5.0, true);
- BOOST_CHECK_EQUAL(c.ulp() == 5.0, true);
+ BOOST_CHECK_EQUAL(c.fpv(), 8.0);
+ BOOST_CHECK_EQUAL(c.re(), 5.0);
+ BOOST_CHECK_EQUAL(c.ulp(), 5.0);
c *= b;
- BOOST_CHECK_EQUAL(c.fpv() == 32.0, true);
- BOOST_CHECK_EQUAL(c.re() == 7.0, true);
- BOOST_CHECK_EQUAL(c.ulp() == 7.0, true);
+ BOOST_CHECK_EQUAL(c.fpv(), 32.0);
+ BOOST_CHECK_EQUAL(c.re(), 7.0);
+ BOOST_CHECK_EQUAL(c.ulp(), 7.0);
}
BOOST_AUTO_TEST_CASE(robust_fpt_div_test1) {
rfpt_type a(2.0, 3.0);
rfpt_type b(4.0);
rfpt_type c = a / b;
- BOOST_CHECK_EQUAL(c.fpv() == 0.5, true);
- BOOST_CHECK_EQUAL(c.re() == 5.0, true);
- BOOST_CHECK_EQUAL(c.ulp() == 5.0, true);
+ BOOST_CHECK_EQUAL(c.fpv(), 0.5);
+ BOOST_CHECK_EQUAL(c.re(), 5.0);
+ BOOST_CHECK_EQUAL(c.ulp(), 5.0);
c /= b;
- BOOST_CHECK_EQUAL(c.fpv() == 0.125, true);
- BOOST_CHECK_EQUAL(c.re() == 7.0, true);
- BOOST_CHECK_EQUAL(c.ulp() == 7.0, true);
+ BOOST_CHECK_EQUAL(c.fpv(), 0.125);
+ BOOST_CHECK_EQUAL(c.re(), 7.0);
+ BOOST_CHECK_EQUAL(c.ulp(), 7.0);
}
BOOST_AUTO_TEST_CASE(robust_dif_constructors_test) {
@@ -225,52 +224,361 @@
BOOST_CHECK_EQUAL(c4[3].dif() / b4[3].dif(), a4[3].dif());
}
-template <typename mpz, typename mpf>
+BOOST_AUTO_TEST_CASE(fpt_exponent_accessor_test) {
+ typedef fpt_exponent_accessor<fpt64> fpt_ea;
+ fpt64 value = 15;
+ BOOST_CHECK_EQUAL(fpt_ea::set_exponent(value, 0), 3);
+ BOOST_CHECK_EQUAL(value, 1.875);
+ value = 0.0625;
+ BOOST_CHECK_EQUAL(fpt_ea::set_exponent(value, 0), -4);
+ BOOST_CHECK_EQUAL(value, 1.0);
+ value = -1.5;
+ BOOST_CHECK_EQUAL(fpt_ea::set_exponent(value, 4), 0);
+ BOOST_CHECK_EQUAL(value, -24.0);
+ value = 0.0;
+ BOOST_CHECK_EQUAL(fpt_ea::set_exponent(value, 4), fpt_ea::kMinExponent);
+ BOOST_CHECK_EQUAL(value, 16.0);
+ value = std::pow(2.0, 2000);
+ BOOST_CHECK_EQUAL(fpt_ea::set_exponent(value, 4), fpt_ea::kMaxExponent);
+ BOOST_CHECK_EQUAL(value, 16.0);
+}
+
+BOOST_AUTO_TEST_CASE(extended_exponent_fpt_test1) {
+ boost::mt19937_64 gen(static_cast<uint32>(time(NULL)));
+ fpt64 b = 0.0;
+ efpt64 eeb(b);
+ for (int i = 0; i < 1000; ++i) {
+ fpt64 a = static_cast<fpt64>(static_cast<int64>(gen()));
+ efpt64 eea(a);
+ efpt64 neg = -eea;
+ efpt64 sum = eea + eeb;
+ efpt64 dif = eea - eeb;
+ efpt64 mul = eea * eeb;
+ BOOST_CHECK_EQUAL(get_d(neg), -a);
+ BOOST_CHECK_EQUAL(get_d(sum), a + b);
+ BOOST_CHECK_EQUAL(get_d(dif), a - b);
+ BOOST_CHECK_EQUAL(get_d(mul), a * b);
+ }
+}
+
+BOOST_AUTO_TEST_CASE(extended_exponent_fpt_test2) {
+ boost::mt19937_64 gen(static_cast<uint32>(time(NULL)));
+ fpt64 a = 0.0;
+ efpt64 eea(a);
+ for (int i = 0; i < 1000; ++i) {
+ fpt64 b = static_cast<fpt64>(static_cast<int64>(gen()));
+ if (b == 0.0) {
+ continue;
+ }
+ efpt64 eeb(b);
+ efpt64 neg = -eea;
+ efpt64 sum = eea + eeb;
+ efpt64 dif = eea - eeb;
+ efpt64 mul = eea * eeb;
+ efpt64 div = eea / eeb;
+ BOOST_CHECK_EQUAL(get_d(neg), -a);
+ BOOST_CHECK_EQUAL(get_d(sum), a + b);
+ BOOST_CHECK_EQUAL(get_d(dif), a - b);
+ BOOST_CHECK_EQUAL(get_d(mul), a * b);
+ BOOST_CHECK_EQUAL(get_d(div), a / b);
+ }
+}
+
+BOOST_AUTO_TEST_CASE(extended_exponent_fpt_test3) {
+ boost::mt19937_64 gen(static_cast<uint32>(time(NULL)));
+ for (int i = 0; i < 1000; ++i) {
+ fpt64 a = static_cast<fpt64>(static_cast<int64>(gen()));
+ fpt64 b = static_cast<fpt64>(static_cast<int64>(gen()));
+ if (b == 0.0) {
+ continue;
+ }
+ efpt64 eea(a);
+ efpt64 eeb(b);
+ efpt64 neg = -eea;
+ efpt64 sum = eea + eeb;
+ efpt64 dif = eea - eeb;
+ efpt64 mul = eea * eeb;
+ efpt64 div = eea / eeb;
+ BOOST_CHECK_EQUAL(get_d(neg), -a);
+ BOOST_CHECK_EQUAL(get_d(sum), a + b);
+ BOOST_CHECK_EQUAL(get_d(dif), a - b);
+ BOOST_CHECK_EQUAL(get_d(mul), a * b);
+ BOOST_CHECK_EQUAL(get_d(div), a / b);
+ }
+}
+
+BOOST_AUTO_TEST_CASE(extended_exponent_fpt_test4) {
+ for (int exp = 0; exp < 64; ++exp)
+ for (int i = 1; i < 100; ++i) {
+ fpt64 a = i;
+ fpt64 b = static_cast<fpt64>(1LL << exp);
+ efpt64 eea(a);
+ efpt64 eeb(b);
+ efpt64 neg = -eea;
+ efpt64 sum = eea + eeb;
+ efpt64 dif = eea - eeb;
+ efpt64 mul = eea * eeb;
+ efpt64 div = eea / eeb;
+ BOOST_CHECK_EQUAL(get_d(neg), -a);
+ BOOST_CHECK_EQUAL(get_d(sum), a + b);
+ BOOST_CHECK_EQUAL(get_d(dif), a - b);
+ BOOST_CHECK_EQUAL(get_d(mul), a * b);
+ BOOST_CHECK_EQUAL(get_d(div), a / b);
+ }
+}
+
+BOOST_AUTO_TEST_CASE(extended_exponent_fpt_test5) {
+ for (int i = 0; i < 100; ++i) {
+ efpt64 a(static_cast<fpt64>(i * i));
+ efpt64 b = a.sqrt();
+ BOOST_CHECK_EQUAL(get_d(b), static_cast<fpt64>(i));
+ }
+}
+
+BOOST_AUTO_TEST_CASE(extended_int_test1) {
+ eint32 e1(0), e2(32), e3(-32);
+ BOOST_CHECK_EQUAL(e1.count(), 0);
+ BOOST_CHECK_EQUAL(e1.chunks()[0], 0);
+ BOOST_CHECK_EQUAL(e1.size(), 0);
+ BOOST_CHECK_EQUAL(e2.count(), 1);
+ BOOST_CHECK_EQUAL(e2.chunks()[0], 32);
+ BOOST_CHECK_EQUAL(e2.size(), 1);
+ BOOST_CHECK_EQUAL(e3.count(), -1);
+ BOOST_CHECK_EQUAL(e3.chunks()[0], 32);
+ BOOST_CHECK_EQUAL(e3.size(), 1);
+}
+
+BOOST_AUTO_TEST_CASE(extended_int_test2) {
+ int64 val64 = 0x7fffffffffffffffLL;
+ eint64 e1(0LL), e2(32LL), e3(-32LL), e4(val64), e5(-val64);
+ BOOST_CHECK_EQUAL(e1.count(), 0);
+ BOOST_CHECK_EQUAL(e1.chunks()[0], 0);
+ BOOST_CHECK_EQUAL(e2.count(), 1);
+ BOOST_CHECK_EQUAL(e2.chunks()[0], 32);
+ BOOST_CHECK_EQUAL(e3.count(), -1);
+ BOOST_CHECK_EQUAL(e3.chunks()[0], 32);
+ BOOST_CHECK_EQUAL(e4.count(), 2);
+ BOOST_CHECK_EQUAL(e4.chunks()[0], 0xffffffff);
+ BOOST_CHECK_EQUAL(e4.chunks()[1], val64 >> 32);
+ BOOST_CHECK_EQUAL(e5.count(), -2);
+ BOOST_CHECK_EQUAL(e5.chunks()[0], 0xffffffff);
+ BOOST_CHECK_EQUAL(e5.chunks()[1], val64 >> 32);
+}
+
+BOOST_AUTO_TEST_CASE(extended_int_test3) {
+ std::vector<uint32> chunks;
+ chunks.push_back(1);
+ chunks.push_back(2);
+ eint64 e1(chunks, true), e2(chunks, false);
+ BOOST_CHECK_EQUAL(e1.count(), 2);
+ BOOST_CHECK_EQUAL(e1.chunks()[0], 2);
+ BOOST_CHECK_EQUAL(e1.chunks()[1], 1);
+ BOOST_CHECK_EQUAL(e2.count(), -2);
+ BOOST_CHECK_EQUAL(e2.chunks()[0], 2);
+ BOOST_CHECK_EQUAL(e2.chunks()[1], 1);
+}
+
+BOOST_AUTO_TEST_CASE(extended_int_test4) {
+ std::vector<uint32> chunks;
+ chunks.push_back(1);
+ chunks.push_back(2);
+ eint64 e1(chunks, true), e2(chunks, false);
+ BOOST_CHECK_EQUAL(e1 == e2, false);
+ BOOST_CHECK_EQUAL(e1 == -e2, true);
+ BOOST_CHECK_EQUAL(e1 != e2, true);
+ BOOST_CHECK_EQUAL(e1 != -e2, false);
+ BOOST_CHECK_EQUAL(e1 < e2, false);
+ BOOST_CHECK_EQUAL(e1 < -e2, false);
+ BOOST_CHECK_EQUAL(e1 <= e2, false);
+ BOOST_CHECK_EQUAL(e1 <= -e2, true);
+ BOOST_CHECK_EQUAL(e1 > e2, true);
+ BOOST_CHECK_EQUAL(e1 > -e2, false);
+ BOOST_CHECK_EQUAL(e1 >= e2, true);
+ BOOST_CHECK_EQUAL(e1 >= -e2, true);
+}
+
+BOOST_AUTO_TEST_CASE(extended_int_test5) {
+ boost::mt19937_64 gen(static_cast<uint32>(time(NULL)));
+ for (int i = 0; i < 100; ++i) {
+ int64 i1 = static_cast<int64>(gen());
+ int64 i2 = static_cast<int64>(gen());
+ eint64 e1(i1), e2(i2);
+ BOOST_CHECK_EQUAL(e1 == e2, i1 == i2);
+ BOOST_CHECK_EQUAL(e1 != e2, i1 != i2);
+ BOOST_CHECK_EQUAL(e1 > e2, i1 > i2);
+ BOOST_CHECK_EQUAL(e1 >= e2, i1 >= i2);
+ BOOST_CHECK_EQUAL(e1 < e2, i1 < i2);
+ BOOST_CHECK_EQUAL(e1 <= e2, i1 <= i2);
+ }
+}
+
+BOOST_AUTO_TEST_CASE(extended_int_test6) {
+ eint32 e1(32);
+ eint32 e2 = -e1;
+ BOOST_CHECK_EQUAL(e2.count(), -1);
+ BOOST_CHECK_EQUAL(e2.size(), 1);
+ BOOST_CHECK_EQUAL(e2.chunks()[0], 32);
+}
+
+BOOST_AUTO_TEST_CASE(extended_int_test7) {
+ boost::mt19937_64 gen(static_cast<uint32>(time(NULL)));
+ for (int i = 0; i < 100; ++i) {
+ int64 i1 = static_cast<int64>(gen()) >> 2;
+ int64 i2 = static_cast<int64>(gen()) >> 2;
+ eint64 e1(i1), e2(i2), e3(i1 + i2), e4(i1 - i2);
+ BOOST_CHECK(e1 + e2 == e3);
+ BOOST_CHECK(e1 - e2 == e4);
+ }
+}
+
+BOOST_AUTO_TEST_CASE(extended_int_test8) {
+ boost::mt19937 gen(static_cast<uint32>(time(NULL)));
+ for (int i = 0; i < 100; ++i) {
+ int64 i1 = static_cast<int32>(gen());
+ int64 i2 = static_cast<int32>(gen());
+ eint64 e1(i1), e2(i2), e3(i1 * i2);
+ BOOST_CHECK(e1 * e2 == e3);
+ }
+}
+
+BOOST_AUTO_TEST_CASE(exnteded_int_test9) {
+ for (int i = -10; i <= 10; ++i) {
+ for (int j = -10; j <= 10; ++j) {
+ eint32 e1(i), e2(j), e3(i+j), e4(i-j), e5(i*j);
+ BOOST_CHECK(e1 + e2 == e3);
+ BOOST_CHECK(e1 - e2 == e4);
+ BOOST_CHECK(e1 * e2 == e5);
+ }
+ }
+}
+
+BOOST_AUTO_TEST_CASE(extended_int_test10) {
+ boost::mt19937_64 gen(static_cast<uint32>(time(NULL)));
+ for (int i = 0; i < 100; ++i) {
+ int64 i1 = static_cast<int64>(gen()) >> 20;
+ int64 i2 = i1 >> 32;
+ eint64 e1(i1), e2(i2);
+ BOOST_CHECK(get_d(e1) == static_cast<fpt64>(i1));
+ BOOST_CHECK(get_d(e2) == static_cast<fpt64>(i2));
+ }
+}
+
+BOOST_AUTO_TEST_CASE(extened_int_test11) {
+ eint32 two(2), one(1);
+ eint2048 value(1);
+ for (int i = 0; i < 1024; ++i) {
+ value = value * two;
+ }
+ BOOST_CHECK_EQUAL(value.count(), 33);
+}
+
+BOOST_AUTO_TEST_CASE(robust_sqrt_expr_test1) {
+ robust_sqrt_expr<int32, fpt64> sqrt_expr;
+ int32 A[1] = {10};
+ int32 B[1] = {100};
+ BOOST_CHECK_EQUAL(sqrt_expr.eval1(A, B), 100.0);
+}
+
+BOOST_AUTO_TEST_CASE(robust_sqrt_expr_test2) {
+ robust_sqrt_expr<int32, fpt64> sqrt_expr;
+ int32 A[2] = {10, 30};
+ int32 B[2] = {400, 100};
+ BOOST_CHECK_EQUAL(sqrt_expr.eval2(A, B), 500.0);
+}
+
+BOOST_AUTO_TEST_CASE(robust_sqrt_expr_test3) {
+ robust_sqrt_expr<int32, fpt64> sqrt_expr;
+ int32 A[2] = {10, -30};
+ int32 B[2] = {400, 100};
+ BOOST_CHECK_EQUAL(sqrt_expr.eval2(A, B), -100.0);
+}
+
+BOOST_AUTO_TEST_CASE(robust_sqrt_expr_test4) {
+ robust_sqrt_expr<int32, fpt64> sqrt_expr;
+ int32 A[3] = {10, 30, 20};
+ int32 B[3] = {4, 1, 9};
+ BOOST_CHECK_EQUAL(sqrt_expr.eval3(A, B), 110.0);
+}
+
+BOOST_AUTO_TEST_CASE(robust_sqrt_expr_test5) {
+ robust_sqrt_expr<int32, fpt64> sqrt_expr;
+ int32 A[3] = {10, 30, -20};
+ int32 B[3] = {4, 1, 9};
+ BOOST_CHECK_EQUAL(sqrt_expr.eval3(A, B), -10.0);
+}
+
+BOOST_AUTO_TEST_CASE(robust_sqrt_expr_test6) {
+ robust_sqrt_expr<int32, fpt64> sqrt_expr;
+ int32 A[4] = {10, 30, 20, 5};
+ int32 B[4] = {4, 1, 9, 16};
+ BOOST_CHECK_EQUAL(sqrt_expr.eval4(A, B), 130.0);
+}
+
+BOOST_AUTO_TEST_CASE(robust_sqrt_expr_test7) {
+ robust_sqrt_expr<int32, fpt64> sqrt_expr;
+ int32 A[4] = {10, 30, -20, -5};
+ int32 B[4] = {4, 1, 9, 16};
+ BOOST_CHECK_EQUAL(sqrt_expr.eval4(A, B), -30.0);
+}
+
+BOOST_AUTO_TEST_CASE(robust_sqrt_expr_test8) {
+ robust_sqrt_expr<eint512, efpt64> sqrt_expr;
+ int32 A[4] = {1000, 3000, -2000, -500};
+ int32 B[4] = {400, 100, 900, 1600};
+ eint512 AA[4], BB[4];
+ for (size_t i = 0; i < 4; ++i) {
+ AA[i] = A[i];
+ BB[i] = B[i];
+ }
+ BOOST_CHECK_EQUAL(get_d(sqrt_expr.eval4(AA, BB)), -30000.0);
+}
+
+template <typename _int, typename _fpt>
class sqrt_expr_tester {
public:
+ static const size_t MX_SQRTS = 4;
+
bool run() {
- static boost::mt19937 gen(static_cast<unsigned int>(time(NULL)));
+ static boost::mt19937 gen(static_cast<uint32>(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;
+ for (size_t i = 0; i < MX_SQRTS; ++i) {
+ a[i] = gen() & 1048575;
+ int64 temp = gen() & 1048575;
+ b[i] = temp * temp;
}
- int mask = (1 << MX_SQRTS);
- for (int i = 0; i < mask; i++) {
- double expected_val = 0.0;
+ uint32 mask = (1 << MX_SQRTS);
+ for (size_t i = 0; i < mask; i++) {
+ fpt64 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]);
+ expected_val += static_cast<fpt64>(a[j]) *
+ sqrt(static_cast<fpt64>(b[j]));
} else {
A[j] = -a[j];
B[j] = b[j];
- expected_val -= a[j] * sqrt(b[j]);
+ expected_val -= static_cast<fpt64>(a[j]) *
+ sqrt(static_cast<fpt64>(b[j]));
}
}
- double received_val = get_d(sqrt_expr_.eval4(A, B));
- ret_val &= almost_equal(expected_val, received_val, 1);
+ fpt64 received_val = get_d(sqrt_expr_.eval4(A, B));
+ ret_val &= almost_equal(expected_val, received_val, 25);
}
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];
+ robust_sqrt_expr<_int, _fpt> sqrt_expr_;
+ _int A[MX_SQRTS];
+ _int B[MX_SQRTS];
+ int64 a[MX_SQRTS];
+ int64 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;
+ sqrt_expr_tester<eint1024, efpt64> tester;
for (int i = 0; i < 2000; ++i) {
- BOOST_CHECK(tester1.run());
- BOOST_CHECK(tester2.run());
+ BOOST_CHECK(tester.run());
}
}
Modified: sandbox/gtl/libs/polygon/test/voronoi_test_helper.hpp
==============================================================================
--- sandbox/gtl/libs/polygon/test/voronoi_test_helper.hpp (original)
+++ sandbox/gtl/libs/polygon/test/voronoi_test_helper.hpp 2011-12-24 12:20:30 EST (Sat, 24 Dec 2011)
@@ -20,7 +20,7 @@
enum kOrientation {
RIGHT = -1,
COLLINEAR = 0,
- LEFT = 1,
+ LEFT = 1
};
template <typename Point2D>
@@ -204,8 +204,10 @@
void save_voronoi_input(const std::vector< point_data<T> > &points, const char *file_name) {
std::ofstream ofs(file_name);
ofs << points.size() << std::endl;
- for (int i = 0; i < (int)points.size(); ++i)
- ofs << points[i].x() << " " << points[i].y() << std::endl;
+ for (typename std::vector< point_data<T> >::iterator it = points.begin();
+ it != points.end(); ++it) {
+ ofs << it->x() << " " << it->y() << std::endl;
+ }
ofs.close();
}
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