diff --git a/extensions/example/Jamfile b/extensions/example/Jamfile index 8110621689..aca687d230 100644 --- a/extensions/example/Jamfile +++ b/extensions/example/Jamfile @@ -13,3 +13,4 @@ project boost-geometry-examples-extensions ; build-project gis ; +build-project generic_robust_predicates ; diff --git a/extensions/example/generic_robust_predicates/Jamfile b/extensions/example/generic_robust_predicates/Jamfile new file mode 100644 index 0000000000..6910d93912 --- /dev/null +++ b/extensions/example/generic_robust_predicates/Jamfile @@ -0,0 +1,12 @@ +# Boost.Geometry (aka GGL, Generic Geometry Library) + +# Use, modification and distribution is subject to 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) + + +project boost-geometry-example-extensions-generic_robust_predicates + : + ; + +exe static_side_2d : static_side_2d.cpp ; diff --git a/extensions/example/generic_robust_predicates/static_side_2d.cpp b/extensions/example/generic_robust_predicates/static_side_2d.cpp new file mode 100644 index 0000000000..04ba5b0515 --- /dev/null +++ b/extensions/example/generic_robust_predicates/static_side_2d.cpp @@ -0,0 +1,79 @@ +// Boost.Geometry (aka GGL, Generic Geometry Library) + +// Copyright (c) 2020 Tinko Bartels, Berlin, Germany. + +// Contributed and/or modified by Tinko Bartels, +// as part of Google Summer of Code 2020 program. + +// Use, modification and distribution is subject to 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) + + +#define BOOST_GEOMETRY_NO_BOOST_TEST + +#include + +#include +#include + +#include +#include +#include + + +namespace bg = boost::geometry; +using point = bg::model::point; + +template +struct side_robust_with_static_filter +{ +private: + using ct = CalculationType; + using expression = bg::detail::generic_robust_predicates::orient2d; + using filter = bg::detail::generic_robust_predicates::stage_a_static + < + expression, + ct + >; + filter m_filter; +public: + side_robust_with_static_filter(ct x_max, ct y_max, ct x_min, ct y_min) + : m_filter(x_max, y_max, x_max, y_max, x_max, y_max, + x_min, y_min, x_min, y_min, x_min, y_min) {}; + + template + < + typename P1, + typename P2, + typename P + > + inline int apply(P1 const& p1, P2 const& p2, P const& p) const + { + int sign = m_filter.apply(bg::get<0>(p1), + bg::get<1>(p1), + bg::get<0>(p2), + bg::get<1>(p2), + bg::get<0>(p), + bg::get<1>(p)); + if (sign != bg::detail::generic_robust_predicates::sign_uncertain) + { + return sign; + } + else + { + //fallback if filter fails. + return bg::strategy::side::side_robust::apply(p1, p2, p); + } + } +}; + +int main() +{ + point p1(0.0, 0.0); + point p2(1.0, 1.0); + point p (0.0, 1.0); + side_robust_with_static_filter static_strategy(2.0, 2.0, 1.0, 1.0); + std::cout << "Side value: " << static_strategy.apply(p1, p2, p) << "\n"; + return 0; +} diff --git a/extensions/test/Jamfile b/extensions/test/Jamfile index 7b0af8bed3..161443eb02 100644 --- a/extensions/test/Jamfile +++ b/extensions/test/Jamfile @@ -28,3 +28,4 @@ build-project gis ; build-project iterators ; build-project nsphere ; build-project triangulation ; +build-project generic_robust_predicates ; diff --git a/extensions/test/generic_robust_predicates/Jamfile b/extensions/test/generic_robust_predicates/Jamfile new file mode 100644 index 0000000000..ef072ff343 --- /dev/null +++ b/extensions/test/generic_robust_predicates/Jamfile @@ -0,0 +1,13 @@ +# Boost.Geometry (aka GGL, Generic Geometry Library) +# +# Use, modification and distribution is subject to 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) + +test-suite boost-geometry-extensions-generic_robust_predicates + : + [ run expression_eval.cpp ] + [ run side3d.cpp ] + [ run staged.cpp ] + ; + diff --git a/extensions/test/generic_robust_predicates/expression_eval.cpp b/extensions/test/generic_robust_predicates/expression_eval.cpp new file mode 100644 index 0000000000..e6a2c5972f --- /dev/null +++ b/extensions/test/generic_robust_predicates/expression_eval.cpp @@ -0,0 +1,62 @@ +// Boost.Geometry (aka GGL, Generic Geometry Library) +// Unit Test + +// Copyright (c) 2020 Tinko Bartels, Berlin, Germany. + +// Contributed and/or modified by Tinko Bartels, +// as part of Google Summer of Code 2020 program. + +// Use, modification and distribution is subject to 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) + +#include + +#include + +#include +#include + +template +void test_all() +{ + using bg::detail::generic_robust_predicates::_1; + using bg::detail::generic_robust_predicates::_2; + using bg::detail::generic_robust_predicates::_3; + using bg::detail::generic_robust_predicates::_4; + using bg::detail::generic_robust_predicates::sum; + using bg::detail::generic_robust_predicates::difference; + using bg::detail::generic_robust_predicates::product; + using bg::detail::generic_robust_predicates::max; + using bg::detail::generic_robust_predicates::abs; + using bg::detail::generic_robust_predicates::post_order; + using bg::detail::generic_robust_predicates::evaluate_expression; + using bg::detail::generic_robust_predicates::evaluate_expressions; + using ct = CalculationType; + ct r1 = evaluate_expression>( + std::array{1.0, 2.0}); + BOOST_CHECK_EQUAL(3.0, r1); + ct r2 = evaluate_expression, abs<_2>>>( + std::array{-10.0, 2.0}); + BOOST_CHECK_EQUAL(10.0, r2); + + using expression = product + < + difference<_1, _2>, + difference<_3, _4> + >; + using evals = post_order; + std::array::value> r; + std::array input {5.0, 3.0, 2.0, 8.0}; + evaluate_expressions(input, r, evals{}); + BOOST_CHECK_EQUAL(2.0, r[0]); + BOOST_CHECK_EQUAL(-6.0, r[1]); + BOOST_CHECK_EQUAL(-12.0, r[2]); +} + + +int test_main(int, char* []) +{ + test_all(); + return 0; +} diff --git a/extensions/test/generic_robust_predicates/side3d.cpp b/extensions/test/generic_robust_predicates/side3d.cpp new file mode 100644 index 0000000000..4b9ac945d7 --- /dev/null +++ b/extensions/test/generic_robust_predicates/side3d.cpp @@ -0,0 +1,76 @@ +// Boost.Geometry (aka GGL, Generic Geometry Library) +// Unit Test + +// Copyright (c) 2020 Tinko Bartels, Berlin, Germany. + +// Contributed and/or modified by Tinko Bartels, +// as part of Google Summer of Code 2020 program. + +// Use, modification and distribution is subject to 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) + +#include + +#include + +#include +#include + +template +void test_all() +{ + using bg::detail::generic_robust_predicates::orient3d; + using bg::detail::generic_robust_predicates::stage_a_semi_static; + using bg::detail::generic_robust_predicates::stage_a_almost_static; + using bg::detail::generic_robust_predicates::stage_a_static; + using ct = CalculationType; + using semi_static = stage_a_semi_static; + BOOST_CHECK_EQUAL(1, + semi_static::apply(1, 0, 0, + 0, 1, 0, + 1, 1, 0, + 0, 0, 1)); + BOOST_CHECK_EQUAL(-1, + semi_static::apply(1, 0, 0, + 0, 1, 0, + 1, 1, 0, + 0, 0, -1)); + BOOST_CHECK_EQUAL(0, + semi_static::apply(1, 0, 0, + 0, 1, 0, + 1, 1, 0, + 0, 0, 0)); + stage_a_static stat(10, 10, 10, + 10, 10, 10, + 10, 10, 10, + 10, 10, 10, + 0, 0, 0, + 0, 0, 0, + 0, 0, 0, + 0, 0, 0); + BOOST_CHECK_EQUAL(1, stat.apply(1, 0, 0, + 0, 1, 0, + 1, 1, 0, + 0, 0, 1)); + + stage_a_static pessimistic(1e40, 1e40, 1e40, + 1e40, 1e40, 1e40, + 1e40, 1e40, 1e40, + 1e40, 1e40, 1e40, + 0, 0, 0, + 0, 0, 0, + 0, 0, 0, + 0, 0, 0); + BOOST_CHECK_EQUAL(bg::detail::generic_robust_predicates::sign_uncertain, + pessimistic.apply(1, 0, 0, + 0, 1, 0, + 1, 1, 0, + 0, 0, 1)); +} + +int test_main(int, char* []) +{ + test_all(); + return 0; +} diff --git a/extensions/test/generic_robust_predicates/staged.cpp b/extensions/test/generic_robust_predicates/staged.cpp new file mode 100644 index 0000000000..5198cec251 --- /dev/null +++ b/extensions/test/generic_robust_predicates/staged.cpp @@ -0,0 +1,43 @@ +// Boost.Geometry (aka GGL, Generic Geometry Library) +// Unit Test + +// Copyright (c) 2020 Tinko Bartels, Berlin, Germany. + +// Contributed and/or modified by Tinko Bartels, +// as part of Google Summer of Code 2020 program. + +// Use, modification and distribution is subject to 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) + +#include + +#include +#include +#include +#include +#include + +using namespace boost::geometry::detail::generic_robust_predicates; + +template +void test_all() +{ + using ct = CalculationType; + using expression = orient2d; + using filter1 = stage_a_static; + using filter2 = stage_a_almost_static; + using filter3 = stage_a_semi_static; + using filter4 = stage_b; + using filter5 = stage_d; + using staged = staged_predicate; + staged s(1e20, 1e20, 1e20, 1e20, 1e20, 1e20, 0., 0., 0., 0., 0., 0.); + s.update(1e20, 1e20, 1e20, 1e20, 1e20, 1e20, 0., 0., 0., 0., 0., 0.); + BOOST_CHECK_EQUAL(1, s.apply(1e-20, 1e-20, 1, 1-1e-10, 1e20, 1e20)); +} + +int test_main(int, char* []) +{ + test_all(); + return 0; +} diff --git a/include/boost/geometry/extensions/generic_robust_predicates/strategies/cartesian/detail/almost_static_filter.hpp b/include/boost/geometry/extensions/generic_robust_predicates/strategies/cartesian/detail/almost_static_filter.hpp new file mode 100644 index 0000000000..6c2e9bdd92 --- /dev/null +++ b/include/boost/geometry/extensions/generic_robust_predicates/strategies/cartesian/detail/almost_static_filter.hpp @@ -0,0 +1,129 @@ +// Boost.Geometry (aka GGL, Generic Geometry Library) + +// Copyright (c) 2020 Tinko Bartels, Berlin, Germany. + +// Contributed and/or modified by Tinko Bartels, +// as part of Google Summer of Code 2020 program. + +// Use, modification and distribution is subject to 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) + +#ifndef BOOST_GEOMETRY_EXTENSIONS_GENERIC_ROBUST_PREDICATES_STRATEGIES_CARTESIAN_DETAIL_ALMOST_STATIC_FILTER_HPP +#define BOOST_GEOMETRY_EXTENSIONS_GENERIC_ROBUST_PREDICATES_STRATEGIES_CARTESIAN_DETAIL_ALMOST_STATIC_FILTER_HPP + +#include +#include +#include + +#include + +namespace boost { namespace geometry +{ + +namespace detail { namespace generic_robust_predicates +{ + +//The almost static filter holds an instance of a static filter and applies it +//when it is called. Its static filter can be updated and applied like this: +// +//almost_static_filter<...> f; +// +//f.update(max1, max2, ..., min1, min2, ...); +// +//f.apply(arg1, arg2, ...); +// +//Unlike with the static filter, global bounds do not need to be known at +//construction time and with incremental algorithms where inputs with higher +//magnitude are added later, the earlier applications of the filter may benefit +//from more optimistic error bounds. + +template +< + typename Expression, + typename CalculationType, + typename StaticFilter +> +class almost_static_filter +{ +private: + using ct = CalculationType; + using extrema_array = std::array>; + extrema_array m_extrema; + StaticFilter m_filter; +public: + static constexpr bool stateful = true; + static constexpr bool updates = true; + + StaticFilter const& filter() const { return m_filter; } + inline almost_static_filter() + { + std::fill(m_extrema.begin(), + m_extrema.begin() + m_extrema.size() / 2, + -std::numeric_limits::infinity()); + std::fill(m_extrema.begin() + m_extrema.size() / 2, + m_extrema.end(), + std::numeric_limits::infinity()); + } + template + int apply(CTs const&... args) const + { + return m_filter.apply(args...); + } + + template + inline void update_extrema(CTs const&... args) + { + std::array input {{ static_cast(args)... }}; + for(unsigned int i = 0; i < m_extrema.size() / 2; ++i) + { + m_extrema[i] = std::max(m_extrema[i], input[i]); + } + for(unsigned int i = m_extrema.size() / 2; i < m_extrema.size(); ++i) + { + m_extrema[i] = std::min(m_extrema[i], input[i]); + } + } + + template + inline bool update_extrema_check(CTs const&... args) + { + bool changed = false; + std::array input {{ static_cast(args)... }}; + for(unsigned int i = 0; i < m_extrema.size() / 2; ++i) + { + if (input[i] > m_extrema[i]) + { + changed = true; + m_extrema[i] = input[i]; + } + } + for(unsigned int i = m_extrema.size() / 2; i < m_extrema.size(); ++i) + { + if (input[i] < m_extrema[i]) + { + changed = true; + m_extrema[i] = input[i]; + } + } + return changed; + } + + template + inline void update(CTs const&... args) + { + update_extrema(args...); + update_filter(); + } + + inline void update_filter() + { + m_filter = StaticFilter(m_extrema); + } +}; + +}} // namespace detail::generic_robust_predicates + +}} // namespace boost::geometry + +#endif // BOOST_GEOMETRY_EXTENSIONS_GENERIC_ROBUST_PREDICATES_STRATEGIES_CARTESIAN_DETAIL_ALMOST_STATIC_FILTER_HPP diff --git a/include/boost/geometry/extensions/generic_robust_predicates/strategies/cartesian/detail/expansion_arithmetic.hpp b/include/boost/geometry/extensions/generic_robust_predicates/strategies/cartesian/detail/expansion_arithmetic.hpp new file mode 100644 index 0000000000..e1a76fdcb9 --- /dev/null +++ b/include/boost/geometry/extensions/generic_robust_predicates/strategies/cartesian/detail/expansion_arithmetic.hpp @@ -0,0 +1,2125 @@ +// Boost.Geometry (aka GGL, Generic Geometry Library) + +// Copyright (c) 2020 Tinko Bartels, Berlin, Germany. + +// Contributed and/or modified by Tinko Bartels, +// as part of Google Summer of Code 2020 program. + +// Use, modification and distribution is subject to 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) + +#ifndef BOOST_GEOMETRY_EXTENSIONS_GENERIC_ROBUST_PREDICATES_STRATEGIES_CARTESIAN_DETAIL_EXPANSION_ARITHMETIC_HPP +#define BOOST_GEOMETRY_EXTENSIONS_GENERIC_ROBUST_PREDICATES_STRATEGIES_CARTESIAN_DETAIL_EXPANSION_ARITHMETIC_HPP + +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include + +namespace boost { namespace geometry +{ + +namespace detail { namespace generic_robust_predicates +{ + +struct abs_comp +{ + template + constexpr bool operator()(Real a, Real b) const + { + return std::abs(a) < std::abs(b); + } +}; + +template +struct negate_impl +{ + static constexpr Real apply(Real a) + { + return a; + } +}; + +template +struct negate_impl +{ + static constexpr Real apply(Real a) + { + return -a; + } +}; + +template +constexpr Real negate(Real a) +{ + return negate_impl::apply(a); +} + +namespace debug_expansion +{ + +constexpr unsigned long long round_to_power_of_two(unsigned long long num) +{ + if (!num) + { + return 0; + } + unsigned long long ret = 1; + while (num >>= 1) + { + ret <<= 1; + } + return ret; +} + +template +inline bool nonoverlapping(Real a, Real b) +{ + int a_exp, b_exp, min_exp, max_exp; + Real a_mant = std::frexp(a, &a_exp); + Real b_mant = std::frexp(b, &b_exp); + unsigned long long scale = 1ULL<<63; + unsigned long long a_mantll = std::abs(a_mant) * scale; + unsigned long long b_mantll = std::abs(b_mant) * scale; + if (a_mantll == 0 || b_mantll == 0) + { + return true; + } + unsigned long long min_mantll, max_mantll; + if (a_exp < b_exp) + { + min_exp = a_exp; + max_exp = b_exp; + min_mantll = a_mantll; + max_mantll = b_mantll; + } + else + { + min_exp = b_exp; + max_exp = a_exp; + min_mantll = b_mantll; + max_mantll = a_mantll; + } + int scale_down = max_exp - min_exp; + if (scale_down > std::numeric_limits::digits) + { + return true; + } + unsigned long long min_mantll_sc = min_mantll >> scale_down; + auto min_mantll_sc_rd = round_to_power_of_two(min_mantll_sc); + bool result = (max_mantll % ( 2 * min_mantll_sc_rd )) == 0; + + return result; +} + +template +inline bool nonadjacent(Real a, Real b) +{ + bool t1 = nonoverlapping(a, b); + bool t2 = nonoverlapping(a, 2*b); + bool t3 = nonoverlapping(2*a, b); + return t1 && t2 && t3; +} + +template +inline bool expansion_nonoverlapping(Iter begin, Iter end) +{ + auto lesser = *begin; + for(auto it = begin + 1; it < end; ++it) + { + if ( *it != 0) + { + if (lesser > std::abs(*it) || !nonoverlapping(lesser, *it) ) + { + return false; + } + lesser = *it; + } + } + return true; +} + +template +inline bool expansion_nonadjacent(Iter begin, Iter end) +{ + auto lesser = *begin; + for(auto it = begin + 1; it < end; ++it) + { + if ( *it != 0) + { + if (lesser > std::abs(*it) || !nonadjacent(lesser, *it)) + { + return false; + } + lesser = *it; + } + } + return true; +} + +template +inline bool expansion_strongly_nonoverlapping(Iter begin, Iter end) +{ + using Real = typename std::iterator_traits::value_type; + Real lesser = *begin; + Real previous = 0.0; + for(auto it = begin + 1; it < end; ++it) + { + if ( *it != 0) + { + if (lesser > std::abs(*it) || !nonoverlapping(lesser, *it)) + { + return false; + } + if ( !nonadjacent(lesser, *it) ) + { + int exp_now, exp_lesser; + std::frexp(lesser, &exp_lesser); + std::frexp(*it, &exp_now); + if ( std::abs(std::frexp(lesser, &exp_lesser)) != 0.5 + || std::abs(std::frexp(*it, &exp_now)) != 0.5 ) + { + return false; + } + if ( !nonadjacent( lesser, previous ) ) + { + return false; + } + } + previous = lesser; + lesser = *it; + } + } + return true; +} + +} // namespace debug_expansion + +template +constexpr Real two_sum_tail(Real a, Real b, Real x) +{ + Real b_virtual = x - a; + Real a_virtual = x - b_virtual; + Real b_rounded = b - b_virtual; + Real a_rounded = a - a_virtual; + Real y = a_rounded + b_rounded; + return y; +} + +template +constexpr Real fast_two_sum_tail(Real a, Real b, Real x) +{ + assert(std::abs(a) >= std::abs(b) || a == 0 || b == 0); + Real b_virtual = x - a; + Real y = b - b_virtual; + return y; +} + +template +constexpr Real two_difference_tail(Real a, Real b, Real x) +{ + Real b_virtual = a - x; + Real a_virtual = x + b_virtual; + Real b_rounded = b_virtual - b; + Real a_rounded = a - a_virtual; + Real y = a_rounded + b_rounded; + return y; +} + +template +constexpr Real fast_two_difference_tail(Real a, Real b, Real x) +{ + assert(std::abs(a) >= std::abs(b) || a == 0); + Real b_virtual = a - x; + Real y = b_virtual - b; + return y; +} + +template +constexpr Real two_product_tail(Real a, Real b, Real x) +{ + Real y = std::fma(a, b, -x); + return y; +} + +template +constexpr Real splitter_helper() +{ + int digits = std::numeric_limits::digits; + int half = digits / 2; + int ceilhalf = half + (half * 2 == digits ? 0 : 1 ); + Real out(1); + for(int i = 0; i < ceilhalf; ++i) + { + out *= Real(2); + } + return out + Real(1); +} + +template +constexpr Real splitter = splitter_helper(); + +template +constexpr std::array split(Real a) +{ + Real c = splitter * a; + Real a_big = c - a; + Real a_hi = c - a_big; + Real a_lo = a - a_hi; + return std::array{ a_hi, a_lo }; +} + +template +constexpr Real two_product_tail_constexpr(Real a, Real b, Real x) +{ + const auto a_split = split(a); + const auto b_split = split(b); + Real a_hi_b_hi = a_split[0] * b_split[0]; + Real err1 = x - a_hi_b_hi; + Real a_lo_b_hi = a_split[1] * b_split[0]; + Real err2 = err1 - a_lo_b_hi; + Real a_hi_b_lo = a_split[0] * b_split[1]; + Real err3 = err2 - a_hi_b_lo; + Real a_lo_b_lo = a_split[1] * b_split[1]; + return a_lo_b_lo - err3; +} + +template +< + bool ZeroElimination, + bool MostSigOnly, + typename OutIter, + typename Real +> +struct insert_ze_impl +{ + static constexpr OutIter apply(OutIter out, Real val) + { + if (val == Real(0)) + { + return out; + } + else + { + *out = val; + ++out; + return out; + } + } +}; + +template +< + typename OutIter, + typename Real +> +struct insert_ze_impl +{ + static constexpr OutIter apply(OutIter out, Real val) + { + *out += val; + return out; + } +}; + +template +< + typename OutIter, + typename Real +> +struct insert_ze_impl +{ + static constexpr OutIter apply(OutIter out, Real val) + { + if (val != 0) + { + *out = val; + } + return out; + } +}; + +template +< + typename OutIter, + typename Real +> +struct insert_ze_impl +{ + static constexpr OutIter apply(OutIter out, Real val) + { + *out = val; + ++out; + return out; + } +}; + +template +< + bool ZeroElimination, + bool MostSigOnly = false, + typename OutIter, + typename Real +> +constexpr OutIter insert_ze(OutIter o, Real r) +{ + return insert_ze_impl + < + ZeroElimination, + MostSigOnly, + OutIter, + Real + >::apply(o, r); +} + +template +< + bool ZeroElimination, + bool MostSigOnly, + typename OutIter, + typename Real +> +struct insert_ze_final_impl +{ + static constexpr OutIter apply(OutIter out, OutIter start, Real val) + { + if (val == Real(0) && out != start) + { + return out; + } + else + { + *out = val; + ++out; + return out; + } + } +}; + +template +< + typename OutIter, + typename Real +> +struct insert_ze_final_impl +{ + static constexpr OutIter apply(OutIter out, OutIter, Real val) + { + *out += val; + ++out; + return out; + } +}; + +template +< + typename OutIter, + typename Real +> +struct insert_ze_final_impl +{ + static constexpr OutIter apply(OutIter out, OutIter, Real val) + { + if (val != 0) + { + *out = val; + } + ++out; + return out; + } +}; + +template +< + typename OutIter, + typename Real +> +struct insert_ze_final_impl +{ + static constexpr OutIter apply(OutIter out, OutIter, Real val) + { + *out = val; + ++out; + return out; + } +}; + +template +< + bool ZeroElimination, + bool MostSigOnly = false, + typename OutIter, + typename Real +> +constexpr OutIter insert_ze_final(OutIter o, OutIter s, Real r) +{ + return insert_ze_final_impl + < + ZeroElimination, + MostSigOnly, + OutIter, + Real + >::apply(o, s, r); +} + +template +< + bool ZeroElimination, + bool MostSigOnly = false, + typename InIter, + typename OutIter, + typename Real, + bool NegateE = false, + bool NegateB = false +> +constexpr OutIter grow_expansion(InIter e_begin, + InIter e_end, + Real b, + OutIter h_begin, + OutIter) +{ + assert(debug_expansion::expansion_nonoverlapping(e_begin, e_end)); + Real Q = negate(b); + auto h_it = h_begin; + for(auto e_it = e_begin; e_it != e_end; ++e_it) + { + Real Q_new = negate(*e_it) + Q; + Real h_new = two_sum_tail(negate(*e_it), Q, Q_new); + Q = Q_new; + h_it = insert_ze(h_it, h_new); + } + h_it = insert_ze_final(h_it, h_begin, Q); + assert(debug_expansion::expansion_nonoverlapping(h_begin, h_it)); + assert( !debug_expansion::expansion_nonadjacent(e_begin, e_end) + || debug_expansion::expansion_nonadjacent(h_begin, h_it) ); + return h_it; +} + +template +struct expansion_sum_advance_impl +{ + template + static constexpr bool apply(Real, Real) + { + return true; + } +}; + +template +struct expansion_sum_advance_impl +{ + template + static constexpr bool apply(Real e, Real b) + { + Real Q = negate(e) + negate(b); + return two_sum_tail(negate(e), + negate(b), Q) + != Real(0); + } +}; + + +template +constexpr bool expansion_sum_advance(Real e, Real b) +{ + return expansion_sum_advance_impl + ::apply(e, b); +} + +template +< + bool ZeroElimination, + typename InIter1, + typename InIter2, + typename OutIter, + bool NegateE = false, + bool NegateF = false +> +constexpr OutIter expansion_sum(InIter1 e_begin, + InIter1 e_end, + InIter2 f_begin, + InIter2 f_end, + OutIter h_begin, + OutIter) +{ + using Real = typename std::iterator_traits::value_type; + assert(debug_expansion::expansion_nonoverlapping(e_begin, e_end)); + assert(debug_expansion::expansion_nonoverlapping(f_begin, f_end)); + const auto elen = std::distance(e_begin, e_end); + const auto flen = std::distance(f_begin, f_end); + auto f_it = f_begin; + auto h_begin_i = h_begin; + bool advance = + expansion_sum_advance(*e_begin, + *f_it); + auto h_it = grow_expansion + < + ZeroElimination, + false, + InIter1, + OutIter, + Real, + NegateE, + NegateF + >(e_begin, e_end, *f_it, h_begin, h_begin + elen + 1); + if (advance) + { + ++h_begin_i; + } + ++f_it; + for(auto i = 1; i < flen; ++i) + { + advance = + expansion_sum_advance(*h_begin_i, + *f_it); + h_it = grow_expansion + < + ZeroElimination, + false, + InIter1, + OutIter, + Real, + false, + NegateF + >(h_begin_i, + h_it, + *f_it, + h_begin_i, + h_it + 1); + if (advance) + { + ++h_begin_i; + } + ++f_it; + } + assert(debug_expansion::expansion_nonoverlapping(h_begin, h_it)); + assert( !debug_expansion::expansion_nonadjacent(e_begin, e_end) + || !debug_expansion::expansion_nonadjacent(f_begin, f_end) + || debug_expansion::expansion_nonadjacent(h_begin, h_it)); + return h_it; +} + +template +< + bool ZeroElimination, + typename InIter1, + typename InIter2, + typename OutIter, + bool NegateE = false, + bool NegateF = false +> +constexpr OutIter expansion_sum_most_sig(InIter1 e_begin, + InIter1 e_end, + InIter2 f_begin, + InIter2 f_end, + OutIter h_begin, + OutIter) +//This procedure destroys the contents of [e_begin, e_end) +{ + using real = typename std::iterator_traits::value_type; + for(auto f_it = f_begin; f_it != f_end - 1; ++f_it) + { + const auto e = *e_begin; + const auto f = negate(*f_it); + const auto Q = e + f; + const auto h = two_sum_tail(e, f, Q); + insert_ze(h_begin, h); + e_end = grow_expansion + < + ZeroElimination, + false, + InIter1, + InIter1, + real, + false, + false + >(e_begin + 1, + e_end, + Q, + e_begin, + e_end); + } + auto h_it = grow_expansion + < + ZeroElimination, + true, + InIter1, + OutIter, + real, + false, + NegateE != NegateF + >(e_begin, + e_end, + *(f_end - 1), + h_begin, + h_begin + 1); + *h_begin = negate(*h_begin); + return h_it; +} + +template +< + bool ZeroElimination, + bool MostSigOnly = false, + typename InIter1, + typename InIter2, + typename OutIter, + bool NegateE = false, + bool NegateF = false, + bool ENoZeros = false, + bool FNoZeros = false +> +constexpr OutIter fast_expansion_sum_inplace(InIter1 e_begin, + InIter1 e_end, + InIter2 f_begin, + InIter2 f_end, + OutIter h_begin, + OutIter h_end) +{ +#ifndef NDEBUG + auto e_approx = std::accumulate(e_begin, e_end, 0.); + auto f_approx = std::accumulate(f_begin, f_end, 0.); +#endif + + assert(e_end <= f_begin); + assert(f_begin != h_begin); + using Real = typename std::iterator_traits::value_type; + assert(debug_expansion::expansion_nonoverlapping(e_begin, e_end)); + assert(debug_expansion::expansion_nonoverlapping(f_begin, f_end)); + if (NegateE) + { + for(auto e_it = e_begin; e_it != e_end; ++e_it) + { + *e_it = -*e_it; + } + } + if (NegateF) + { + for(auto f_it = f_begin; f_it != f_end; ++f_it) + { + *f_it = -*f_it; + } + } + auto e_old_end = e_end; + typename std::iterator_traits::difference_type e_shortened = 0; + if (!ENoZeros) + { + e_end = std::remove(e_begin, e_end, Real(0)); + e_shortened = std::distance(e_end, e_old_end); + } + + std::rotate(e_end, f_begin, f_end); + auto f_old_end = f_end; + f_end = f_end - e_shortened - std::distance(e_old_end, f_begin); + f_begin = e_end; + + if (!FNoZeros) + { + f_end = std::remove(f_begin, f_end, Real(0)); + } + if (!ZeroElimination) + { + std::fill(f_end, f_old_end, Real(0)); + } + + std::inplace_merge(e_begin, e_end, f_end, abs_comp{}); + + if (std::distance(e_begin, f_end) == 1) + { + *h_begin = *e_begin; + if ( !ZeroElimination ) + { + std::fill(h_begin + 1, h_end, 0.); + return h_end; + } + return h_begin + 1; + } + + InIter1 g_it = e_begin; + InIter2 g_end = f_end; + OutIter h_it = h_begin; + Real Q = *g_it + *(g_it + 1); + Real h_new = fast_two_sum_tail(*(g_it + 1), *(g_it), Q); + h_it = insert_ze(h_it, h_new); + g_it += 2; + for(; g_it < g_end; ++g_it) + { + Real Q_new = Q + *g_it; + h_new = two_sum_tail(Q, *g_it, Q_new); + h_it = insert_ze(h_it, h_new); + Q = Q_new; + } + h_it = insert_ze_final(h_it, h_begin, Q); + + assert(debug_expansion::expansion_nonoverlapping(h_begin, h_it)); + if ( !ZeroElimination ) + { + h_it = h_end; + } + +#ifndef NDEBUG + auto h_approx = std::accumulate(h_begin, h_it, 0.); + assert( + MostSigOnly + || std::abs( (NegateE ? -1 : 1) * e_approx + + (NegateF ? -1 : 1) * f_approx - h_approx ) + <= 1e-5 * ( std::abs(e_approx) + std::abs(f_approx)) + ); +#endif + return h_it; +} + +template +< + bool ZeroElimination, + bool MostSigOnly = false, + typename InIter1, + typename InIter2, + typename OutIter, + bool NegateE = false, + bool NegateF = false +> +constexpr OutIter fast_expansion_sum_not_inplace(InIter1 e_begin, + InIter1 e_end, + InIter2 f_begin, + InIter2 f_end, + OutIter h_begin, + OutIter h_end) +{ + assert(e_begin != h_begin); + assert(f_begin != h_begin); + using Real = typename std::iterator_traits::value_type; + assert(debug_expansion::expansion_nonoverlapping(e_begin, e_end)); + assert(debug_expansion::expansion_nonoverlapping(f_begin, f_end)); + + auto e_it = e_begin; + auto f_it = f_begin; + Real Q; + if (std::abs(*f_it) > std::abs(*e_it)) + { + Q = negate(*e_it); + ++e_it; + } + else + { + Q = negate(*f_it); + ++f_it; + } + auto h_it = h_begin; + if ((e_it != e_end) && (f_it != f_end)) + { + Real Q_new; + Real h_new; + if (std::abs(*f_it) > std::abs(*e_it)) + { + Q_new = negate(*e_it) + Q; + h_new = fast_two_sum_tail(negate(*e_it), Q, Q_new); + ++e_it; + } + else + { + Q_new = negate(*f_it) + Q; + h_new = fast_two_sum_tail(negate(*f_it), Q, Q_new); + ++f_it; + } + Q = Q_new; + h_it = insert_ze(h_it, h_new); + while((e_it != e_end) && (f_it != f_end)) + { + if (std::abs(*f_it) > std::abs(*e_it)) + { + Q_new = negate(*e_it) + Q; + h_new = two_sum_tail(negate(*e_it), Q, Q_new); + ++e_it; + } + else + { + Q_new = negate(*f_it) + Q; + h_new = two_sum_tail(negate(*f_it), Q, Q_new); + ++f_it; + } + Q = Q_new; + h_it = insert_ze(h_it, h_new); + } + } + while(e_it != e_end) + { + Real Q_new = negate(*e_it) + Q; + Real h_new = two_sum_tail(negate(*e_it), Q, Q_new); + h_it = insert_ze(h_it, h_new); + Q = Q_new; + ++e_it; + } + while(f_it != f_end) + { + Real Q_new = negate(*f_it) + Q; + Real h_new = two_sum_tail(negate(*f_it), Q, Q_new); + h_it = insert_ze(h_it, h_new); + Q = Q_new; + ++f_it; + } + h_it = insert_ze_final(h_it, h_begin, Q); + assert(debug_expansion::expansion_nonoverlapping(h_begin, h_it)); + return h_it; +} + +template +< + bool ZeroElimination, + bool MostSigOnly, + typename InIter1, + typename InIter2, + typename OutIter, + bool Inplace, + bool NegateE, + bool NegateF +> +struct fast_expansion_sum_impl +{ + static constexpr OutIter apply(InIter1 e_begin, + InIter1 e_end, + InIter2 f_begin, + InIter2 f_end, + OutIter h_begin, + OutIter h_end) + { + return fast_expansion_sum_not_inplace + < + ZeroElimination, + MostSigOnly, + InIter1, + InIter2, + OutIter, + NegateE, + NegateF + >(e_begin, e_end, f_begin, f_end, h_begin, h_end); + } +}; + +template +< + bool ZeroElimination, + bool MostSigOnly, + typename InIter1, + typename InIter2, + typename OutIter, + bool NegateE, + bool NegateF +> +struct fast_expansion_sum_impl + < + ZeroElimination, + MostSigOnly, + InIter1, + InIter2, + OutIter, + true, + NegateE, + NegateF + > +{ + static constexpr OutIter apply(InIter1 e_begin, + InIter1 e_end, + InIter2 f_begin, + InIter2 f_end, + OutIter h_begin, + OutIter h_end) + { + return fast_expansion_sum_inplace + < + ZeroElimination, + MostSigOnly, + InIter1, + InIter2, + OutIter, + NegateE, + NegateF + >(e_begin, e_end, f_begin, f_end, h_begin, h_end); + } +}; + +template +< + bool ZeroElimination, + bool MostSigOnly = false, + typename InIter1, + typename InIter2, + typename OutIter, + bool Inplace, + bool NegateE, + bool NegateF +> +constexpr OutIter fast_expansion_sum(InIter1 e_begin, + InIter1 e_end, + InIter2 f_begin, + InIter2 f_end, + OutIter h_begin, + OutIter h_end) +{ + return fast_expansion_sum_impl + < + ZeroElimination, + MostSigOnly, + InIter1, + InIter2, + OutIter, + Inplace, + NegateE, + NegateF + >::apply(e_begin, e_end, f_begin, f_end, h_begin, h_end); +} + +template +< + bool ZeroElimination, + typename InIter, + typename Real, + typename OutIter +> +constexpr OutIter scale_expansion(InIter e_begin, + InIter e_end, + Real b, + OutIter h_begin, + OutIter) +{ + assert(debug_expansion::expansion_nonoverlapping(e_begin, e_end)); + + auto e_it = e_begin; + auto h_it = h_begin; + Real Q = *e_it * b; + auto h_new = two_product_tail(*e_it, b, Q); + h_it = insert_ze(h_it, h_new); + ++e_it; + for(; e_it != e_end; ++e_it) + { + Real product_1 = *e_it * b; + Real product_0 = two_product_tail(*e_it, b, product_1); + Real sum = Q + product_0; + h_new = two_sum_tail(Q, product_0, sum); + h_it = insert_ze(h_it, h_new); + Q = product_1 + sum; + h_new = two_sum_tail(product_1, sum, Q); + h_it = insert_ze(h_it, h_new); + } + h_it = insert_ze_final(h_it, h_begin, Q); + + assert( debug_expansion::expansion_nonoverlapping(h_begin, h_it) ); + assert( !debug_expansion::expansion_nonadjacent(e_begin, e_end) + || debug_expansion::expansion_nonadjacent(h_begin, h_it) ); + assert( !debug_expansion::expansion_strongly_nonoverlapping(e_begin, e_end) + || debug_expansion::expansion_strongly_nonoverlapping(h_begin, h_it)); + return h_it; +} + +template +struct greater_than_or_equal +{ + template + using fn = boost::mp11::mp_bool< Lhs::value >= Rhs::value >; +}; + +constexpr int expansion_sum_length(int s1, int s2) +{ + return s1 != -1 && s2 != -1 ? s1 + s2 : -1; +} + +constexpr int expansion_product_length(int s1, int s2, bool same = false) +{ + if (s1 == -1 || s2 == -1) + { + return -1; + } + else if (same && s1 == 2 && s2 == 2) + { + return 6; + } + return 2 * s1 * s2; +} + +template +using default_zero_elimination_policy = + boost::mp11::mp_bool<(Length > 16) || Length == -1>; + +template +using default_fast_expansion_sum_policy = + boost::mp11::mp_bool< (Length1 > 2 || Length1 == -1) + && (Length2 > 2 || Length2 == -2)>; + +template +< + int ELength, + int FLength, + bool Inplace = false, + bool ENegate = false, + bool FNegate = false, + template class ZeroElimination = default_zero_elimination_policy, + bool FastExpansion = default_fast_expansion_sum_policy + < + ELength, + FLength + >::value, + bool MostSigOnly = false, + int HLength = expansion_sum_length(ELength, FLength) +> +struct expansion_plus_impl +{ + template + static constexpr OutIter apply( + InIter1 e_begin, + InIter1 e_end, + InIter2 f_begin, + InIter2 f_end, + OutIter h_begin, + OutIter h_end) + { + return fast_expansion_sum + < + ZeroElimination::value, + MostSigOnly, + InIter1, + InIter2, + OutIter, + Inplace, + ENegate, + FNegate + >(e_begin, e_end, f_begin, f_end, h_begin, h_end); + } +}; + +template +< + bool Inplace, + bool ENegate, + bool FNegate, + template class ZEP, + bool MostSigOnly, + int HLength +> +struct expansion_plus_impl + < + 1, + 1, + Inplace, + ENegate, + FNegate, + ZEP, + false, + MostSigOnly, + HLength + > +{ + template + static constexpr OutIter apply(InIter1 e_begin, + InIter1, + InIter2 f_begin, + InIter2, + OutIter h_begin, + OutIter) + { + auto x = negate(*e_begin) + negate(*f_begin); + auto y = two_sum_tail(negate(*e_begin), + negate(*f_begin), x); + auto h_it = h_begin; + h_it = insert_ze::value>(h_it, y); + h_it = insert_ze_final::value>(h_it, h_begin, x); + return h_it; + } +}; + +template +< + int ELength, + bool Inplace, + bool ENegate, + bool FNegate, + template class ZEP, + bool MostSigOnly, + int HLength +> +struct expansion_plus_impl + < + ELength, + 1, + Inplace, + ENegate, + FNegate, + ZEP, + false, + MostSigOnly, + HLength + > +{ + template + static constexpr OutIter apply(InIter1 e_begin, + InIter1 e_end, + InIter2 f_begin, + InIter2, + OutIter h_begin, + OutIter h_end) + { + return grow_expansion + < + ZEP::value, + MostSigOnly, + InIter1, + OutIter, + typename std::iterator_traits::value_type, + ENegate, + FNegate + >(e_begin, e_end, *f_begin, h_begin, h_end); + } +}; + +template +< + int FLength, + bool Inplace, + bool ENegate, + bool FNegate, + template class ZEP, + bool FE, + bool MostSigOnly, + int HLength +> +struct expansion_plus_impl + < + 1, + FLength, + Inplace, + ENegate, + FNegate, + ZEP, + FE, + MostSigOnly, + HLength + > +{ + template + static constexpr OutIter apply(InIter1 e_begin, + InIter1, + InIter2 f_begin, + InIter2 f_end, + OutIter h_begin, + OutIter h_end) + { + return grow_expansion + < + ZEP::value, + MostSigOnly, + InIter2, + OutIter, + typename std::iterator_traits::value_type, + ENegate, + FNegate + >(f_begin, f_end, *e_begin, h_begin, h_end); + } +}; + +template +< + int FLength, + int ELength, + bool Inplace, + bool ENegate, + bool FNegate, + template class ZEP, + int HLength +> +struct expansion_plus_impl + < + FLength, + ELength, + Inplace, + ENegate, + FNegate, + ZEP, + false, + false, + HLength + > +{ + template + static constexpr OutIter apply(InIter1 e_begin, + InIter1 e_end, + InIter2 f_begin, + InIter2 f_end, + OutIter h_begin, + OutIter h_end) + { + return expansion_sum + < + ZEP::value, + InIter1, + InIter2, + OutIter, + ENegate, + FNegate + >(e_begin, e_end, f_begin, f_end, h_begin, h_end); + } +}; + +template +< + int FLength, + int ELength, + bool Inplace, + bool ENegate, + bool FNegate, + template class ZEP, + int HLength +> +struct expansion_plus_impl + < + FLength, + ELength, + Inplace, + ENegate, + FNegate, + ZEP, + false, + true, + HLength + > +{ + template + static constexpr OutIter apply(InIter1 e_begin, + InIter1 e_end, + InIter2 f_begin, + InIter2 f_end, + OutIter h_begin, + OutIter h_end) + { + return expansion_sum_most_sig + < + ZEP::value, + InIter1, + InIter2, + OutIter, + ENegate, + FNegate + >(e_begin, e_end, f_begin, f_end, h_begin, h_end); + } +}; + + + +template +struct zero_init +{ + template + static constexpr void apply(Real& a) + { + a = 0; + } +}; + +template <> struct zero_init +{ + template + static constexpr void apply(Real& a) {} +}; + + +template +< + int ELength, + int FLength, + bool Inplace, + template class ZEP = default_zero_elimination_policy, + template class FEP = default_fast_expansion_sum_policy, + bool MostSigOnly = false, + typename InIter1, + typename InIter2, + typename OutIter, + int Result = expansion_sum_length(ELength, FLength) +> +constexpr OutIter expansion_plus(InIter1 e_begin, + InIter1 e_end, + InIter2 f_begin, + InIter2 f_end, + OutIter h_begin, + OutIter h_end) +{ + zero_init::apply(*h_begin); + return expansion_plus_impl + < + ELength, + FLength, + Inplace, + false, + false, + ZEP, + FEP::value, + MostSigOnly + >::apply(e_begin, e_end, f_begin, f_end, h_begin, h_end); +} + +template +< + int ELength, + int FLength, + bool Inplace, + template class ZEP = default_zero_elimination_policy, + template class = default_fast_expansion_sum_policy, + bool MostSigOnly = false, + typename InIter, + typename Real, + typename OutIter, + int Result = expansion_sum_length(ELength, FLength) +> +constexpr OutIter expansion_plus(InIter e_begin, + InIter e_end, + Real f, + OutIter h_begin, + OutIter h_end) +{ + zero_init::apply(*h_begin); + static_assert(FLength == 1, + "FLength must be 1 if f is a single component."); + return grow_expansion + < + ZEP::value, + MostSigOnly + >(e_begin, e_end, f, h_begin, h_end); +} + +template +< + int ELength, + int FLength, + bool Inplace, + template class ZEP = default_zero_elimination_policy, + template class = default_fast_expansion_sum_policy, + bool MostSigOnly = false, + typename InIter, + typename Real, + typename OutIter, + int Result = expansion_sum_length(ELength, FLength) +> +constexpr OutIter expansion_plus(Real e, + InIter f_begin, + InIter f_end, + OutIter h_begin, + OutIter h_end) +{ + zero_init::apply(*h_begin); + static_assert(ELength == 1, + "ELength must be 1 if e is a single component."); + return grow_expansion + < + ZEP::value, + MostSigOnly + >(f_begin, f_end, e, h_begin, h_end); +} + +template +< + int ELength, + int FLength, + bool Inplace, + template class ZEP = default_zero_elimination_policy, + template class = default_fast_expansion_sum_policy, + bool MostSigOnly = false, + typename Real, + typename OutIter, + int Result = expansion_sum_length(ELength, FLength) +> +constexpr OutIter expansion_plus( + Real e, + Real f, + OutIter h_begin, + OutIter) +{ + zero_init::apply(*h_begin); + static_assert(FLength == 1 && ELength == 1, + "ELength and FLength must be 1 if they are single components."); + Real x = e + f; + if (MostSigOnly) + { + *h_begin = x; + return h_begin + 1; + } + Real y = two_sum_tail(e, f, *(h_begin + 1)); + OutIter h_it = h_begin; + h_it = insert_ze::value>(h_it, y); + h_it = insert_ze_final::value>(h_it, x); + return h_it; +} + +template +< + int ELength, + int FLength, + bool Inplace, + bool StageB, + template class ZEP = default_zero_elimination_policy, + template class FEP = default_fast_expansion_sum_policy, + bool MostSigOnly = false, + typename InIter1, + typename InIter2, + typename OutIter, + int Result = expansion_sum_length(ELength, FLength) +> +constexpr OutIter expansion_minus(InIter1 e_begin, + InIter1 e_end, + InIter2 f_begin, + InIter2 f_end, + OutIter h_begin, + OutIter h_end) +{ + zero_init::apply(*h_begin); + return expansion_plus_impl + < + ELength, + FLength, + Inplace, + false, + true, + ZEP, + FEP::value, + MostSigOnly + >::apply(e_begin, e_end, f_begin, f_end, h_begin, h_end); +} + +template +< + int ELength, + int FLength, + bool Inplace, + bool StageB, + template class ZEP = default_zero_elimination_policy, + template class FEP = default_fast_expansion_sum_policy, + bool MostSigOnly = false, + typename InIter, + typename Real, + typename OutIter, + int Result = expansion_sum_length(ELength, FLength) +> +constexpr OutIter expansion_minus(InIter e_begin, + InIter e_end, + Real f, + OutIter h_begin, + OutIter h_end) +{ + static_assert(FLength == 1, + "FLength must be 1 if f is a single component."); + return expansion_plus + < + ELength, + FLength, + Inplace, + false, + true, + ZEP, + FEP, + MostSigOnly + >(e_begin, e_end, f, h_begin, h_end); +} + +template +< + int ELength, + int FLength, + bool Inplace, + bool StageB, + template class ZEP = default_zero_elimination_policy, + template class = default_fast_expansion_sum_policy, + bool MostSigOnly = false, + typename Real, + typename InIter, + typename OutIter, + int Result = expansion_sum_length(ELength, FLength) +> +constexpr OutIter expansion_minus(Real e, + InIter f_begin, + InIter f_end, + OutIter h_begin, + OutIter h_end) +{ + zero_init::apply(*h_begin); + static_assert(ELength == 1, + "ELength must be 1 if e is a single component."); + return grow_expansion + < + ZEP::value, + MostSigOnly, + InIter, + OutIter, + Real, + false, + true + >(f_begin, f_end, e, h_begin, h_end); +} + +template +< + int ELength, + int FLength, + bool Inplace, + bool StageB, + template class ZEP = default_zero_elimination_policy, + template class = default_fast_expansion_sum_policy, + bool MostSigOnly = false, + typename Real, + typename OutIter, + int Result = expansion_sum_length(ELength, FLength) +> +constexpr std::enable_if_t expansion_minus(Real e, + Real f, + OutIter h_begin, + OutIter) +{ + static_assert(ELength == 1 && FLength == 1, + "ELength and FLength must be 1 if they are single components."); + Real x = e - f; + if (MostSigOnly) + { + *h_begin = x; + return h_begin + 1; + } + Real y = two_difference_tail(e, f, x); + OutIter h_it = h_begin; + h_it = insert_ze::value>(h_it, y); + h_it = insert_ze_final::value>(h_it, h_begin, x); + return h_it; +} + +template +< + int ELength, + int FLength, + bool Inplace, + bool StageB, + template class = default_zero_elimination_policy, + template class = default_fast_expansion_sum_policy, + bool MostSigOnly = false, + typename Real, + typename OutIter, + int Result = expansion_sum_length(ELength, FLength) +> +constexpr std::enable_if_t expansion_minus(Real e, + Real f, + OutIter h_begin, + OutIter) +{ + static_assert(ELength == 1 && FLength == 1, + "ELength and FLength must be 1 if they are single components."); + *(h_begin) = e - f; + return h_begin + 1; +} + +template +< + int ELength, + int FLength, + bool Inplace, + template class ZEP = default_zero_elimination_policy, + template class = default_fast_expansion_sum_policy, + bool LeftEqualsRight, + typename InIter, + typename Real, + typename OutIter, + int Result = expansion_product_length(ELength, FLength, LeftEqualsRight) +> +constexpr OutIter expansion_times(InIter e_begin, + InIter e_end, + Real f, + OutIter h_begin, + OutIter h_end) +{ + static_assert(FLength == 1, + "FLength must be 1 if f is a single component."); + return scale_expansion + < + ZEP::value + >(e_begin, e_end, f, h_begin, h_end); +} + +template +< + int ELength, + int FLength, + bool Inplace, + template class ZEP = default_zero_elimination_policy, + template class = default_fast_expansion_sum_policy, + bool LeftEqualsRight, + typename InIter, + typename Real, + typename OutIter, + int Result = expansion_product_length(ELength, FLength, LeftEqualsRight) +> +constexpr OutIter expansion_times(Real e, + InIter f_begin, + InIter f_end, + OutIter h_begin, + OutIter h_end) +{ + static_assert(ELength == 1, + "ELength must be 1 if e is a single component."); + auto h_it = scale_expansion + < + ZEP::value + >(f_begin, f_end, e, h_begin, h_end); + return h_it; +} + +template +< + int ELength, + int FLength, + bool Inplace, + template class ZEP = default_zero_elimination_policy, + template class = default_fast_expansion_sum_policy, + bool LeftEqualsRight, + typename Real, + typename OutIter, + int Result = expansion_product_length(ELength, FLength, LeftEqualsRight) +> +constexpr OutIter expansion_times(Real e, + Real f, + OutIter h_begin, + OutIter) +{ + static_assert(ELength == 1 && FLength == 1, + "ELength and FLength must be 1 if they are single components."); + Real x = e * f; + Real y = two_product_tail(e, f, x); + OutIter h_it = h_begin; + h_it = insert_ze::value>(h_it, y); + h_it = insert_ze_final::value>(h_it, h_begin, x); + return h_it; +} + +template +< + int ELength, + int FLength, + bool, + template class = default_zero_elimination_policy, + template class = default_fast_expansion_sum_policy, + bool LeftEqualsRight, + typename In1, + typename In2, + typename Out, + int = expansion_product_length(ELength, FLength, LeftEqualsRight) +> +constexpr Out expansion_times(In1, In1, In2, In2, Out, Out); + +#if defined(BOOST_CLANG) +#pragma clang diagnostic push +#pragma clang diagnostic ignored "-Wunused-variable" +#endif + +#if defined(BOOST_GCC) && (BOOST_GCC >= 40900) +#pragma GCC diagnostic push +#pragma GCC diagnostic ignored "-Wunused-variable" +#endif + +template +< + int ELength, + int FLength, + template class ZeroEliminationPolicy = + default_zero_elimination_policy, + template class FastExpansionSumPolicy = + default_fast_expansion_sum_policy, + bool LeftEqualsRight = false, + int Result = expansion_product_length(ELength, FLength, LeftEqualsRight), + typename ESmaller = boost::mp11::mp_bool +> +struct expansion_times_impl +{ + template + static constexpr OutIter apply(InIter1 e_begin, + InIter1 e_end, + InIter2 f_begin, + InIter2 f_end, + OutIter h_begin, + OutIter h_end) + { + assert(debug_expansion::expansion_nonoverlapping(e_begin, e_end)); + assert(debug_expansion::expansion_nonoverlapping(f_begin, f_end)); + const auto e_dyn_length = std::distance(e_begin, e_end); + const auto h_old_end = h_end; + if (e_dyn_length == 1) + { + auto h_it = expansion_times + < + 1, + FLength, + false, + ZeroEliminationPolicy, + FastExpansionSumPolicy, + LeftEqualsRight + >(*e_begin, f_begin, f_end, h_begin, h_end); + assert(debug_expansion::expansion_nonoverlapping(h_begin, h_it)); + assert(std::abs( + std::accumulate(e_begin, e_end, 0.) + * std::accumulate(f_begin, f_end, 0.) + - std::accumulate(h_begin, h_it, 0.)) <= + std::abs(std::accumulate(h_begin, h_it, 0.) * 1e-10)); + return h_it; + } + else if (e_dyn_length > 1 ) + { + constexpr int e_length1 = ELength == -1 ? -1 : ELength / 2; + auto e_mid = e_begin + e_dyn_length / 2; + auto h_mid = h_begin + + std::distance(e_begin, e_mid) + * std::distance(f_begin, f_end) * 2; + h_mid = expansion_times + < + e_length1, + FLength, + false, + ZeroEliminationPolicy, + FastExpansionSumPolicy, + LeftEqualsRight + >(e_begin, e_mid, f_begin, f_end, h_begin, h_mid); + constexpr int e_length2 = ELength == -1 ? + -1 + : ELength - e_length1; + h_end = expansion_times + < + e_length2, + FLength, + false, + ZeroEliminationPolicy, + FastExpansionSumPolicy, + LeftEqualsRight + >(e_mid, e_end, f_begin, f_end, h_mid, h_end); + + constexpr int summand_length1 = + expansion_product_length(e_length1, FLength, false); + constexpr int summand_length2 = + expansion_product_length(e_length2, FLength, false); + OutIter h_it = expansion_plus + < + summand_length1, + summand_length2, + true, + ZeroEliminationPolicy, + FastExpansionSumPolicy + >(h_begin, h_mid, h_mid, h_end, h_begin, h_end); + + assert(debug_expansion::expansion_nonoverlapping(h_begin, h_it)); + assert(std::abs( + std::accumulate(e_begin, e_end, 0.) + * std::accumulate(f_begin, f_end, 0.) + - std::accumulate(h_begin, h_it, 0.)) <= + std::abs(std::accumulate(h_begin, h_it, 0.) * 1e-10)); + return h_it; + } + else if (e_dyn_length == 0) + { + return h_begin; + } + assert(false); + return h_begin; + } +}; + +#if defined(BOOST_CLANG) +#pragma clang diagnostic pop +#endif + +#if defined(BOOST_GCC) && (BOOST_GCC >= 40900) +#pragma GCC diagnostic pop +#endif + +template +< + bool ZeroElimination +> +struct advance_ze_impl +{ + template + Iter apply(Iter it) + { + return *it == 0 ? it : it + 1; + } +}; + +template <> +struct advance_ze_impl +{ + template + Iter apply(Iter it) + { + return it + 1; + } +}; + +template +struct constant_ze +{ + template + using fn = boost::mp11::mp_bool; +}; + +template +< + bool ZE, + bool EZeroEliminated +> +struct two_square_impl +{ + template + static constexpr OutIter apply(InIter e_begin, + InIter e_end, + OutIter h_begin, + OutIter h_end) + { + std::array, 5> cache; + cache[2] = *(e_begin) * *(e_begin); + auto h_it = h_begin; + h_it = insert_ze(h_it, + two_product_tail(*e_begin, *e_begin, cache[2])); + cache[1] = *(e_begin + 1) * 2.0 * *(e_begin); + cache[0] = two_product_tail( *( e_begin + 1 ), + 2.0 * *(e_begin), + cache[1] ); + expansion_plus<2, 1, false>(cache.cbegin(), + cache.cbegin() + 2, + cache[2], + cache.begin() + 2, + cache.begin() + 5); + h_it = insert_ze(h_it, cache[2]); + cache[1] = *(e_begin + 1) * *(e_begin + 1); + cache[0] = two_product_tail( *(e_begin + 1), + *(e_begin + 1), + cache[1]); + h_it = expansion_plus + < + 2, + 2, + false, + constant_ze::template fn + >(cache.cbegin(), + cache.cbegin() + 2, + cache.cbegin() + 3, + cache.cbegin() + 5, + h_it, + h_it + 4); + return h_it; + } +}; + +template <> +struct two_square_impl +{ + template + static constexpr OutIter apply(InIter e_begin, + InIter e_end, + OutIter h_begin, + OutIter h_end) + { + if (e_begin + 1 == e_end) + { + *(h_begin + 1) = *e_begin * *e_begin; + *h_begin = two_product_tail(*e_begin, *e_begin, *(h_begin + 1)); + return h_begin + 2; + } + else + { + return two_square_impl + ::apply(e_begin, e_end, h_begin, h_end); + } + } +}; + +template +< + template class ZEP, + template class FEP +> +struct expansion_times_impl + < + 2, + 2, + ZEP, + FEP, + true, + 6, + boost::mp11::mp_true + > +{ + template + static constexpr OutIter apply(InIter1 e_begin, + InIter1 e_end, + InIter2, + InIter2, + OutIter h_begin, + OutIter h_end) + { + return two_square_impl + < + ZEP<6>::value, + ZEP<2>::value + >::apply(e_begin, e_end, h_begin, h_end); + } +}; + + +template +< + int ELength, + int FLength, + template class ZEP, + template class FEP, + bool LeftEqualsRight, + int Result +> +struct expansion_times_impl + < + ELength, + FLength, + ZEP, + FEP, + LeftEqualsRight, + Result, + boost::mp11::mp_false + > +{ + template + static constexpr OutIter apply(InIter1 e_begin, + InIter1 e_end, + InIter2 f_begin, + InIter2 f_end, + OutIter h_begin, + OutIter h_end) + { + return expansion_times_impl + < + FLength, + ELength, + ZEP, + FEP, + LeftEqualsRight + >::apply(f_begin, + f_end, + e_begin, + e_end, + h_begin, + h_end); + } +}; + +template +< + template class ZEP, + template class FEP, + bool LeftEqualsRight +> +struct expansion_times_impl + < + 1, + 1, + ZEP, + FEP, + LeftEqualsRight, + 2 + > +{ + template + static constexpr OutIter apply(InIter1 e_begin, + InIter1, + InIter2 f_begin, + InIter2, + OutIter h_begin, + OutIter) + { + auto x = *e_begin * *f_begin; + auto y = two_product_tail(*e_begin, *f_begin, x); + OutIter h_it = h_begin; + h_it = insert_ze::value>(h_it, y); + h_it = insert_ze_final::value>(h_it, h_begin, x); + return h_it; + } +}; + +template +< + int ELength, + int FLength, + bool Inplace, + template class ZEP, + template class FEP, + bool LeftEqualsRight, + typename InIter1, + typename InIter2, + typename OutIter, + int Result +> +constexpr OutIter expansion_times(InIter1 e_begin, + InIter1 e_end, + InIter2 f_begin, + InIter2 f_end, + OutIter h_begin, + OutIter h_end) +{ + return expansion_times_impl + < + ELength, + FLength, + ZEP, + FEP, + LeftEqualsRight, + Result + >::apply(e_begin, e_end, f_begin, f_end, h_begin, h_end); +} + +template +< + typename Iter +> +constexpr Iter compress(Iter e_begin, Iter e_end) +{ + if (std::distance(e_begin, e_end) < 2) return e_end; + Iter bottom_it = std::prev(e_end); + auto Q = *bottom_it; + Iter e_it; + for (e_it = std::prev(bottom_it);; --e_it) + { + auto Q_next = Q + *e_it; + auto q = two_sum_tail(Q, *e_it, Q_next); + Q = Q_next; + if (q != 0) + { + *bottom_it = Q; + --bottom_it; + Q = q; + } + if (e_it == e_begin) break; + } + *bottom_it = Q; + Iter top_it = e_begin; + for(e_it = std::next(bottom_it); e_it != e_end; ++e_it) + { + auto Q_next = *e_it + Q; + auto q = fast_two_sum_tail(*e_it, Q, Q_next); + Q = Q_next; + if (q != 0) + { + *top_it = q; + ++top_it; + } + } + *top_it = Q; + ++top_it; + return top_it; +} + +}} // namespace detail::generic_robust_predicates + +}} // namespace boost::geometry + +#endif // BOOST_GEOMETRY_EXTENSIONS_GENERIC_ROBUST_PREDICATES_STRATEGIES_CARTESIAN_DETAIL_EXPANSION_ARITHMETIC_HPP diff --git a/include/boost/geometry/extensions/generic_robust_predicates/strategies/cartesian/detail/expansion_eval.hpp b/include/boost/geometry/extensions/generic_robust_predicates/strategies/cartesian/detail/expansion_eval.hpp new file mode 100644 index 0000000000..e9b282f1af --- /dev/null +++ b/include/boost/geometry/extensions/generic_robust_predicates/strategies/cartesian/detail/expansion_eval.hpp @@ -0,0 +1,853 @@ +// Boost.Geometry (aka GGL, Generic Geometry Library) + +// Copyright (c) 2020 Tinko Bartels, Berlin, Germany. + +// Contributed and/or modified by Tinko Bartels, +// as part of Google Summer of Code 2020 program. + +// Use, modification and distribution is subject to 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) + +#ifndef BOOST_GEOMETRY_EXTENSIONS_GENERIC_ROBUST_PREDICATES_STRATEGIES_CARTESIAN_DETAIL_EXPANSION_EVAL_HPP +#define BOOST_GEOMETRY_EXTENSIONS_GENERIC_ROBUST_PREDICATES_STRATEGIES_CARTESIAN_DETAIL_EXPANSION_EVAL_HPP + +#include +#include + +#include +#include +#include + +#include +#include + +namespace boost { namespace geometry +{ + +namespace detail { namespace generic_robust_predicates +{ + +template +< + typename Expression, + bool StageB = false, + operator_types Op = Expression::operator_type +> +struct expansion_size_impl {}; + +template +struct expansion_size_impl +{ + static constexpr std::size_t value = 1; +}; + +template +struct expansion_size_impl +{ +private: + static constexpr std::size_t left_size = + expansion_size_impl::value; + static constexpr std::size_t right_size = + expansion_size_impl::value; +public: + static constexpr std::size_t value = left_size + right_size; +}; + +template +struct expansion_size_impl +{ +private: + static constexpr std::size_t left_size = + expansion_size_impl::value; + static constexpr std::size_t right_size = + expansion_size_impl::value; +public: + static constexpr std::size_t value = + StageB && Expression::left::is_leaf && Expression::right::is_leaf ? + 1 + : left_size + right_size; +}; + +template +struct expansion_size_impl +{ +private: + static constexpr std::size_t left_size = + expansion_size_impl::value; + static constexpr std::size_t right_size = + expansion_size_impl::value; +public: + static constexpr std::size_t value = + std::is_same + < + typename Expression::left, + typename Expression::right + >::value && left_size == 2 && right_size == 2 ? + 6 + : 2 * left_size * right_size; +}; + +template using expansion_size = + boost::mp11::mp_size_t< expansion_size_impl::value >; + +template using expansion_size_stage_b = + boost::mp11::mp_size_t< expansion_size_impl::value >; + +template +< + typename Expression, + std::size_t LeftLength, + std::size_t RightLength, + operator_types Op = Expression::operator_type +> +constexpr std::size_t final_expansion_size() +{ + if (Op == operator_types::product) + { + return 2 * LeftLength * RightLength; + } + return 1; +}; + +template using no_zero_elimination_policy = boost::mp11::mp_false; + +template +< + operator_types Op, + std::size_t LeftLength, + std::size_t RightLength, + bool Inplace, + typename Iter, + template class ZEPolicy, + template class FEPolicy, + bool StageB = false, + bool LeftEqualsRight = false, + bool MostSigOnly = false +> +struct perform_op_impl {}; + +template +< + std::size_t LeftLength, + std::size_t RightLength, + bool Inplace, + typename Iter, + template class ZEPolicy, + template class FEPolicy, + bool StageB, + bool LeftEqualsRight, + bool MostSigOnly +> +struct perform_op_impl + < + operator_types::sum, + LeftLength, + RightLength, + Inplace, + Iter, + ZEPolicy, + FEPolicy, + StageB, + LeftEqualsRight, + MostSigOnly + > +{ + template + static constexpr Iter apply(Args...args) + { + return expansion_plus + < + LeftLength, + RightLength, + Inplace, + ZEPolicy, + FEPolicy, + MostSigOnly + >(args...); + } +}; + +template +< + std::size_t LeftLength, + std::size_t RightLength, + bool Inplace, + typename Iter, + template class ZEPolicy, + template class FEPolicy, + bool StageB, + bool LeftEqualsRight, + bool MostSigOnly +> +struct perform_op_impl + < + operator_types::difference, + LeftLength, + RightLength, + Inplace, + Iter, + ZEPolicy, + FEPolicy, + StageB, + LeftEqualsRight, + MostSigOnly + > +{ + template + static constexpr Iter apply(Args...args) + { + return expansion_minus + < + LeftLength, + RightLength, + Inplace, + StageB, + ZEPolicy, + FEPolicy, + MostSigOnly + >(args...); + } +}; + +template +< + std::size_t LeftLength, + std::size_t RightLength, + bool Inplace, + typename Iter, + template class ZEPolicy, + template class FEPolicy, + bool StageB, + bool LeftEqualsRight, + bool MostSigOnly +> +struct perform_op_impl + < + operator_types::product, + LeftLength, + RightLength, + Inplace, + Iter, + ZEPolicy, + FEPolicy, + StageB, + LeftEqualsRight, + MostSigOnly + > +{ + template + static constexpr Iter apply(Args...args) + { + return expansion_times + < + LeftLength, + RightLength, + Inplace, + ZEPolicy, + FEPolicy, + LeftEqualsRight + >(args...); + } +}; + +template +< + typename Evals, + typename Eval, + typename Sizes, + typename AccumulatedSizes, + typename ZEEvals, + typename Iter, + typename Real, + template class ZEPolicy, + template class FEPolicy, + bool StageB = false, + operator_types Op = Eval::operator_type, + bool LeftLeaf = Eval::left::is_leaf, + bool RightLeaf = Eval::right::is_leaf, + bool LeftEqualsRight = false, + bool MostSigOnly = false +> +struct eval_expansion_impl {}; + +template +struct set_exp_end +{ + template + static inline void apply(Iter const&, IterArr&) {} +}; + +template +struct set_exp_end +{ + template + static inline void apply(Iter const& it, IterArr& ze_ends) + { + ze_ends[I] = it; + } +}; + +template +< + typename Evals, + typename Eval, + typename Sizes, + typename AccumulatedSizes, + typename ZEEvals, + typename Iter, + typename Real, + template class ZEPolicy, + template class FEPolicy, + bool StageB, + operator_types Op, + bool LeftEqualsRight, + bool MostSigOnly +> +struct eval_expansion_impl + < + Evals, + Eval, + Sizes, + AccumulatedSizes, + ZEEvals, + Iter, + Real, + ZEPolicy, + FEPolicy, + StageB, + Op, + true, + true, + LeftEqualsRight, + MostSigOnly + > +{ +private: + using left = typename Eval::left; + using right = typename Eval::right; + using eval_index = boost::mp11::mp_find; + static constexpr std::size_t size = + boost::mp11::mp_at::value; + static constexpr std::size_t start = + boost::mp11::mp_at::value; +public: + template + static constexpr Iter apply(Iter begin, + Iter, + IterArr& ze_ends, + Reals const&... args) + { + std::array input + {{ static_cast(args)... }}; + Real left_val = input[left::argn - 1]; + Real right_val = input[right::argn - 1]; + auto end = perform_op_impl + < + Op, + 1, + 1, + false, + Iter, + ZEPolicy, + FEPolicy, + StageB, + LeftEqualsRight, + MostSigOnly + >::apply(left_val, + right_val, + begin + start, + begin + start + size); + set_exp_end + < + ( boost::mp11::mp_find::value + < std::tuple_size::value), + boost::mp11::mp_find::value + >::apply(end, ze_ends); + return end; + } +}; + +template +struct get_exp_end +{ + template + static inline Iter apply(Iter default_end, + IterArr const& ze_ends) + { + return default_end; + } +}; + +template struct get_exp_end +{ + template + static inline Iter apply(Iter, IterArr const& ze_ends) + { + return ze_ends[I]; + } +}; + +template +< + typename Evals, + typename Eval, + typename Sizes, + typename AccumulatedSizes, + typename ZEEvals, + typename Iter, + typename Real, + template class ZEPolicy, + template class FEPolicy, + bool StageB, + operator_types Op, + bool LeftEqualsRight, + bool MostSigOnly +> +struct eval_expansion_impl + < + Evals, + Eval, + Sizes, + AccumulatedSizes, + ZEEvals, + Iter, + Real, + ZEPolicy, + FEPolicy, + StageB, + Op, + true, + false, + LeftEqualsRight, + MostSigOnly + > +{ +private: + using left = typename Eval::left; + using right = typename Eval::right; + using eval_index = boost::mp11::mp_find; + static constexpr std::size_t size = + boost::mp11::mp_at::value; + static constexpr std::size_t start = + boost::mp11::mp_at::value; + using right_eval_index = boost::mp11::mp_find; + static constexpr std::size_t right_size = + boost::mp11::mp_at::value; + static constexpr std::size_t right_start = + boost::mp11::mp_at::value; +public: + template + static constexpr Iter apply(Iter begin, + Iter, + IterArr& ze_ends, + Reals const&... args) + { + std::array input + {{ static_cast(args)... }}; + Real left_val = input[left::argn - 1]; + Iter right_end = + get_exp_end + < + ZEPolicy::value, + boost::mp11::mp_find::value + >::apply(begin + right_start + right_size, ze_ends); + Iter end = perform_op_impl + < + Op, + 1, + right_size, + false, + Iter, + ZEPolicy, + FEPolicy, + StageB, + LeftEqualsRight + >::apply( + left_val, + begin + right_start, + right_end, + begin + start, + begin + start + size); + set_exp_end + < + ( boost::mp11::mp_find::value + < std::tuple_size::value), + boost::mp11::mp_find::value + >::apply(end, ze_ends); + return end; + } +}; + +template +< + typename Evals, + typename Eval, + typename Sizes, + typename AccumulatedSizes, + typename ZEEvals, + typename Iter, + typename Real, + template class ZEPolicy, + template class FEPolicy, + bool StageB, + operator_types Op, + bool LeftEqualsRight, + bool MostSigOnly +> +struct eval_expansion_impl + < + Evals, + Eval, + Sizes, + AccumulatedSizes, + ZEEvals, + Iter, + Real, + ZEPolicy, + FEPolicy, + StageB, + Op, + false, + true, + LeftEqualsRight, + MostSigOnly + > +{ +private: + using left = typename Eval::left; + using right = typename Eval::right; + using eval_index = boost::mp11::mp_find; + static constexpr std::size_t size = + boost::mp11::mp_at::value; + static constexpr std::size_t start = + boost::mp11::mp_at::value; + using left_eval_index = boost::mp11::mp_find; + static constexpr std::size_t left_size = + boost::mp11::mp_at::value; + static constexpr std::size_t left_start = + boost::mp11::mp_at::value; +public: + template + static constexpr Iter apply(Iter begin, + Iter, + IterArr& ze_ends, + Reals const&... args) + { + std::array input + {{ static_cast(args)... }}; + Real right_val = input[right::argn - 1]; + Iter left_end = + get_exp_end + < + ZEPolicy::value, + boost::mp11::mp_find::value + >::apply(begin + left_start + left_size, ze_ends); + auto end = perform_op_impl + < + Op, + left_size, + 1, + false, + Iter, + ZEPolicy, + FEPolicy, + StageB, + LeftEqualsRight + >::apply( + begin + left_start, + left_end, + right_val, + begin + start, + begin + start + size); + set_exp_end + < + ( boost::mp11::mp_find::value + < std::tuple_size::value), + boost::mp11::mp_find::value + >::apply(end, ze_ends); + return end; + } +}; + +template +< + typename Evals, + typename Eval, + typename Sizes, + typename AccumulatedSizes, + typename ZEEvals, + typename Iter, + typename Real, + template class ZEPolicy, + template class FEPolicy, + bool StageB, + operator_types Op, + bool LeftEqualsRight, + bool MostSigOnly +> +struct eval_expansion_impl + < + Evals, + Eval, + Sizes, + AccumulatedSizes, + ZEEvals, + Iter, + Real, + ZEPolicy, + FEPolicy, + StageB, + Op, + false, + false, + LeftEqualsRight, + MostSigOnly + > +{ +private: + using left = typename Eval::left; + using right = typename Eval::right; + using eval_index = boost::mp11::mp_find; + static constexpr std::size_t size = + boost::mp11::mp_at::value; + static constexpr std::size_t start = + boost::mp11::mp_at::value; + using left_eval_index = + boost::mp11::mp_find; + static constexpr std::size_t left_size = + boost::mp11::mp_at::value; + static constexpr std::size_t left_start = + boost::mp11::mp_at::value; + using right_eval_index = boost::mp11::mp_find; + static constexpr std::size_t right_size = + boost::mp11::mp_at::value; + static constexpr std::size_t right_start = + boost::mp11::mp_at::value; +public: + template + static constexpr Iter apply(Iter begin, Iter, + IterArr& ze_ends, + Reals const&...) + { + Iter left_end = + get_exp_end + < + ( boost::mp11::mp_find::value + < std::tuple_size::value), + boost::mp11::mp_find::value + >::apply(begin + left_start + left_size, ze_ends); + Iter right_end = + get_exp_end + < + ( boost::mp11::mp_find::value + < std::tuple_size::value), + boost::mp11::mp_find::value + >::apply(begin + right_start + right_size, ze_ends); + auto end = perform_op_impl + < + Op, + left_size, + right_size, + false, + Iter, + ZEPolicy, + FEPolicy, + StageB, + LeftEqualsRight, + MostSigOnly + >::apply( + begin + left_start, + left_end, + begin + right_start, + right_end, + begin + start, + begin + start + size); + set_exp_end + < + ( boost::mp11::mp_find::value + < std::tuple_size::value), + boost::mp11::mp_find::value + >::apply(end, ze_ends); + return end; + } +}; + +template +< + typename Evals, + typename RemainingEvals, + typename Sizes, + typename AccumulatedSizes, + typename ZEEvals, + typename Iter, + typename Real, + template class ZEPolicy, + template class FEPolicy, + bool StageB = false, + std::size_t RemainingSize = boost::mp11::mp_size::value +> +struct eval_expansions_impl +{ + template + static constexpr Iter apply(Iter begin, + Iter end, + IterArr& ze_ends, + Reals const&... args) + { + using eval = boost::mp11::mp_front; + eval_expansion_impl + < + Evals, + eval, + Sizes, + AccumulatedSizes, + ZEEvals, + Iter, + Real, + ZEPolicy, + FEPolicy, + StageB, + eval::operator_type, + eval::left::is_leaf, + eval::right::is_leaf, + std::is_same::value, + false + >::apply(begin, end, ze_ends, args...); + return eval_expansions_impl + < + Evals, + boost::mp11::mp_pop_front, + Sizes, + AccumulatedSizes, + ZEEvals, + Iter, + Real, + ZEPolicy, + FEPolicy, + StageB + >::apply(begin, end, ze_ends, args...); + } +}; + +template +using force_zero_elimination_policy = boost::mp11::mp_true; + +template +< + typename Evals, + typename RemainingEvals, + typename Sizes, + typename AccumulatedSizes, + typename ZEEvals, + typename Iter, + typename Real, + template class ZEPolicy, + template class FEPolicy, + bool StageB +> +struct eval_expansions_impl + < + Evals, + RemainingEvals, + Sizes, + AccumulatedSizes, + ZEEvals, + Iter, + Real, + ZEPolicy, + FEPolicy, + StageB, + 1 + > +{ + template + static constexpr Iter apply(Iter begin, + Iter end, + IterArr& ze_ends, + Reals const&... args) + { + using eval = boost::mp11::mp_front; + + auto new_end = eval_expansion_impl + < + Evals, + eval, + Sizes, + AccumulatedSizes, + ZEEvals, + Iter, + Real, + ZEPolicy, + FEPolicy, + StageB, + eval::operator_type, + eval::left::is_leaf, + eval::right::is_leaf, + std::is_same + < + typename eval::left, + typename eval::right + >::value, + true + >::apply(begin, end, ze_ends, args...); + return new_end; + } +}; + +template +< + template class ZEPolicy, + bool StageB +> +struct is_zero_elim_q +{ + template + using fn = ZEPolicy + < + expansion_size_impl::value + >; +}; + +template +< + typename Evals, + typename Sizes, + typename AccumulatedSizes, + typename Iter, + typename Real, + bool StageB, + template class ZEPolicy, + template class FEPolicy, + typename ...Reals +> +static constexpr Iter eval_expansions(Iter begin, + Iter end, + Reals const&... args) +{ + using ze_evals = boost::mp11::mp_copy_if_q + < + Evals, + is_zero_elim_q + >; + std::array::value> ze_ends; + return eval_expansions_impl + < + Evals, + Evals, + Sizes, + AccumulatedSizes, + ze_evals, + Iter, + Real, + ZEPolicy, + FEPolicy, + StageB + >::apply(begin, end, ze_ends, args...); +} + +}} // namespace detail::generic_robust_predicates + +}} // namespace boost::geometry + +#endif // BOOST_GEOMETRY_EXTENSIONS_GENERIC_ROBUST_PREDICATES_STRATEGIES_CARTESIAN_DETAIL_EXPANSION_EVAL_HPP diff --git a/include/boost/geometry/extensions/generic_robust_predicates/strategies/cartesian/detail/expression_eval.hpp b/include/boost/geometry/extensions/generic_robust_predicates/strategies/cartesian/detail/expression_eval.hpp new file mode 100644 index 0000000000..7da252d64f --- /dev/null +++ b/include/boost/geometry/extensions/generic_robust_predicates/strategies/cartesian/detail/expression_eval.hpp @@ -0,0 +1,281 @@ +// Boost.Geometry (aka GGL, Generic Geometry Library) + +// Copyright (c) 2020 Tinko Bartels, Berlin, Germany. + +// Use, modification and distribution is subject to 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) + +#ifndef BOOST_GEOMETRY_EXTENSIONS_GENERIC_ROBUST_PREDICATES_STRATEGIES_CARTESIAN_DETAIL_EXPRESSION_EVAL_HPP +#define BOOST_GEOMETRY_EXTENSIONS_GENERIC_ROBUST_PREDICATES_STRATEGIES_CARTESIAN_DETAIL_EXPRESSION_EVAL_HPP + +#include +#include +#include + +#include +#include +#include + +#include + +namespace boost { namespace geometry +{ + +namespace detail { namespace generic_robust_predicates +{ + +//The templates in this file are meant to be used for the evaluation of +//expressions with floating-point precision and floating-point rounding. + +template +struct evaluate_expression_unary_impl {}; + +template <> +struct evaluate_expression_unary_impl +{ + template + static constexpr void apply(CT const& c, CT& out) + { + out = std::abs(c); + } +}; + +template +struct evaluate_expression_binary_impl {}; + +template <> +struct evaluate_expression_binary_impl +{ + template + static constexpr void apply(CT const& l, CT const& r, CT& out) + { + out = l + r; + } +}; + +template <> +struct evaluate_expression_binary_impl +{ + template + static constexpr void apply(CT const& l, CT const& r, CT& out) + { + out = l * r; + } +}; + +template <> +struct evaluate_expression_binary_impl +{ + template + static constexpr void apply(CT const& l, CT const& r, CT& out) + { + out = l - r; + } +}; + +template <> +struct evaluate_expression_binary_impl +{ + template + static constexpr void apply(CT const& l, CT const& r, CT& out) + { + out = std::min(l, r); + } +}; + +template <> +struct evaluate_expression_binary_impl +{ + template + static constexpr void apply(CT const& l, CT const& r, CT& out) + { + out = std::max(l, r); + } +}; + +template +< + std::size_t I, + typename Expression, + bool IsLeaf = Expression::is_leaf +> +struct get_arg_or_out_impl {}; + +template +struct get_arg_or_out_impl +{ + template + static constexpr auto apply(In const&, Out const& out) + { + return out[I]; + } +}; + +template +struct get_arg_or_const +{ + template + static constexpr auto apply(In const& in) + { + return in[Argn - 1]; + } +}; + +template +struct get_arg_or_const +{ + template + static constexpr auto apply(In const&) + { + return Expression::value; + } +}; + +template +struct get_arg_or_out_impl +{ + template + static constexpr auto apply(In const& in, Out const&) + { + return get_arg_or_const::apply(in); + } +}; + +template +< + typename InputArr, + typename OutputArr, + typename Expressions, + typename Expression, + operator_arities Op = Expression::operator_arity +> +struct evaluate_expression_impl {}; + +template +< + typename InputArr, + typename OutputArr, + typename Expressions, + typename Expression +> +struct evaluate_expression_impl + < + InputArr, + OutputArr, + Expressions, + Expression, + operator_arities::binary + > +{ + using left = typename Expression::left; + using right = typename Expression::right; + static constexpr std::size_t i_out = + boost::mp11::mp_find::value; + static constexpr std::size_t i_left = + boost::mp11::mp_find::value; + static constexpr std::size_t i_right = + boost::mp11::mp_find::value; + static constexpr void apply( + InputArr const& input, + OutputArr& output + ) + { + auto& out = output[i_out]; + auto const& l = + get_arg_or_out_impl::apply(input, output); + auto const& r = + get_arg_or_out_impl::apply(input, output); + evaluate_expression_binary_impl + ::apply(l, r, out); + } +}; + +template +< + typename InputArr, + typename OutputArr, + typename Expressions, + typename Expression +> +struct evaluate_expression_impl + < + InputArr, + OutputArr, + Expressions, + Expression, + operator_arities::unary + > +{ + using child = typename Expression::child; + static constexpr std::size_t i_out = + boost::mp11::mp_find::value; + static constexpr std::size_t i_child = + boost::mp11::mp_find::value; + static constexpr void apply( + InputArr const& input, + OutputArr& output + ) + { + auto& out = output[i_out]; + auto const& c = get_arg_or_out_impl + < + i_child, + child + >::apply(input, output); + evaluate_expression_unary_impl + ::apply(c, out); + } +}; + +//Expects a list of variadic list of expressions that are to be evaluated. +//The results are stored in output. + +template +< + typename InputArr, + typename OutputArr, + typename ...Expressions +> +constexpr void evaluate_expressions(InputArr const& input, + OutputArr& output, + boost::mp11::mp_list) +{ + using exp_list = boost::mp11::mp_list; + //Below we use a bit of a hack to do a comma-fold in C++14. + using dummy = int[]; + (void)dummy{ + 0, + (evaluate_expression_impl + < + InputArr, + OutputArr, + exp_list, + Expressions + >::apply(input, output), 0)... + }; +} + +template +< + typename Expression, + typename InputArr +> +constexpr auto evaluate_expression(InputArr const& input) +{ + using stack = typename boost::mp11::mp_unique>; + using evals = typename boost::mp11::mp_remove_if; + std::array + < + typename InputArr::value_type, + boost::mp11::mp_size::value + > results; + evaluate_expressions(input, results, evals{}); + return results.back(); +} + +}} // namespace detail::generic_robust_predicates + +}} // namespace boost::geometry + +#endif // BOOST_GEOMETRY_EXTENSIONS_GENERIC_ROBUST_PREDICATES_STRATEGIES_CARTESIAN_DETAIL_EXPRESSION_EVAL_HPP diff --git a/include/boost/geometry/extensions/generic_robust_predicates/strategies/cartesian/detail/expression_tree.hpp b/include/boost/geometry/extensions/generic_robust_predicates/strategies/cartesian/detail/expression_tree.hpp new file mode 100644 index 0000000000..064604a501 --- /dev/null +++ b/include/boost/geometry/extensions/generic_robust_predicates/strategies/cartesian/detail/expression_tree.hpp @@ -0,0 +1,268 @@ +// Boost.Geometry (aka GGL, Generic Geometry Library) + +// Copyright (c) 2020 Tinko Bartels, Berlin, Germany. + +// Contributed and/or modified by Tinko Bartels, +// as part of Google Summer of Code 2020 program. + +// Use, modification and distribution is subject to 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) + +#ifndef BOOST_GEOMETRY_EXTENSIONS_GENERIC_ROBUST_PREDICATES_STRATEGIES_CARTESIAN_DETAIL_EXPRESSION_TREE_HPP +#define BOOST_GEOMETRY_EXTENSIONS_GENERIC_ROBUST_PREDICATES_STRATEGIES_CARTESIAN_DETAIL_EXPRESSION_TREE_HPP + +#include +#include +#include + +#include +#include +#include + +namespace boost { namespace geometry +{ + +namespace detail { namespace generic_robust_predicates +{ + +//The templates in this file are designed to represent arithmetic expressions +//In the type system so that they can be processed at compile-time for the +//automatic creation of floating-point filters. +// +//E.g. the arithmetic expression a * b + c can be represented by the type +// +//sum< product < argument<1>, argument<2> >, argument<3> > +// +//The arguments represent the input variables and are the leaves of the +//expression tree (other possible leaves are constants). Because they represent +//placeholders their indexing is one-based rather than zero-based, similar to +//the placeholders for std::bind. + +enum class operator_types { + sum, difference, product, abs, no_op, max, min +}; + +enum class operator_arities { nullary, unary, binary }; + +constexpr int sign_uncertain = -2; + +struct sum_error_type {}; +struct product_error_type {}; +struct no_error_type {}; + +template +struct internal_node +{ + static constexpr bool is_leaf = false; + static constexpr bool non_negative = false; + using all_children = boost::mp11::mp_list; //for convenience +}; + +template +struct internal_binary_node : internal_node +{ + using left = Left; + using right = Right; + static constexpr operator_arities operator_arity = + operator_arities::binary; +}; + +template +struct internal_unary_node : internal_node +{ + using child = Child; + static constexpr operator_arities operator_arity = operator_arities::unary; +}; + +template +struct sum : public internal_binary_node +{ + static constexpr bool sign_exact = + (Left::is_leaf && Right::is_leaf) + || ( Left::sign_exact && Right::sign_exact + && Left::non_negative && Right::non_negative); + static constexpr bool non_negative = Left::non_negative + && Right::non_negative; + static constexpr operator_types operator_type = operator_types::sum; + using error_type = sum_error_type; +}; + +template +struct difference : public internal_binary_node +{ + static constexpr bool sign_exact = Left::is_leaf && Right::is_leaf; + static constexpr bool non_negative = false; + static constexpr operator_types operator_type = operator_types::difference; + using error_type = sum_error_type; +}; + +template +struct product : public internal_binary_node +{ + static constexpr bool sign_exact = Left::sign_exact && Right::sign_exact; + static constexpr bool non_negative = + (Left::non_negative && Right::non_negative) + || std::is_same::value; + static constexpr operator_types operator_type = operator_types::product; + using error_type = product_error_type; +}; + +template +struct max : public internal_binary_node +{ + static constexpr bool sign_exact = Left::sign_exact && Right::sign_exact; + static constexpr bool non_negative = + Left::non_negative || Right::non_negative; + static constexpr operator_types operator_type = operator_types::max; + using error_type = no_error_type; +}; + +template +struct min : public internal_binary_node +{ + static constexpr bool sign_exact = Left::sign_exact && Right::sign_exact; + static constexpr bool non_negative = + Left::non_negative && Right::non_negative; + static constexpr operator_types operator_type = operator_types::min; + using error_type = no_error_type; +}; + +template +struct abs : public internal_unary_node +{ + using error_type = no_error_type; + static constexpr operator_types operator_type = operator_types::abs; + static constexpr bool sign_exact = Child::sign_exact; + static constexpr bool non_negative = true; +}; + +struct leaf +{ + static constexpr bool is_leaf = true; + static constexpr bool sign_exact = true; + static constexpr bool non_negative = false; + static constexpr operator_types operator_type = operator_types::no_op; + static constexpr operator_arities operator_arity = operator_arities::nullary; +}; + +template +struct argument : public leaf +{ + static constexpr std::size_t argn = Argn; +}; + +template +struct static_constant_interface : public leaf +{ + using value_type = NumberType; + static constexpr NumberType value = 0; //override + static constexpr std::size_t argn = 0; +}; + +template +using is_leaf = boost::mp11::mp_bool; + +//post_order_impl and post_order are templates for compile-time traversion of +//expression trees. The post_order traversal was chosen because it guarantees +//that children (subexpressions) are evaluated before there parents as is +//necessary for the evaluation of the arithmetic expressions that these +//expression trees represent. + +template +< + typename In, + typename Out, + template class Anchor = is_leaf, + bool IsBinary = In::operator_arity == operator_arities::binary, + bool AtAnchor = Anchor::value +> +struct post_order_impl; + +template +< + typename In, + typename Out, + template class Anchor, + bool IsBinary +> +struct post_order_impl +{ + using type = Out; +}; + +template class Anchor> +struct post_order_impl +{ + using leftl = typename post_order_impl + < + typename In::left, + boost::mp11::mp_list<>, + Anchor + >::type; + using rightl = typename post_order_impl + < + typename In::right, + boost::mp11::mp_list<>, + Anchor + >::type; + using merged = boost::mp11::mp_append; + using type = + boost::mp11::mp_unique>; +}; + +template class Anchor> +struct post_order_impl +{ + using childl = typename post_order_impl + < + typename In::child, + boost::mp11::mp_list<>, + Anchor + >::type; + using merged = boost::mp11::mp_append; + using type = + boost::mp11::mp_unique>; +}; + +template class Anchor = is_leaf> +using post_order = + typename post_order_impl, Anchor>::type; + +template +constexpr std::size_t max_argn = 0; + +template +constexpr std::size_t max_argn = + std::max(max_argn, + max_argn); + +template +constexpr std::size_t max_argn = + max_argn; + +template +constexpr std::size_t max_argn = + Expression::argn; + +using _1 = argument<1>; +using _2 = argument<2>; +using _3 = argument<3>; +using _4 = argument<4>; +using _5 = argument<5>; +using _6 = argument<6>; +using _7 = argument<7>; +using _8 = argument<8>; +using _9 = argument<9>; +using _10 = argument<10>; +using _11 = argument<11>; +using _12 = argument<12>; +using _13 = argument<13>; +using _14 = argument<14>; +using _15 = argument<15>; + +}} // namespace detail::generic_robust_predicates + +}} // namespace boost::geometry + +#endif // BOOST_GEOMETRY_EXTENSIONS_GENERIC_ROBUST_PREDICATES_STRATEGIES_CARTESIAN_DETAIL_EXPRESSION_TREE_HPP diff --git a/include/boost/geometry/extensions/generic_robust_predicates/strategies/cartesian/detail/expressions.hpp b/include/boost/geometry/extensions/generic_robust_predicates/strategies/cartesian/detail/expressions.hpp new file mode 100644 index 0000000000..b4b60cb37e --- /dev/null +++ b/include/boost/geometry/extensions/generic_robust_predicates/strategies/cartesian/detail/expressions.hpp @@ -0,0 +1,507 @@ +// Boost.Geometry (aka GGL, Generic Geometry Library) + +// Copyright (c) 2020 Tinko Bartels, Berlin, Germany. + +// Contributed and/or modified by Tinko Bartels, +// as part of Google Summer of Code 2020 program. + +// Use, modification and distribution is subject to 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) + +#ifndef BOOST_GEOMETRY_EXTENSIONS_GENERIC_ROBUST_PREDICATES_STRATEGIES_CARTESIAN_DETAIL_EXPRESSIONS_HPP +#define BOOST_GEOMETRY_EXTENSIONS_GENERIC_ROBUST_PREDICATES_STRATEGIES_CARTESIAN_DETAIL_EXPRESSIONS_HPP + +#include +#include +#include +#include + +#include + +namespace boost { namespace geometry +{ + +namespace detail { namespace generic_robust_predicates +{ + +static constexpr std::size_t sizet_sqrt(std::size_t n) +{ + std::size_t out = 0; + while (out * out < n) + { + ++out; + } + return out; +} + +template +struct indexed_exp +{ + using index = I; + using exp = Exp; +}; + +template +using extract_exp = typename IExp::exp; + +template +using extract_exps = boost::mp11::mp_transform; + +template +struct remove_from_minor_q +{ + template + using fn = boost::mp11::mp_bool + < + IExp::index::value / N::value == I::value + || IExp::index::value % N::value == 0 + >; +}; + +template +struct minor_q +{ + template + using fn = boost::mp11::mp_remove_if_q>; +}; + +template +struct minor_factor_q +{ + template + using fn = typename boost::mp11::mp_at_c::exp; +}; + +template +struct alternating_sum_impl {}; + +template +struct alternating_sum_impl +{ + using type = sum + < + difference, + typename alternating_sum_impl::type + >; +}; + +template +struct alternating_sum_impl +{ + using type = difference; +}; + +template +struct alternating_sum_impl +{ + using type = Exp; +}; + +template +using alternating_sum = typename alternating_sum_impl::type; + +template +struct det_impl; + +template +using make_det = + typename boost::mp11::mp_apply::type; + +template +struct det_impl +{ +private: + static constexpr std::size_t N = sizet_sqrt(sizeof...(Exps)); + static_assert(N * N == sizeof...(Exps), + "Arguments must be a square matrix."); + using indices = boost::mp11::mp_iota_c; + using minor_indices = boost::mp11::mp_iota_c; + using indexed_exps = boost::mp11::mp_transform + < + indexed_exp, + indices, + boost::mp11::mp_list + >; + using indexed_minors = boost::mp11::mp_transform_q + < + minor_q, indexed_exps>, + minor_indices + >; + using minors = boost::mp11::mp_transform; + using minor_dets = boost::mp11::mp_transform; + using minor_factors = boost::mp11::mp_transform_q + < + minor_factor_q, indexed_exps>, + minor_indices + >; + using alt_summands = boost::mp11::mp_transform + < + product, + minor_factors, + minor_dets + >; +public: + using type = boost::mp11::mp_apply + < + alternating_sum, + alt_summands + >; +}; + +template +struct det_impl +{ + using type = Exp; +}; + +template +using det = typename det_impl::type; + +template +struct orient_entry +{ + template + using fn = difference + < + argument, + argument + >; +}; + +template +using collinear_entry = argument; + +template +using orient = boost::mp11::mp_apply + < + det, + boost::mp11::mp_transform_q + < + orient_entry, + boost::mp11::mp_iota_c< N * N > + > + >; + +template +using collinear = boost::mp11::mp_apply + < + det, + boost::mp11::mp_transform + < + collinear_entry, + boost::mp11::mp_iota_c + > + >; + +template +struct incircle_entry; + +template +struct incircle_entry_offset +{ + template + using fn = typename incircle_entry::template fn + < + boost::mp11::mp_size_t + >; +}; + +template +struct multi_sum_impl; + +template +struct multi_sum_impl +{ + using type = sum::type>; +}; + +template +struct multi_sum_impl +{ + using type = Exp; +}; + +template <> +struct multi_sum_impl<> +{ + using type = void; +}; + +template +struct size_t_offset_q +{ + template + using fn = boost::mp11::mp_size_t; +}; + +template +struct incircle_entry +{ +private: + template + using diff = difference + < + argument + < + N * ((I::value) / (N + 1)) + + (I::value) % (N + 1) + 1 + >, + argument + < + N * (N + 1) + (I::value) % (N + 1) + 1 + > + >; + template + using lift_indices = boost::mp11::mp_transform_q + < + size_t_offset_q<(N+ 1) * (I / (N + 1))>, + boost::mp11::mp_iota_c + >; + + template + using line_diffs = + boost::mp11::mp_transform + < + diff, + lift_indices + >; + + template + using line_squares = + boost::mp11::mp_transform + < + product, + line_diffs, + line_diffs + >; + + template + using lift = typename boost::mp11::mp_apply + < + multi_sum_impl, + line_squares + >::type; +public: + template + using fn = + boost::mp11::mp_if_c + < + (I::value % (N + 1) != 0), + diff>, + lift + >; +}; + +template +using innsphere = boost::mp11::mp_apply + < + det, + boost::mp11::mp_transform_q + < + incircle_entry, + boost::mp11::mp_iota_c< (N + 1) * (N + 1) > + > + >; + +using orient2d = orient<2>; +using orient3d = orient<3>; +using incircle = innsphere<2>; +using insphere = innsphere<3>; + +struct incircle_simplified_impl +{ +private: + using qpx = difference, argument<7>>; + using qpy = difference, argument<8>>; + using rpx = difference, argument<7>>; + using rpy = difference, argument<8>>; + using tpx = difference, argument<7>>; + using tpy = difference, argument<8>>; +public: + using type = det + < + difference + < + product, + product + >, + sum + < + product + < + tpx, + difference < argument<5>, argument<3> > + >, + product + < + tpy, + difference < argument<6>, argument<4> > + > + >, + difference + < + product, + product + >, + sum + < + product + < + rpx, + difference < argument<1>, argument<3> > + >, + product + < + rpy, + difference < argument<2>, argument<4> > + > + > + >; +}; + +using incircle_simplified = incircle_simplified_impl::type; + +struct orient2d_no_translation_impl +{ +private: + using adet = det + < + argument<3>, argument<4>, + argument<5>, argument<6> + >; + using bdet = det + < + argument<1>, argument<2>, + argument<5>, argument<6> + >; + using cdet = det + < + argument<1>, argument<2>, + argument<3>, argument<4> + >; +public: + using type = sum< difference, cdet >; +}; + +using orient2d_no_translation = orient2d_no_translation_impl::type; + +struct orient3d_no_translation_impl +{ +private: + using adet = det + < + argument<4>, argument<5>, argument<6>, + argument<7>, argument<8>, argument<9>, + argument<10>,argument<11>,argument<12> + >; + using bdet = det + < + argument<1>, argument<2>, argument<3>, + argument<7>, argument<8>, argument<9>, + argument<10>,argument<11>,argument<12> + >; + using cdet = det + < + argument<1>, argument<2>, argument<3>, + argument<4>, argument<5>, argument<6>, + argument<10>,argument<11>,argument<12> + >; + using ddet = det + < + argument<1>, argument<2>, argument<3>, + argument<4>, argument<5>, argument<6>, + argument<7>, argument<8>, argument<9> + >; +public: + using type = sum< difference, difference >; +}; + +using orient3d_no_translation = orient3d_no_translation_impl::type; + +struct incircle_no_translation_impl +{ +private: + using alift = sum< product, argument<1>>, + product, argument<2>> >; + using blift = sum< product, argument<3>>, + product, argument<4>> >; + using clift = sum< product, argument<5>>, + product, argument<6>> >; + using dlift = sum< product, argument<7>>, + product, argument<8>> >; + using adet = det < argument<3>, argument<4>, blift, + argument<5>, argument<6>, clift, + argument<7>, argument<8>, dlift >; + using bdet = det < argument<1>, argument<2>, alift, + argument<5>, argument<6>, clift, + argument<7>, argument<8>, dlift >; + using cdet = det < argument<1>, argument<2>, alift, + argument<3>, argument<4>, blift, + argument<7>, argument<8>, dlift >; + using ddet = det < argument<1>, argument<2>, alift, + argument<3>, argument<4>, blift, + argument<5>, argument<6>, clift >; +public: + using type = sum < difference, difference >; +}; + +using incircle_no_translation = incircle_no_translation_impl::type; + +struct insphere_no_translation_impl +{ +private: + using alift = sum< sum, argument<1>>, + product, argument<2>>>, + product, argument<3>> >; + using blift = sum< sum, argument<4>>, + product, argument<5>>>, + product, argument<6>> >; + using clift = sum< sum, argument<7>>, + product, argument<8>>>, + product, argument<9>> >; + using dlift = sum< sum, argument<10>>, + product, argument<11>>>, + product, argument<12>> >; + using elift = sum< sum, argument<13>>, + product, argument<14>>>, + product, argument<15>> >; + using adet = det < argument<4>, argument<5>, argument<6>, blift, + argument<7>, argument<8>, argument<9>, clift, + argument<10>, argument<11>, argument<12>, dlift, + argument<13>, argument<14>, argument<15>, elift>; + using bdet = det < argument<1>, argument<2>, argument<3>, alift, + argument<7>, argument<8>, argument<9>, clift, + argument<10>, argument<11>, argument<12>, dlift, + argument<13>, argument<14>, argument<15>, elift>; + using cdet = det < argument<1>, argument<2>, argument<3>, alift, + argument<4>, argument<5>, argument<6>, blift, + argument<10>, argument<11>, argument<12>, dlift, + argument<13>, argument<14>, argument<15>, elift>; + using ddet = det < argument<1>, argument<2>, argument<3>, alift, + argument<4>, argument<5>, argument<6>, blift, + argument<7>, argument<8>, argument<9>, clift, + argument<13>, argument<14>, argument<15>, elift>; + using edet = det < argument<1>, argument<2>, argument<3>, alift, + argument<4>, argument<5>, argument<6>, blift, + argument<7>, argument<8>, argument<9>, clift, + argument<10>, argument<11>, argument<12>, dlift>; +public: + using type = sum + < + difference, + sum + < + difference, + edet + > + >; +}; + +using insphere_no_translation = insphere_no_translation_impl::type; + +}} // namespace detail::generic_robust_predicates + +}} // namespace boost::geometry + +#endif // BOOST_GEOMETRY_EXTENSIONS_GENERIC_ROBUST_PREDICATES_STRATEGIES_CARTESIAN_DETAIL_EXPRESSIONS_HPP diff --git a/include/boost/geometry/extensions/generic_robust_predicates/strategies/cartesian/detail/interval_error_bound.hpp b/include/boost/geometry/extensions/generic_robust_predicates/strategies/cartesian/detail/interval_error_bound.hpp new file mode 100644 index 0000000000..5716017015 --- /dev/null +++ b/include/boost/geometry/extensions/generic_robust_predicates/strategies/cartesian/detail/interval_error_bound.hpp @@ -0,0 +1,260 @@ +// Boost.Geometry (aka GGL, Generic Geometry Library) + +// Copyright (c) 2020 Tinko Bartels, Berlin, Germany. + +// Contributed and/or modified by Tinko Bartels, +// as part of Google Summer of Code 2020 program. + +// Use, modification and distribution is subject to 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) + +#ifndef BOOST_GEOMETRY_EXTENSIONS_GENERIC_ROBUST_PREDICATES_STRATEGIES_CARTESIAN_DETAIL_INTERVAL_ERROR_BOUND_HPP +#define BOOST_GEOMETRY_EXTENSIONS_GENERIC_ROBUST_PREDICATES_STRATEGIES_CARTESIAN_DETAIL_INTERVAL_ERROR_BOUND_HPP + +#include + +namespace boost { namespace geometry +{ + +namespace detail { namespace generic_robust_predicates +{ + +//This file contains methods to manipulate expression trees based on the ideas +//found in https://en.wikipedia.org/wiki/Interval_arithmetic. The purpose of +//this is to transform error expressions for semi-static filters, i.e. error +//expressions that compute error bounds for specific inputs, to error +//expressions for static filters that compute error bounds for upper and lower +//bounds on input values. + +template +struct interval_min_impl +{ + using type = Expression; +}; + +template +struct interval_max_impl +{ + using type = Expression; +}; + +template +struct interval_min_impl, MaxArgn> +{ +private: + using min_left = typename interval_min_impl::type; + using max_right = typename interval_max_impl::type; +public: + using type = difference; +}; + +template +struct interval_min_impl, MaxArgn> +{ +private: + using min_left = typename interval_min_impl::type; + using min_right = typename interval_min_impl::type; +public: + using type = sum; +}; + +template +struct interval_max_impl, MaxArgn> +{ +private: + using max_left = typename interval_max_impl::type; + using min_right = typename interval_min_impl::type; +public: + using type = difference; +}; + +template +struct interval_max_impl, MaxArgn> +{ +private: + using max_left = typename interval_max_impl::type; + using max_right = typename interval_max_impl::type; +public: + using type = sum; +}; + +template +struct interval_min_impl, MaxArgn> +{ +private: + using min_left = typename interval_min_impl::type; + using max_left = typename interval_max_impl::type; + using min_right = typename interval_min_impl::type; + using max_right = typename interval_max_impl::type; +public: + using type = min + < + min + < + product, + product + >, + min + < + product, + product + > + >; +}; + +template +struct interval_min_impl, MaxArgn> +{ +private: + using min_child = typename interval_min_impl::type; + using max_child = typename interval_max_impl::type; +public: + using type = min + < + product, + product + >; +}; + +template +struct interval_max_impl, MaxArgn> +{ +private: + using min_left = typename interval_min_impl::type; + using max_left = typename interval_max_impl::type; + using min_right = typename interval_min_impl::type; + using max_right = typename interval_max_impl::type; +public: + using type = max + < + max + < + product, + product + >, + max + < + product, + product + > + >; +}; + +template +struct interval_max_impl, MaxArgn> +{ +private: + using min_child = typename interval_min_impl::type; + using max_child = typename interval_max_impl::type; +public: + using type = max + < + product, + product + >; +}; + +template +struct interval_min_impl, MaxArgn> +{ +private: + using min_child = typename interval_min_impl::type; + using max_child = typename interval_max_impl::type; +public: + using type = min, abs>; +}; + +template +struct interval_max_impl, MaxArgn> +{ +private: + using min_child = typename interval_min_impl::type; + using max_child = typename interval_max_impl::type; +public: + using type = max, abs>; +}; + +template +struct interval_min_impl, MaxArgn> +{ + using type = argument; +}; + +template +struct interval_max_impl, MaxArgn> +{ + using type = argument; +}; + +template +struct interval_impl +{ + using type = Expression; +}; + +template +struct interval_impl, MaxArgn> +{ +private: + using min_child = typename interval_min_impl::type; + using max_child = typename interval_max_impl::type; +public: + using type = max, abs>; +}; + +template +struct interval_impl, MaxArgn> +{ +private: + using left = typename interval_impl::type; + using right = typename interval_impl::type; +public: + using type = max; +}; + +template +struct interval_impl, MaxArgn> +{ +private: + using left = typename interval_impl::type; + using right = typename interval_impl::type; +public: + using type = difference; +}; + +template +struct interval_impl, MaxArgn> +{ +private: + using left = typename interval_impl::type; + using right = typename interval_impl::type; +public: + using type = sum; +}; + +template +struct interval_impl, MaxArgn> +{ +private: + using left = typename interval_impl::type; + using right = typename interval_impl::type; +public: + using type = product; +}; + +template +using interval_min = typename interval_min_impl::type; + +template +using interval_max = typename interval_max_impl::type; + +template +using interval = + typename interval_impl>::type; + +}} // namespace detail::generic_robust_predicates + +}} // namespace boost::geometry + +#endif // BOOST_GEOMETRY_EXTENSIONS_GENERIC_ROBUST_PREDICATES_STRATEGIES_CARTESIAN_DETAIL_INTERVAL_ERROR_BOUND_HPP diff --git a/include/boost/geometry/extensions/generic_robust_predicates/strategies/cartesian/detail/semi_static_filter.hpp b/include/boost/geometry/extensions/generic_robust_predicates/strategies/cartesian/detail/semi_static_filter.hpp new file mode 100644 index 0000000000..e91d3cc477 --- /dev/null +++ b/include/boost/geometry/extensions/generic_robust_predicates/strategies/cartesian/detail/semi_static_filter.hpp @@ -0,0 +1,120 @@ +// Boost.Geometry (aka GGL, Generic Geometry Library) + +// Copyright (c) 2020 Tinko Bartels, Berlin, Germany. + +// Contributed and/or modified by Tinko Bartels, +// as part of Google Summer of Code 2020 program. + +// Use, modification and distribution is subject to 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) + +#ifndef BOOST_GEOMETRY_EXTENSIONS_GENERIC_ROBUST_PREDICATES_STRATEGIES_CARTESIAN_DETAIL_SEMI_STATIC_FILTER_HPP +#define BOOST_GEOMETRY_EXTENSIONS_GENERIC_ROBUST_PREDICATES_STRATEGIES_CARTESIAN_DETAIL_SEMI_STATIC_FILTER_HPP + +#include + +#include +#include +#include + +#include +#include + +namespace boost { namespace geometry +{ + +namespace detail { namespace generic_robust_predicates +{ + +//The semi static filter evaluates an Expression and an ErrorExpression. On +//application it evaluates both expressions and verifies the sign of the +//Expression by the following rules: +//If the absolute value of the number resulting from the evaluation of the +//expression is larger than or equal to the result for the error expression, +//then the sign can be returned. If not, then a constant is returned to +//represent an uncertain sign. +// +//The filters that are build from this template are meant to be semi-static in +//the sense that the error expression including constants that depend only the +//epsilon of the calculation type are known statically at compile-time but the +//final value of the error bound depends on the specific inputs for each call +//to the filter and parts of it need to be computed at each call. +// +//This filter is stateless. + +template +< + typename Expression, + typename CalculationType, + typename ErrorExpression +> +struct semi_static_filter +{ +private: + using stack = typename boost::mp11::mp_unique>; + using evals = typename boost::mp11::mp_remove_if; + using error_eval_stack = boost::mp11::mp_unique + < + post_order + >; + using error_eval_stack_remainder = boost::mp11::mp_set_difference + < + error_eval_stack, + evals + >; + using all_evals = boost::mp11::mp_append + < + evals, + error_eval_stack_remainder + >; + using ct = CalculationType; +public: + static constexpr bool stateful = false; + static constexpr bool updates = false; + + template + static inline ct error_bound(CTs const&... args) + { + std::array input + {{ static_cast(args)... }}; + return evaluate_expression(input); + } + + template + static inline int apply(CTs const&... args) + { + std::array input + {{ static_cast(args)... }}; + std::array::value> results; + evaluate_expressions(input, results, all_evals{}); + constexpr std::size_t i_eb = + boost::mp11::mp_find::value; + const ct error_bound = results[i_eb]; + constexpr std::size_t i_e = + boost::mp11::mp_find::value; + const ct det = results[i_e]; + if (det > error_bound) + { + return 1; + } + else if (det < -error_bound) + { + return -1; + } + else if (error_bound == 0 && det == 0) + { + return 0; + } + else + { + return sign_uncertain; + } + } +}; + +}} // namespace detail::generic_robust_predicates + +}} // namespace boost::geometry + +#endif // BOOST_GEOMETRY_EXTENSIONS_GENERIC_ROBUST_PREDICATES_STRATEGIES_CARTESIAN_DETAIL_SEMI_STATIC_FILTER_HPP diff --git a/include/boost/geometry/extensions/generic_robust_predicates/strategies/cartesian/detail/stage_a.hpp b/include/boost/geometry/extensions/generic_robust_predicates/strategies/cartesian/detail/stage_a.hpp new file mode 100644 index 0000000000..e1e50ad501 --- /dev/null +++ b/include/boost/geometry/extensions/generic_robust_predicates/strategies/cartesian/detail/stage_a.hpp @@ -0,0 +1,95 @@ +// Boost.Geometry (aka GGL, Generic Geometry Library) + +// Copyright (c) 2020 Tinko Bartels, Berlin, Germany. + +// Contributed and/or modified by Tinko Bartels, +// as part of Google Summer of Code 2020 program. + +// Use, modification and distribution is subject to 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) + +#ifndef BOOST_GEOMETRY_EXTENSIONS_GENERIC_ROBUST_PREDICATES_STRATEGIES_CARTESIAN_DETAIL_STAGE_A_HPP +#define BOOST_GEOMETRY_EXTENSIONS_GENERIC_ROBUST_PREDICATES_STRATEGIES_CARTESIAN_DETAIL_STAGE_A_HPP + +#include +#include +#include + +#include +#include +#include + +namespace boost { namespace geometry +{ + +namespace detail { namespace generic_robust_predicates +{ + +//The following templates apply the ideas implemented in error_bound.hpp and +//interval_error_bound.hpp to generate error expressions from expressions. + +template +< + typename Expression, + typename CalculationType +> +using stage_a_error_bound_expression = + typename stage_a_error_bound_expression_impl + < + Expression, + CalculationType + >::type; + +template +< + typename Expression, + typename CalculationType +> +using stage_a_semi_static = semi_static_filter + < + Expression, + CalculationType, + stage_a_error_bound_expression + < + Expression, + CalculationType + > + >; + +template +< + typename Expression, + typename CalculationType +> +using stage_a_static = static_filter + < + Expression, + CalculationType, + interval + < + stage_a_error_bound_expression + < + Expression, + CalculationType + > + > + >; + +template +< + typename Expression, + typename CalculationType +> +using stage_a_almost_static = almost_static_filter + < + Expression, + CalculationType, + stage_a_static + >; + +}} // namespace detail::generic_robust_predicates + +}} // namespace boost::geometry + +#endif // BOOST_GEOMETRY_EXTENSIONS_GENERIC_ROBUST_PREDICATES_STRATEGIES_CARTESIAN_DETAIL_STAGE_A_HPP diff --git a/include/boost/geometry/extensions/generic_robust_predicates/strategies/cartesian/detail/stage_a_error_bound.hpp b/include/boost/geometry/extensions/generic_robust_predicates/strategies/cartesian/detail/stage_a_error_bound.hpp new file mode 100644 index 0000000000..1f404f22eb --- /dev/null +++ b/include/boost/geometry/extensions/generic_robust_predicates/strategies/cartesian/detail/stage_a_error_bound.hpp @@ -0,0 +1,235 @@ +#ifndef BOOST_GEOMETRY_EXTENSIONS_GENERIC_ROBUST_PREDICATES_STRATEGIES_CARTESIAN_DETAIL_STAGE_A_ERROR_BOUND_HPP +#define BOOST_GEOMETRY_EXTENSIONS_GENERIC_ROBUST_PREDICATES_STRATEGIES_CARTESIAN_DETAIL_STAGE_A_ERROR_BOUND_HPP + +#include +#include + +#include + +#include + +namespace boost { namespace geometry +{ + +namespace detail { namespace generic_robust_predicates +{ + +enum class stage_a_error_propagation_cases { + exact, op_on_exacts, sum_or_diff, product +}; + +template +< + typename Expression, + operator_arities Arity = Expression::operator_arity +> +constexpr stage_a_error_propagation_cases stage_a_error_propagation_case = + stage_a_error_propagation_cases::exact; + +template +constexpr stage_a_error_propagation_cases stage_a_error_propagation_case + < + Expression, + operator_arities::binary + > + = Expression::left::is_leaf && Expression::right::is_leaf ? + stage_a_error_propagation_cases::op_on_exacts : + (Expression::operator_type == operator_types::product ? + stage_a_error_propagation_cases::product : + stage_a_error_propagation_cases::sum_or_diff); + + +template +< + typename Expression, + stage_a_error_propagation_cases = + stage_a_error_propagation_case +> +struct stage_a_error_bound {}; + +template +struct stage_a_error_bound +{ + using magnitude = abs; + static constexpr std::array a {0, 0, 0}; +}; + +template +struct stage_a_error_bound + < + Expression, + stage_a_error_propagation_cases::op_on_exacts + > +{ + using magnitude = abs; + static constexpr std::array a {1, 0, 0}; +}; + +constexpr std::array coeff_max(const std::array a, + const std::array b) +{ + bool a_bigger = + a[0] > b[0] + || (a[0] == b[0] && a[1] > b[1]) + || (a[0] == b[0] && a[1] == b[1] && a[2] > b[2]); + return a_bigger ? a : b; +} + +constexpr std::array coeff_product(const std::array a, + const std::array b) +{ + return std::array { + a[0] + b[0], + a[1] + b[1] + a[0] * b[0], + a[2] + b[2] + a[0] * b[1] + a[1] * b[0] + }; +} + +constexpr std::array +coeff_mult_by_1_plus_eps(const std::array a) +{ + return std::array { + a[0], + a[1] + a[0], + a[2] + a[1] + }; +} + +constexpr std::array coeff_inc_first(const std::array a) +{ + return std::array { a[0] + 1, a[1], a[2] }; +} + +constexpr std::array +coeff_div_by_1_minus_eps(const std::array a) +{ + return std::array { + a[0], + a[0] + a[1], + a[0] + a[1] + a[2] + 1 + }; +} + +template +constexpr int round_to_next_2_pow() { + static_assert(N >= 1, "Expects positive integer."); + int out = 1; + while(out < N) + { + out *= 2; + } + return out; +} + +template +struct stage_a_error_bound + < + Expression, + stage_a_error_propagation_cases::sum_or_diff + > +{ +private: + using leb = stage_a_error_bound; + using reb = stage_a_error_bound; + static constexpr auto max_a = coeff_max(leb::a, reb::a); +public: + using magnitude = sum; + static constexpr std::array a = + coeff_inc_first(coeff_mult_by_1_plus_eps(max_a)); +}; + +template +struct stage_a_error_bound + < + Expression, + stage_a_error_propagation_cases::product + > +{ +private: + using leb = stage_a_error_bound; + using reb = stage_a_error_bound; + static constexpr std::array a_prod = coeff_product(leb::a, reb::a); +public: + using magnitude = + product; + static constexpr std::array a = + coeff_inc_first(coeff_mult_by_1_plus_eps(a_prod)); +}; + +template +< + typename Expression, + stage_a_error_propagation_cases = + stage_a_error_propagation_case +> +struct stage_a_condition {}; + +template +struct stage_a_condition + < + Expression, + stage_a_error_propagation_cases::sum_or_diff + > +{ +private: + using leb = stage_a_error_bound; + using reb = stage_a_error_bound; + static constexpr auto max_a = coeff_max(leb::a, reb::a); + static constexpr std::array a = + coeff_mult_by_1_plus_eps(coeff_mult_by_1_plus_eps(coeff_div_by_1_minus_eps(max_a))); + static constexpr int c = round_to_next_2_pow(); + static constexpr int eps_square_coeff = + a[2] > 0 ? c * ((a[1] + 1) / c + 1) : c * (a[1] / c + 1); +public: + using magnitude = sum; + static constexpr std::array coefficients {a[0], eps_square_coeff}; +}; + +template +struct stage_a_condition + < + Expression, + stage_a_error_propagation_cases::product + > +{ +private: + using leb = stage_a_error_bound; + using reb = stage_a_error_bound; + static constexpr auto a_prod = coeff_product(leb::a, reb::a); + static constexpr std::array a = + coeff_mult_by_1_plus_eps(coeff_mult_by_1_plus_eps(coeff_div_by_1_minus_eps(a_prod))); + static constexpr int c = round_to_next_2_pow(); + static constexpr int eps_square_coeff = + a[2] > 0 ? c * ((a[1] + 1) / c + 1) : c * (a[1] / c + 1); +public: + using magnitude = product; + static constexpr std::array coefficients {a[0], eps_square_coeff}; +}; + +template +< + typename Expression, + typename CalculationType +> +struct stage_a_error_bound_expression_impl +{ +private: + using ct = CalculationType; + using stage_a_cond = stage_a_condition; + static constexpr ct eps = std::numeric_limits::epsilon() / 2.0; + struct constant : public static_constant_interface + { + static constexpr ct value = + stage_a_cond::coefficients[0] * eps + + stage_a_cond::coefficients[1] * eps * eps; + static constexpr bool non_negative = true; + }; +public: + using type = product; +}; + +}} // namespace detail::generic_robust_predicates + +}} // namespace boost::geometry + +#endif // BOOST_GEOMETRY_EXTENSIONS_GENERIC_ROBUST_PREDICATES_STRATEGIES_CARTESIAN_DETAIL_STAGE_A_ERROR_BOUND_HPP diff --git a/include/boost/geometry/extensions/generic_robust_predicates/strategies/cartesian/detail/stage_b.hpp b/include/boost/geometry/extensions/generic_robust_predicates/strategies/cartesian/detail/stage_b.hpp new file mode 100644 index 0000000000..f1ecd505d4 --- /dev/null +++ b/include/boost/geometry/extensions/generic_robust_predicates/strategies/cartesian/detail/stage_b.hpp @@ -0,0 +1,253 @@ +// Boost.Geometry (aka GGL, Generic Geometry Library) + +// Copyright (c) 2020 Tinko Bartels, Berlin, Germany. + +// Contributed and/or modified by Tinko Bartels, +// as part of Google Summer of Code 2020 program. + +// Use, modification and distribution is subject to 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) + +#ifndef BOOST_GEOMETRY_EXTENSIONS_GENERIC_ROBUST_PREDICATES_STRATEGIES_CARTESIAN_DETAIL_STAGE_B_HPP +#define BOOST_GEOMETRY_EXTENSIONS_GENERIC_ROBUST_PREDICATES_STRATEGIES_CARTESIAN_DETAIL_STAGE_B_HPP + +#include +#include + +#include +#include + +#include + +namespace boost { namespace geometry +{ + +namespace detail { namespace generic_robust_predicates +{ + +template +< + typename Expression, + bool IsDifference = Expression::operator_type == operator_types::difference +> +struct is_leaf_difference_impl +{ + using type = boost::mp11::mp_false; +}; + +template +< + typename Expression +> +struct is_leaf_difference_impl +{ + using type = boost::mp11::mp_bool + < + Expression::left::is_leaf && Expression::right::is_leaf + >; +}; + +template +using is_leaf_difference = typename is_leaf_difference_impl::type; + +template +< + typename Evals, + typename LeafDifferences, + typename AccumulatedSizes, + typename Iter, + typename CT, + std::size_t RemainingDifferences = + boost::mp11::mp_size::value +> +struct all_differences_zero_tail +{ + template + static bool apply(Iter begin, Iter end, CTs const&... args) + { + using eval = boost::mp11::mp_front; + using left = typename eval::left; + using right = typename eval::right; + std::array input + {{ static_cast(args)... }}; + CT left_val = input[left::argn - 1]; + CT right_val = input[right::argn - 1]; + using eval_index = boost::mp11::mp_find; + constexpr std::size_t start = + boost::mp11::mp_at::value; + *(begin + start) = left_val - right_val; + if (two_difference_tail(left_val, right_val, *(begin + start)) + == CT(0)) + { + return all_differences_zero_tail + < + Evals, + boost::mp11::mp_pop_front, + AccumulatedSizes, + Iter, + CT + >::apply(begin, end, args...); + } + return false; + } +}; + +template +< + typename Evals, + typename LeafDifferences, + typename AccumulatedSizes, + typename Iter, + typename CT +> +struct all_differences_zero_tail + < + Evals, + LeafDifferences, + AccumulatedSizes, + Iter, + CT, + 0 + > +{ + template + static bool apply(Iter, Iter, CTs const&...) + { + return true; + } +}; + +template +< + typename Expression, + typename CalculationType, + template class ZEPolicy = default_zero_elimination_policy, + template class FEPolicy = default_fast_expansion_sum_policy +> +struct stage_b +{ +private: + using ct = CalculationType; +public: + static constexpr bool stateful = false; + static constexpr bool updates = false; + + template + static inline int apply(CTs const&... args) + { + using stack = typename boost::mp11::mp_unique>; + using evals = typename boost::mp11::mp_remove_if; + using sizes_pre = boost::mp11::mp_transform + < + expansion_size_stage_b, + boost::mp11::mp_pop_back + >; + using sizes = boost::mp11::mp_push_back + < + sizes_pre, + boost::mp11::mp_size_t + < + final_expansion_size + < + Expression, + expansion_size_stage_b + < + typename Expression::left + >::value, + expansion_size_stage_b + < + typename Expression::right + >::value + >() + > + >; + + using accumulated_sizes = boost::mp11::mp_push_front + < + boost::mp11::mp_partial_sum + < + sizes, + boost::mp11::mp_size_t<0>, + boost::mp11::mp_plus + >, + boost::mp11::mp_size_t<0> + >; + + using result_array = + std::array + < + ct, + boost::mp11::mp_back::value + >; + result_array results; + + using leaf_differences = + boost::mp11::mp_copy_if; + + bool all_zero = all_differences_zero_tail + < + evals, + leaf_differences, + accumulated_sizes, + typename result_array::iterator, + ct + >::apply(results.begin(), results.end(), args...); + + if ( !all_zero ) + { + return sign_uncertain; + } + + using remainder = + typename boost::mp11::mp_remove_if + < + evals, + is_leaf_difference + >; + + using ze_evals = boost::mp11::mp_copy_if_q + < + remainder, + is_zero_elim_q + >; + std::array + < + typename result_array::iterator, + boost::mp11::mp_size::value + > ze_ends; + + auto most_significant = eval_expansions_impl + < + evals, + remainder, + sizes, + accumulated_sizes, + ze_evals, + decltype(results.begin()), + ct, + ZEPolicy, + FEPolicy, + true + >::apply(results.begin(), results.end(), ze_ends, args...) - 1; + + if ( *most_significant == 0) + { + return 0; + } + else if ( *most_significant > 0 ) + { + return 1; + } + else + { + return -1; + } + } +}; + +}} // namespace detail::generic_robust_predicates + +}} // namespace boost::geometry + +#endif // BOOST_GEOMETRY_EXTENSIONS_GENERIC_ROBUST_PREDICATES_STRATEGIES_CARTESIAN_DETAIL_STAGE_B_HPP diff --git a/include/boost/geometry/extensions/generic_robust_predicates/strategies/cartesian/detail/stage_d.hpp b/include/boost/geometry/extensions/generic_robust_predicates/strategies/cartesian/detail/stage_d.hpp new file mode 100644 index 0000000000..0e7f488a46 --- /dev/null +++ b/include/boost/geometry/extensions/generic_robust_predicates/strategies/cartesian/detail/stage_d.hpp @@ -0,0 +1,116 @@ +// Boost.Geometry (aka GGL, Generic Geometry Library) + +// Copyright (c) 2020 Tinko Bartels, Berlin, Germany. + +// Contributed and/or modified by Tinko Bartels, +// as part of Google Summer of Code 2020 program. + +// Use, modification and distribution is subject to 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) + +#ifndef BOOST_GEOMETRY_EXTENSIONS_GENERIC_ROBUST_PREDICATES_STRATEGIES_CARTESIAN_DETAIL_STAGE_D_HPP +#define BOOST_GEOMETRY_EXTENSIONS_GENERIC_ROBUST_PREDICATES_STRATEGIES_CARTESIAN_DETAIL_STAGE_D_HPP + +#include +#include +#include + +#include +#include +#include +#include + +#include +#include + +namespace boost { namespace geometry +{ + +namespace detail { namespace generic_robust_predicates +{ + +template +< + typename Expression, + typename CalculationType, + template class ZEPolicy = default_zero_elimination_policy, + template class FEPolicy = default_fast_expansion_sum_policy +> +struct stage_d +{ +private: + using ct = CalculationType; +public: + static constexpr bool stateful = false; + static constexpr bool updates = false; + + template + static inline int apply(CTs const&... args) + { + using stack = typename boost::mp11::mp_unique>; + using evals = typename boost::mp11::mp_remove_if; + using sizes_pre = boost::mp11::mp_transform + < + expansion_size, + boost::mp11::mp_pop_back + >; + using sizes = boost::mp11::mp_push_back + < + sizes_pre, + boost::mp11::mp_size_t + < + final_expansion_size + < + Expression, + expansion_size::value, + expansion_size::value + >() + > + >; + using accumulated_sizes = boost::mp11::mp_push_front + < + boost::mp11::mp_partial_sum + < + sizes, + boost::mp11::mp_size_t<0>, + boost::mp11::mp_plus + >, + boost::mp11::mp_size_t<0> + >; + + using result_array = + std::array::value>; + result_array results; + + auto most_significant = eval_expansions + < + evals, + sizes, + accumulated_sizes, + decltype(results.begin()), + ct, + false, + ZEPolicy, + FEPolicy + >(results.begin(), results.end(), args...) - 1; + if ( *most_significant == 0) + { + return 0; + } + else if ( *most_significant > 0 ) + { + return 1; + } + else + { + return -1; + } + } +}; + +}} // namespace detail::generic_robust_predicates + +}} // namespace boost::geometry + +#endif // BOOST_GEOMETRY_EXTENSIONS_GENERIC_ROBUST_PREDICATES_STRATEGIES_CARTESIAN_DETAIL_STAGE_D_HPP diff --git a/include/boost/geometry/extensions/generic_robust_predicates/strategies/cartesian/detail/staged_predicate.hpp b/include/boost/geometry/extensions/generic_robust_predicates/strategies/cartesian/detail/staged_predicate.hpp new file mode 100644 index 0000000000..29c3dbd531 --- /dev/null +++ b/include/boost/geometry/extensions/generic_robust_predicates/strategies/cartesian/detail/staged_predicate.hpp @@ -0,0 +1,238 @@ +// Boost.Geometry (aka GGL, Generic Geometry Library) + +// Copyright (c) 2020 Tinko Bartels, Berlin, Germany. + +// Contributed and/or modified by Tinko Bartels, +// as part of Google Summer of Code 2020 program. + +// Use, modification and distribution is subject to 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) + +#ifndef BOOST_GEOMETRY_EXTENSIONS_GENERIC_ROBUST_PREDICATES_STRATEGIES_CARTESIAN_DETAIL_STAGED_PREDICATE_HPP +#define BOOST_GEOMETRY_EXTENSIONS_GENERIC_ROBUST_PREDICATES_STRATEGIES_CARTESIAN_DETAIL_STAGED_PREDICATE_HPP + +#include +#include +#include + +#include +#include +#include + +#include + +namespace boost { namespace geometry +{ + +namespace detail { namespace generic_robust_predicates +{ + + +template +constexpr bool has_next_and_is_stateful = + boost::mp11::mp_front::stateful; + +template <> +constexpr bool has_next_and_is_stateful> = false; + +template +< + typename StatefulStages, + typename RemainingStages, + bool next_is_stateful = has_next_and_is_stateful +> +struct next_stage +{ + template + static inline int apply(StatefulStages const& stages, + CTs const&... args) + { + using stage = boost::mp11::mp_front; + int sign = stage::template apply<>(args...); + if (sign == sign_uncertain) + { + return next_stage + < + StatefulStages, + boost::mp11::mp_pop_front + >::apply(stages, args...); + } + else + { + return sign; + } + } +}; + +template +< + typename StatefulStages, + typename RemainingStages +> +struct next_stage + < + StatefulStages, + RemainingStages, + true + > +{ + template + static inline int apply(StatefulStages const& stages, + CTs const&... args) + { + using stage = boost::mp11::mp_front; + int sign = std::get(stages).apply(args...); + if (sign == sign_uncertain) + { + return next_stage + < + StatefulStages, + boost::mp11::mp_pop_front + >::apply(stages, args...); + } + else + { + return sign; + } + } +}; + +template +< + typename StatefulStages +> +struct next_stage + < + StatefulStages, + boost::mp11::mp_list<>, + false + > +{ + template + static inline int apply(StatefulStages const&, + CTs const&...) + { + return sign_uncertain; + } +}; + +template using is_stateful = boost::mp11::mp_bool; + +template +struct construct_stage_impl +{ + template static inline Stage apply(Array const& a) + { + Stage out(a); + return out; + } +}; + +template +struct construct_stage_impl +{ + template static inline Stage apply(Array const&) + { + Stage out; + return out; + } +}; + +template +struct construct_stages_impl +{ + template + static inline boost::mp11::mp_rename + apply(Array const& a) + { + using stage = boost::mp11::mp_front; + std::tuple first(construct_stage_impl::apply(a)); + return std::tuple_cat(first, + construct_stages_impl + < + boost::mp11::mp_pop_front + >::apply(a)); + } +}; + +template <> +struct construct_stages_impl> +{ + template + static inline std::tuple<> apply(Array const&) + { + return std::tuple<>{}; + } +}; + +template +using is_updatable = boost::mp11::mp_bool; + +template +struct update_all_impl +{ + template + static void apply(Stages& stages, CTs const&... args) + { + using current_stage = boost::mp11::mp_front; + std::get(stages).update(args...); + update_all_impl> + ::apply(stages, args...); + } +}; + +template <> +struct update_all_impl> +{ + template + static void apply(Stages&, CTs const&...) {} +}; + +template +< + typename CalculationType, + typename ...Stages +> +struct staged_predicate +{ +private: + using ct = CalculationType; + using stages = std::tuple; + using stateful_stages = boost::mp11::mp_copy_if; + using updatable_stages = boost::mp11::mp_copy_if; + stateful_stages m_stages; + using all_stages = boost::mp11::mp_list; +public: + static constexpr bool stateful = + boost::mp11::mp_any_of::value; + static constexpr bool updates = + boost::mp11::mp_any_of::value; + template + inline staged_predicate(CTs const&... args) : m_stages( + construct_stages_impl::apply( + std::array{ static_cast(args)... } + )) {} + + template + inline void update(CTs const&... args) + { + update_all_impl::apply(m_stages, args...); + } + + template + inline int apply(CTs const&... args) const + { + return next_stage + < + stateful_stages, + all_stages + >::apply(m_stages, args...); + } +}; + +}} // namespace detail::generic_robust_predicates + +}} // namespace boost::geometry + +#endif // BOOST_GEOMETRY_EXTENSIONS_GENERIC_ROBUST_PREDICATES_STRATEGIES_CARTESIAN_DETAIL_STAGED_PREDICATE_HPP diff --git a/include/boost/geometry/extensions/generic_robust_predicates/strategies/cartesian/detail/static_filter.hpp b/include/boost/geometry/extensions/generic_robust_predicates/strategies/cartesian/detail/static_filter.hpp new file mode 100644 index 0000000000..dce3f2ba58 --- /dev/null +++ b/include/boost/geometry/extensions/generic_robust_predicates/strategies/cartesian/detail/static_filter.hpp @@ -0,0 +1,109 @@ +// Boost.Geometry (aka GGL, Generic Geometry Library) + +// Copyright (c) 2020 Tinko Bartels, Berlin, Germany. + +// Contributed and/or modified by Tinko Bartels, +// as part of Google Summer of Code 2020 program. + +// Use, modification and distribution is subject to 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) + +#ifndef BOOST_GEOMETRY_EXTENSIONS_GENERIC_ROBUST_PREDICATES_STRATEGIES_CARTESIAN_DETAIL_STATIC_FILTER_HPP +#define BOOST_GEOMETRY_EXTENSIONS_GENERIC_ROBUST_PREDICATES_STRATEGIES_CARTESIAN_DETAIL_STATIC_FILTER_HPP + +#include +#include + +#include +#include + +#include +#include + +namespace boost { namespace geometry +{ + +namespace detail { namespace generic_robust_predicates +{ + +//The static filter works similar to the semi static filter with the exception +//that the error bound is only computed once at construction time. For this +//purpose, the filter is stateful and it requires upper and lower bounds on the +//inputs at compile time like this: +// +//static_filter<...> f(max1, max2, ..., min1, min2, ...); +// +//for each call to the filter, i.e. f.apply(arg1, arg2, ...); it must hold that +//min1 <= arg1 <= max1 and so on. + +template +< + typename Expression, + typename CalculationType, + typename ErrorExpression +> +class static_filter +{ +private: + using ct = CalculationType; + ct m_error_bound; + static constexpr std::size_t extrema_count = max_argn; +public: + static constexpr bool stateful = true; + static constexpr bool updates = false; + + inline static_filter() : m_error_bound(std::numeric_limits::max()) {} + + inline ct error_bound() const { return m_error_bound; } + + template + inline static_filter(CTs const&... args) + : m_error_bound(evaluate_expression( + std::array + {static_cast(args)...})) + { + static_assert(sizeof...(CTs) == extrema_count, + "Number of constructor arguments is incompatible with error expression."); + } + + inline static_filter(std::array const& extrema) + : m_error_bound(evaluate_expression(extrema)) {} + + inline static_filter(ct const& b) : m_error_bound(b) {} + + template + inline int apply(CTs const&... args) const + { + std::array input {static_cast(args)...}; + const ct det = evaluate_expression(input); + if (det > m_error_bound) + { + return 1; + } + else if (det < -m_error_bound) + { + return -1; + } + else if (m_error_bound == 0 && det == 0) + { + return 0; + } + else + { + return sign_uncertain; + } + } + + template + inline int operator()(CTs const&... args) const + { + return apply(args...); + } +}; + +}} // namespace detail::generic_robust_predicates + +}} // namespace boost::geometry + +#endif // BOOST_GEOMETRY_EXTENSIONS_GENERIC_ROBUST_PREDICATES_STRATEGIES_CARTESIAN_DETAIL_STATIC_FILTER_HPP