From c1c75b57c5848b5d3d1171d78a149f8d3572abad Mon Sep 17 00:00:00 2001 From: Tinko Bartels Date: Tue, 28 Jul 2020 13:35:26 +0200 Subject: [PATCH 1/6] [extensions] Automatically generated floating point filters. --- extensions/example/Jamfile | 1 + .../example/generic_robust_predicates/Jamfile | 12 + .../static_side_2d.cpp | 79 +++ extensions/test/Jamfile | 1 + .../test/generic_robust_predicates/Jamfile | 12 + .../generic_robust_predicates/approximate.cpp | 66 ++ .../test/generic_robust_predicates/side3d.cpp | 76 +++ .../cartesian/detail/almost_static_filter.hpp | 133 ++++ .../cartesian/detail/approximate.hpp | 478 ++++++++++++++ .../cartesian/detail/coefficient_list.hpp | 471 ++++++++++++++ .../cartesian/detail/error_bound.hpp | 605 ++++++++++++++++++ .../cartesian/detail/expression_tree.hpp | 250 ++++++++ .../cartesian/detail/expressions.hpp | 103 +++ .../cartesian/detail/interval_error_bound.hpp | 243 +++++++ .../cartesian/detail/semi_static_filter.hpp | 91 +++ .../strategies/cartesian/detail/stage_a.hpp | 141 ++++ .../cartesian/detail/static_filter.hpp | 98 +++ .../cartesian/detail/static_util.hpp | 130 ++++ 18 files changed, 2990 insertions(+) create mode 100644 extensions/example/generic_robust_predicates/Jamfile create mode 100644 extensions/example/generic_robust_predicates/static_side_2d.cpp create mode 100644 extensions/test/generic_robust_predicates/Jamfile create mode 100644 extensions/test/generic_robust_predicates/approximate.cpp create mode 100644 extensions/test/generic_robust_predicates/side3d.cpp create mode 100644 include/boost/geometry/extensions/generic_robust_predicates/strategies/cartesian/detail/almost_static_filter.hpp create mode 100644 include/boost/geometry/extensions/generic_robust_predicates/strategies/cartesian/detail/approximate.hpp create mode 100644 include/boost/geometry/extensions/generic_robust_predicates/strategies/cartesian/detail/coefficient_list.hpp create mode 100644 include/boost/geometry/extensions/generic_robust_predicates/strategies/cartesian/detail/error_bound.hpp create mode 100644 include/boost/geometry/extensions/generic_robust_predicates/strategies/cartesian/detail/expression_tree.hpp create mode 100644 include/boost/geometry/extensions/generic_robust_predicates/strategies/cartesian/detail/expressions.hpp create mode 100644 include/boost/geometry/extensions/generic_robust_predicates/strategies/cartesian/detail/interval_error_bound.hpp create mode 100644 include/boost/geometry/extensions/generic_robust_predicates/strategies/cartesian/detail/semi_static_filter.hpp create mode 100644 include/boost/geometry/extensions/generic_robust_predicates/strategies/cartesian/detail/stage_a.hpp create mode 100644 include/boost/geometry/extensions/generic_robust_predicates/strategies/cartesian/detail/static_filter.hpp create mode 100644 include/boost/geometry/extensions/generic_robust_predicates/strategies/cartesian/detail/static_util.hpp 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..c18ef3c4f1 --- /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(int argc, char** argv) +{ + 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..0595681dce --- /dev/null +++ b/extensions/test/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) + +test-suite boost-geometry-extensions-generic_robust_predicates + : + [ run approximate.cpp ] + [ run side3d.cpp : : : off ] + ; + diff --git a/extensions/test/generic_robust_predicates/approximate.cpp b/extensions/test/generic_robust_predicates/approximate.cpp new file mode 100644 index 0000000000..a7a4b76487 --- /dev/null +++ b/extensions/test/generic_robust_predicates/approximate.cpp @@ -0,0 +1,66 @@ +// 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::approximate_value; + using bg::detail::generic_robust_predicates::get_approx; + using bg::detail::generic_robust_predicates::approximate_interim; + using ct = CalculationType; + ct r1 = approximate_value, ct>( + std::array{1.0, 2.0}); + BOOST_CHECK_EQUAL(3.0, r1); + ct r2 = approximate_value, abs<_2>>, ct>( + 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}; + approximate_interim(r, input); + ct r3 = get_approx(r, input); + BOOST_CHECK_EQUAL(2.0, r3); + ct r4 = get_approx(r, input); + BOOST_CHECK_EQUAL(-6.0, r4); + ct r5 = get_approx(r, input); + BOOST_CHECK_EQUAL(-12.0, r5); +} + + +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/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..6dcfd615dc --- /dev/null +++ b/include/boost/geometry/extensions/generic_robust_predicates/strategies/cartesian/detail/almost_static_filter.hpp @@ -0,0 +1,133 @@ +// 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 +{ + +template +struct make_filter_impl +{ + template + static Filter apply(const ExtremaArray& extrema, const Reals&... args) + { + return make_filter_impl + ::apply(extrema, args..., extrema[N]); + } +}; + +template +struct make_filter_impl +{ + template + static Filter apply(const ExtremaArray& extrema, const Reals&... args) + { + Filter f(args...); + return f; + } +}; + +template +< + typename Expression, + typename CalculationType, + typename StaticFilter +> +class almost_static_filter +{ +private: + using expression_max_argn = max_argn; + using ct = CalculationType; + using extrema_array = std::array; + extrema_array m_extrema; + StaticFilter m_filter; +public: + const StaticFilter& 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(const Reals&... args) const + { + return m_filter.apply(args...); + } + + template + inline void update_extrema(const Reals&... args) + { + std::array input {{ static_cast(args)... }}; + for(int i = 0; i < m_extrema.size() / 2; ++i) + { + m_extrema[i] = std::max(m_extrema[i], input[i]); + } + for(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(const Reals&... args) + { + bool changed = false; + std::array input {{ static_cast(args)... }}; + for(int i = 0; i < m_extrema.size() / 2; ++i) + { + if(input[i] > m_extrema[i]) + { + changed = true; + m_extrema[i] = input[i]; + } + } + for(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; + } + + inline void update_filter() + { + m_filter = make_filter_impl + < + StaticFilter, + 0, + 2 * expression_max_argn::value + >::apply(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/approximate.hpp b/include/boost/geometry/extensions/generic_robust_predicates/strategies/cartesian/detail/approximate.hpp new file mode 100644 index 0000000000..65d4f100bf --- /dev/null +++ b/include/boost/geometry/extensions/generic_robust_predicates/strategies/cartesian/detail/approximate.hpp @@ -0,0 +1,478 @@ +// 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_APPROXIMATE_HPP +#define BOOST_GEOMETRY_EXTENSIONS_GENERIC_ROBUST_PREDICATES_STRATEGIES_CARTESIAN_DETAIL_APPROXIMATE_HPP + +#include +#include +#include + +#include +#include +#include + +#include + +namespace boost { namespace geometry +{ + +namespace detail { namespace generic_robust_predicates +{ + +template +< + typename Node, + std::size_t N, + typename Real, + typename InputArr +> +struct get_nth_real_impl +{ + static inline Real apply(const InputArr& input) + { + return input[N - 1]; + } +}; + +template +< + typename Node, + typename Real, + typename InputArr +> +struct get_nth_real_impl +{ + static inline Real apply(const InputArr&) + { + return Node::value; + } +}; + +template +< + typename Node, + std::size_t N, + typename Real, + typename InputArr +> +inline Real get_nth_real(const InputArr& input) +{ + return get_nth_real_impl::apply(input); +} + +template +< + typename All, + typename Node, + typename Real, + typename Arr, + bool is_leaf, + typename InputArr +> +struct get_approx_impl +{ + static inline Real apply(Arr& interim_results, const InputArr& input) + { + return interim_results[boost::mp11::mp_find::value]; + } +}; + +template +< + typename All, + typename Node, + typename Real, + typename Arr, + typename InputArr +> +struct get_approx_impl +{ + static inline Real apply(Arr& interim_results, const InputArr& input) + { + return get_nth_real(input); + } +}; + +template +< + typename All, + typename Node, + typename Real, + typename Arr, + typename InputArr +> +inline Real get_approx(Arr& interim_results, const InputArr& input) +{ + return get_approx_impl + < + All, + Node, + Real, + Arr, + is_leaf::value, + InputArr + >::apply(interim_results, input); +} + +template +< + typename All, + typename Remaining, + typename Real, + typename Arr, + operator_types Op, + typename InputArr +> +struct approximate_interim_impl {}; + +template +< + typename All, + typename Remaining, + typename Real, + typename Arr, + bool Empty, + typename InputArr +> +struct approximate_remainder_impl +{ + static inline void apply(Arr& interim_results, const InputArr& input) + { + using node = boost::mp11::mp_front; + approximate_interim_impl + < + All, + Remaining, + Real, + Arr, + node::operator_type, + InputArr + >::apply(interim_results, input); + } +}; + +template +< + typename All, + typename Remaining, + typename Real, + typename Arr, + typename InputArr +> +struct approximate_remainder_impl + < + All, + Remaining, + Real, + Arr, + true, + InputArr + > +{ + static inline void apply(Arr& interim_results, const InputArr& input) {} +}; + +template +< + typename All, + typename Remaining, + typename Real, + typename Arr, + typename InputArr +> +inline void approximate_remainder(Arr& interim_results, const InputArr& input) +{ + approximate_remainder_impl + < + All, + Remaining, + Real, + Arr, + boost::mp11::mp_empty::value, + InputArr + >::apply(interim_results, input); +} + +template +< + typename All, + typename Remaining, + typename Real, + typename Arr, + typename InputArr +> +struct approximate_interim_impl + < + All, + Remaining, + Real, + Arr, + operator_types::product, + InputArr + > +{ + static inline void apply(Arr& interim_results, const InputArr& input) + { + using node = boost::mp11::mp_front; + interim_results[boost::mp11::mp_find::value] = + get_approx(interim_results, + input) + * get_approx(interim_results, + input); + approximate_remainder + < + All, + boost::mp11::mp_pop_front, + Real + >(interim_results, input); + } +}; + +template +< + typename All, + typename Remaining, + typename Real, + typename Arr, + typename InputArr +> +struct approximate_interim_impl + < + All, + Remaining, + Real, + Arr, + operator_types::max, + InputArr + > +{ + static inline void apply(Arr& interim_results, const InputArr& input) + { + using node = boost::mp11::mp_front; + interim_results[boost::mp11::mp_find::value] = std::max( + get_approx(interim_results, + input), + get_approx(interim_results, + input)); + approximate_remainder + < + All, + boost::mp11::mp_pop_front, + Real + >(interim_results, input); + } +}; + +template +< + typename All, + typename Remaining, + typename Real, + typename Arr, + typename InputArr +> +struct approximate_interim_impl + < + All, + Remaining, + Real, + Arr, + operator_types::min, + InputArr + > +{ + static inline void apply(Arr& interim_results, const InputArr& input) + { + using node = boost::mp11::mp_front; + interim_results[boost::mp11::mp_find::value] = std::min( + get_approx(interim_results, + input), + get_approx(interim_results, + input)); + approximate_remainder + < + All, + boost::mp11::mp_pop_front, + Real + >(interim_results, input); + } +}; + +template +< + typename All, + typename Remaining, + typename Real, + typename Arr, + typename InputArr +> +struct approximate_interim_impl + < + All, + Remaining, + Real, + Arr, + operator_types::sum, + InputArr + > +{ + static inline void apply(Arr& interim_results, const InputArr& input) + { + using node = boost::mp11::mp_front; + interim_results[boost::mp11::mp_find::value] = + get_approx(interim_results, + input) + + get_approx(interim_results, + input); + approximate_remainder + < + All, + boost::mp11::mp_pop_front, + Real + >(interim_results, input); + } +}; + +template +< + typename All, + typename Remaining, + typename Real, + typename Arr, + typename InputArr +> +struct approximate_interim_impl + < + All, + Remaining, + Real, + Arr, + operator_types::difference, + InputArr + > +{ + static inline void apply(Arr& interim_results, const InputArr& input) + { + using node = boost::mp11::mp_front; + interim_results[boost::mp11::mp_find::value] = + get_approx(interim_results, + input) + - get_approx(interim_results, + input); + approximate_remainder + < + All, + boost::mp11::mp_pop_front, + Real + >(interim_results, input); + } +}; + +template +< + typename All, + typename Remaining, + typename Real, + typename Arr, + typename InputArr +> +struct approximate_interim_impl + < + All, + Remaining, + Real, + Arr, + operator_types::abs, + InputArr + > +{ + static inline void apply(Arr& interim_results, const InputArr& input) + { + using node = boost::mp11::mp_front; + interim_results[boost::mp11::mp_find::value] = + std::abs(get_approx + < + All, + typename node::child, + Real + >(interim_results, input)); + approximate_remainder + < + All, + boost::mp11::mp_pop_front, + Real + >(interim_results, input); + } +}; + +template +< + typename All, + typename Remaining, + typename Real, + typename Arr, + typename InputArr +> +struct approximate_interim_impl + < + All, + Remaining, + Real, + Arr, + operator_types::no_op, + InputArr + > +{ + static inline void apply(Arr& interim_results, const InputArr& input) + { + approximate_remainder + < + All, + boost::mp11::mp_pop_front, + Real + >(interim_results, input); + } +}; + +template +< + typename All, + typename Remaining, + typename Real, + typename Arr, + typename InputArr +> +inline void approximate_interim(Arr& interim_results, const InputArr& input) +{ + approximate_remainder + < + All, + Remaining, + Real + >(interim_results, input); +} + +template +inline Real approximate_value(const InputArr& input) +{ + using stack = typename boost::mp11::mp_unique>; + using evals = typename boost::mp11::mp_remove_if; + std::array::value> results; + approximate_interim(results, input); + return results.back(); +} + +}} // namespace detail::generic_robust_predicates + +}} // namespace boost::geometry + +#endif // BOOST_GEOMETRY_EXTENSIONS_GENERIC_ROBUST_PREDICATES_STRATEGIES_CARTESIAN_DETAIL_APPROXIMATE_HPP diff --git a/include/boost/geometry/extensions/generic_robust_predicates/strategies/cartesian/detail/coefficient_list.hpp b/include/boost/geometry/extensions/generic_robust_predicates/strategies/cartesian/detail/coefficient_list.hpp new file mode 100644 index 0000000000..46e39a0acf --- /dev/null +++ b/include/boost/geometry/extensions/generic_robust_predicates/strategies/cartesian/detail/coefficient_list.hpp @@ -0,0 +1,471 @@ +// 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_COEFFICIENT_LIST_HPP +#define BOOST_GEOMETRY_EXTENSIONS_GENERIC_ROBUST_PREDICATES_STRATEGIES_CARTESIAN_DETAIL_COEFFICIENT_LIST_HPP + +#include + +#include +#include + +#include + +namespace boost { namespace geometry +{ + +namespace detail { namespace generic_robust_predicates +{ + +//only meant to be used in assertions +template using is_mp_int = + boost::mp11::mp_same>; + +template using is_coeff_list = boost::mp11::mp_and + < + boost::mp11::mp_all_of, + boost::mp11::mp_similar> + >; + +template +struct coeff_truncate_impl +{ +private: + static_assert(is_coeff_list::value, "L must be a coefficient list"); + using first_nz = boost::mp11::mp_find_if; + using after_second_nz = + boost::mp11::mp_min< + boost::mp11::mp_size_t, + boost::mp11::mp_size + >; + using tail = boost::mp11::mp_erase_c; + using head = boost::mp11::mp_erase + < + L, + after_second_nz, + boost::mp11::mp_size + >; + using round_up = boost::mp11::mp_less + < + boost::mp11::mp_find_if, + boost::mp11::mp_size + >; +public: + using type = boost::mp11::mp_if + < + round_up, + boost::mp11::mp_push_back>, + head + >; + static_assert(is_coeff_list::value, "type should be a coefficient list"); +}; + +template using coeff_truncate = typename coeff_truncate_impl::type; + +template using app_zero_b = + boost::mp11::mp_push_front>; +template using app_zero_f = + boost::mp11::mp_push_back>; +template using mult_by_1_p_eps = coeff_truncate + < + boost::mp11::mp_transform + < + boost::mp11::mp_plus, + app_zero_b, + app_zero_f + > + >; + +template +< + typename L, + typename N, + typename done = boost::mp11::mp_bool +> +struct mult_by_1_p_eps_pow_impl +{ +private: + using next = mult_by_1_p_eps; +public: + using type = + typename mult_by_1_p_eps_pow_impl + < + next, + boost::mp11::mp_int + >::type; +}; + +template +struct mult_by_1_p_eps_pow_impl +{ +public: + using type = L; +}; + +template using mult_by_1_p_eps_pow = coeff_truncate + < + typename mult_by_1_p_eps_pow_impl::type + >; + +template using div_by_1_m_eps_helper = + boost::mp11::mp_partial_sum + < + L, + boost::mp11::mp_int<0>, + boost::mp11::mp_plus + >; + +template using div_by_1_m_eps = coeff_truncate + < + boost::mp11::mp_push_back + < + boost::mp11::mp_pop_back>, + inc>> + > + >; + +template +< + typename L1, + typename L2, + typename L, + typename L1_Empty = typename empty_or_void::type, + typename L2_Empty = typename empty_or_void::type +> +struct coeff_merge_impl {}; + +template +< + typename L1, + typename L2, + typename L +> +struct coeff_merge_impl +< + L1, + L2, + L, + boost::mp11::mp_false, + boost::mp11::mp_false +> +{ + using type = typename coeff_merge_impl + < + boost::mp11::mp_pop_front, + boost::mp11::mp_pop_front, + boost::mp11::mp_push_back + < + L, + boost::mp11::mp_plus + < + boost::mp11::mp_front, + boost::mp11::mp_front + > + > + >::type; +}; + +template +< + typename L1, + typename L2, + typename L +> +struct coeff_merge_impl + < + L1, + L2, + L, + boost::mp11::mp_true, + boost::mp11::mp_true + > +{ + using type = L; +}; + +template +< + typename L1, + typename L2, + typename L +> +struct coeff_merge_impl + < + L1, + L2, + L, + boost::mp11::mp_true, + boost::mp11::mp_false + > +{ + using type = boost::mp11::mp_append; +}; + +template +< + typename L1, + typename L2, + typename L +> +struct coeff_merge_impl + < + L1, + L2, + L, + boost::mp11::mp_false, + boost::mp11::mp_true + > +{ + using type = boost::mp11::mp_append; +}; + +template using coeff_merge = + coeff_truncate + < + typename coeff_merge_impl + < + L1, + L2, + boost::mp11::mp_list<> + >::type + >; + +template +< + typename L1, + typename L2, + typename L, + typename L1_Empty = typename empty_or_void::type, + typename L2_Empty = typename empty_or_void::type +> +struct coeff_max_impl {}; + +template +< + typename L1, + typename L2, + typename L +> +struct coeff_max_impl +{ + using type = typename coeff_max_impl + < + boost::mp11::mp_if + < + boost::mp11::mp_less + < + boost::mp11::mp_front, + boost::mp11::mp_front + >, + boost::mp11::mp_list<>, + boost::mp11::mp_pop_front + >, + boost::mp11::mp_if + < + boost::mp11::mp_less + < + boost::mp11::mp_front, + boost::mp11::mp_front + >, + boost::mp11::mp_list<>, + boost::mp11::mp_pop_front + >, + boost::mp11::mp_push_back + < + L, + boost::mp11::mp_max + < + boost::mp11::mp_front, + boost::mp11::mp_front + > + > + >::type; +}; + +template +< + typename L1, + typename L2, + typename L +> +struct coeff_max_impl + < + L1, + L2, + L, + boost::mp11::mp_true, + boost::mp11::mp_true + > +{ + using type = L; +}; + +template +< + typename L1, + typename L2, + typename L +> +struct coeff_max_impl + < + L1, + L2, + L, + boost::mp11::mp_true, + boost::mp11::mp_false + > +{ + using type = boost::mp11::mp_append; +}; + +template +< + typename L1, + typename L2, + typename L +> +struct coeff_max_impl + < + L1, + L2, + L, + boost::mp11::mp_false, + boost::mp11::mp_true + > +{ + using type = boost::mp11::mp_append; +}; + +template using coeff_max = + coeff_truncate + < + typename coeff_max_impl + < + L1, + L2, + boost::mp11::mp_list<> + >::type + >; + +template +< + typename L, + std::size_t Tail_Size = + boost::mp11::mp_size::value + - boost::mp11::mp_find_if::value +> +struct coeff_round_impl +{ +private: + using first_nz = boost::mp11::mp_find_if; + using tail = boost::mp11::mp_erase_c; + using zero_tail = boost::mp11::mp_same + < + boost::mp11::mp_find_if, + boost::mp11::mp_size + >; + using head = boost::mp11::mp_erase_c + < + L, + first_nz::value + 2, + boost::mp11::mp_size::value + >; + using minor = boost::mp11::mp_if + < + zero_tail, + boost::mp11::mp_back, + inc> + >; + using major = boost::mp11::mp_at; + using major_rounded = boost::mp11::mp_int< 1 << log_2_ceil::value >; + using minor_rounded = boost::mp11::mp_int + < + (minor::value / major_rounded::value) * major_rounded::value < minor::value ? + (minor::value / major_rounded::value + 1) * major_rounded::value + : (minor::value / major_rounded::value) * major_rounded::value + >; +public: + using type = boost::mp11::mp_push_back + < + boost::mp11::mp_pop_back, + minor_rounded + >; +}; + +template +struct coeff_round_impl { using type = L; }; + +template +struct coeff_round_impl { using type = L; }; + +template using coeff_round = typename coeff_round_impl::type; + +template using indexed_value_product = + boost::mp11::mp_list + < + boost::mp11::mp_plus + < + boost::mp11::mp_first, + boost::mp11::mp_first + >, + boost::mp11::mp_int + < + boost::mp11::mp_second::value + * boost::mp11::mp_second::value + > + >; + +template struct add_nested +{ + template using fn = + boost::mp11::mp_plus; +}; + +template using list_product_fold = + boost::mp11::mp_map_update_q + < + M, + IV, + add_nested< boost::mp11::mp_second > + >; + +template using index_less = + boost::mp11::mp_less + < + boost::mp11::mp_first, + boost::mp11::mp_first + >; + +template using list_product = + strip_index + < + boost::mp11::mp_sort + < + boost::mp11::mp_fold + < + boost::mp11::mp_product + < + indexed_value_product, + indexed, + indexed + >, + boost::mp11::mp_list<>, + list_product_fold + >, + index_less + > + >; + +}} // namespace detail::generic_robust_predicates + +}} // namespace boost::geometry + +#endif // BOOST_GEOMETRY_EXTENSIONS_GENERIC_ROBUST_PREDICATES_STRATEGIES_CARTESIAN_DETAIL_COEFFICIENT_LIST_HPP diff --git a/include/boost/geometry/extensions/generic_robust_predicates/strategies/cartesian/detail/error_bound.hpp b/include/boost/geometry/extensions/generic_robust_predicates/strategies/cartesian/detail/error_bound.hpp new file mode 100644 index 0000000000..6df2bafda6 --- /dev/null +++ b/include/boost/geometry/extensions/generic_robust_predicates/strategies/cartesian/detail/error_bound.hpp @@ -0,0 +1,605 @@ +// 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_ERROR_BOUND_HPP +#define BOOST_GEOMETRY_EXTENSIONS_GENERIC_ROBUST_PREDICATES_STRATEGIES_CARTESIAN_DETAIL_ERROR_BOUND_HPP + +#include + +#include +#include +#include + +#include + +namespace boost { namespace geometry +{ + +namespace detail { namespace generic_robust_predicates +{ + +template using second_is_coeff_list = + is_coeff_list< boost::mp11::mp_second >; + +template using is_error_map = boost::mp11::mp_and + < + boost::mp11::mp_is_map, + boost::mp11::mp_all_of, + boost::mp11::mp_similar + > + //TODO: It would be desirable to also validate that the keys are good +>; + +template +struct error_map_insert_impl +{ +private: + using key = boost::mp11::mp_front; + using value = boost::mp11::mp_second; + using other_value = typename mp_map_at_second_or_void::type; + using merged_value = coeff_merge; + using nkv = boost::mp11::mp_list; +public: + using type = boost::mp11::mp_map_replace; +}; + +template using error_map_insert = + typename error_map_insert_impl::type; + +template +struct add_fold_operator +{ +public: + template using fn = boost::mp11::mp_map_insert + < + M, + boost::mp11::mp_list + < + K, + coeff_merge + < + typename mp_map_at_second_or_void::type, + typename mp_map_at_second_or_void::type + > + > + >; +}; + +template using add_children = boost::mp11::mp_fold + < + boost::mp11::mp_set_union + < + boost::mp11::mp_map_keys, + boost::mp11::mp_map_keys + >, + boost::mp11::mp_list<>, + add_fold_operator::template fn + >; + +template +< + typename Exp, + typename EM, + typename Out, + typename LErr = typename val_or_empty_list::type, + typename RErr = typename val_or_empty_list::type, + typename skip_decompose = + boost::mp11::mp_or + < + boost::mp11::mp_not + < + boost::mp11::mp_same + < + typename Exp::error_type, + sum_error_type + > + >, + boost::mp11::mp_same + < + boost::mp11::mp_list<>, + LErr, + RErr + > + > +> +struct decompose_add_impl +{ +private: + using decomp_left = typename decompose_add_impl + < + typename Exp::left, + EM, + Out + >::type; + using decomp_right = typename decompose_add_impl + < + typename Exp::right, + EM, + decomp_left + >::type; +public: + using type = decomp_right; +}; + +template +< + typename Exp, + typename EM, + typename Out, + typename LErr, + typename RErr +> +struct decompose_add_impl + < + Exp, + EM, + Out, + LErr, + RErr, + boost::mp11::mp_true + > +{ + using type = error_map_insert + < + boost::mp11::mp_list + < + Exp, + boost::mp11::mp_list> + >, + Out + >; +}; + +template +< + typename Exp, + typename EM, + typename Out +> +using decompose_add = typename decompose_add_impl::type; + +template +< + typename Exp, + typename EM, + typename LErr = typename val_or_empty_list::type, + typename RErr = typename val_or_empty_list::type, + typename Children_Empty = + boost::mp11::mp_and + < + typename empty_or_void::type, + typename empty_or_void::type + > +> +struct sum_err_impl +{ +private: + static_assert(is_error_map::value, "LErr needs to be a valid error map."); + static_assert(is_error_map::value, "RErr needs to be a valid error map."); + using children = add_children; + static_assert( + is_error_map::value, "children needs to be a valid error map."); +public: + /* + using type = boost::mp11::mp_map_insert< + children, + boost::mp11::mp_list>> + >;*/ + using type = decompose_add; + static_assert(is_error_map::value, "type needs to be a valid error map."); +}; + +template +< + typename Exp, + typename EM, + typename LErr, + typename RErr +> +struct sum_err_impl + < + Exp, + EM, + LErr, + RErr, + boost::mp11::mp_true + > +{ + static_assert(is_error_map::value, "LErr needs to be a valid error map."); + static_assert(is_error_map::value, "RErr needs to be a valid error map."); + using type = boost::mp11::mp_list + < + boost::mp11::mp_list + < + Exp, + boost::mp11::mp_list> + > + >; + static_assert(is_error_map::value, "type needs to be a valid error map."); +}; + +template using sum_err = + typename sum_err_impl::type; + +template using pad_second = boost::mp11::mp_list + < + boost::mp11::mp_front, + boost::mp11::mp_push_front + < + boost::mp11::mp_second, + boost::mp11::mp_int<0> + > + >; + +template using pop_front_second = boost::mp11::mp_list + < + boost::mp11::mp_front, + boost::mp11::mp_pop_front> + >; + +template using increment_first_of_second = + boost::mp11::mp_transform_front; + +template using prod_entry_merge = + boost::mp11::mp_list + < + boost::mp11::mp_flatten + < + boost::mp11::mp_list + < + boost::mp11::mp_first, + boost::mp11::mp_first + > + >, + list_product + < + boost::mp11::mp_second, + boost::mp11::mp_second + > + >; + +template +< + typename Exp, + typename EM, + typename LErr = typename val_or_empty_list::type, + typename RErr = typename val_or_empty_list::type +> +struct prod_children_impl +{ +private: + static_assert(is_error_map::value, "LErr needs to be a valid error map."); + static_assert(is_error_map::value, "RErr needs to be a valid error map."); + using left = typename Exp::left; + using right = typename Exp::right; + using padded_lerr = boost::mp11::mp_transform; + using added_left_lerr = decompose_add; + using padded_rerr = boost::mp11::mp_transform; + using added_right_rerr = decompose_add; + using prod = boost::mp11::mp_product + < + prod_entry_merge, + added_left_lerr, + added_right_rerr + >; + using stripped_prod = boost::mp11::mp_transform; +public: + using type = stripped_prod; + static_assert(is_error_map::value, "type needs to be a valid error map."); +}; + +template using prod_children = + typename prod_children_impl::type; + +template +< + typename Exp, + typename EM, + typename LErr = typename val_or_empty_list::type, + typename RErr = typename val_or_empty_list::type +> +struct product_err_impl +{ +private: + static_assert(is_error_map::value, "LErr needs to be a valid error map."); + static_assert(is_error_map::value, "RErr needs to be a valid error map."); + using children = prod_children; + static_assert(is_error_map::value, "children needs to be a valid error map."); +public: + using type = boost::mp11::mp_map_insert + < + children, + boost::mp11::mp_list + < + Exp, + boost::mp11::mp_list< boost::mp11::mp_int<1> > + > + >; + static_assert(is_error_map::value, "type needs to be a valid error map."); +}; + +template using product_err = + typename product_err_impl::type; + +template +< + typename Errors, + typename Exp, + typename IsSum = boost::mp11::mp_same + < + typename Exp::error_type, + sum_error_type + > +> +struct sum_or_product_err_impl +{ + using type = sum_err; +}; + +template +< + typename Errors, + typename Exp +> +struct sum_or_product_err_impl +{ + using type = product_err; +}; + +template using sum_or_product_err = + typename sum_or_product_err_impl::type; + +template +< + typename Errors, + typename Exp +> +struct error_fold_impl +{ +private: + using lerr = typename val_or_empty_list::type; + using rerr = typename val_or_empty_list::type; + using err = sum_or_product_err; +public: + using type = boost::mp11::mp_map_insert + < + Errors, + boost::mp11::mp_list + >; +}; + + +template using error_fold = + typename error_fold_impl::type; + +template using evals_error = + boost::mp11::mp_fold + < + Evals, + boost::mp11::mp_list<>, + error_fold + >; + +template using is_mp_list = boost::mp11::mp_similar + < + boost::mp11::mp_list<>, + T + >; + +template +struct list_to_product_impl +{ +private: + using key = boost::mp11::mp_front; + using value = boost::mp11::mp_second; + using multiplications = boost::mp11::mp_int::value - 1>; + using nvalue = mult_by_1_p_eps_pow; + using nkey = boost::mp11::mp_reverse_fold< + boost::mp11::mp_pop_back, + boost::mp11::mp_back, + product + >; +public: + using type = boost::mp11::mp_list; +}; + +template using list_to_product = typename list_to_product_impl::type; + +template +< + typename M, + typename KV, + typename KeyMPList = is_mp_list< boost::mp11::mp_front > +> +struct error_map_list_to_product_fold_impl +{ + using type = error_map_insert, M>; +}; + +template +struct error_map_list_to_product_fold_impl + < + M, + KV, + boost::mp11::mp_false + > +{ + using type = error_map_insert; +}; + +template using error_map_list_to_product_fold = + typename error_map_list_to_product_fold_impl::type; + +template +struct error_map_list_to_product_impl +{ + using type = boost::mp11::mp_fold + < + M, + boost::mp11::mp_list<>, + error_map_list_to_product_fold + >; +}; + +template using error_map_list_to_product = + typename error_map_list_to_product_impl::type; + +template +struct abs_unless_non_negative +{ + using type = boost::mp11::mp_if_c + < + Expression::non_negative, + Expression, + abs + >; +}; + +template using abs_fold = + boost::mp11::mp_push_back + < + EM, + boost::mp11::mp_list + < + typename abs_unless_non_negative> + ::type, + boost::mp11::mp_second + > + >; + +template using abs_all = + boost::mp11::mp_fold + < + EM, + boost::mp11::mp_list<>, + abs_fold + >; + +template +< + typename KV1, + typename KV2 +> +struct abs_sum_error_term_impl +{ +private: + using key1 = boost::mp11::mp_front; + using key2 = boost::mp11::mp_front; + using nkey = sum; + using val1 = boost::mp11::mp_second; + using val2 = boost::mp11::mp_second; + using mval = coeff_max; + static_assert(is_coeff_list::value, "merged coefficient list fails coefficient list check"); + using nval = mult_by_1_p_eps; +public: + using type = boost::mp11::mp_list; +}; + +template using abs_sum_error_term = typename abs_sum_error_term_impl::type; + +//TODO: The following could be probably optimized for potentially produce better error bounds in some cases +// if the error map is treated as a minheap by ordering of epsilon-polynomial. +template using error_map_sum_up = + boost::mp11::mp_fold + < + boost::mp11::mp_pop_front, + boost::mp11::mp_first, + abs_sum_error_term + >; + +template::value> +struct error_map_balanced_sum_up_impl +{ +private: + static constexpr std::size_t mid = MS/2; + using left = typename error_map_balanced_sum_up_impl + < + boost::mp11::mp_erase_c + >::type; + using right = typename error_map_balanced_sum_up_impl + < + boost::mp11::mp_erase_c + >::type; +public: + using type = abs_sum_error_term; +}; + +template +struct error_map_balanced_sum_up_impl +{ + using type = boost::mp11::mp_front; +}; + +template using error_map_balanced_sum_up = + typename error_map_balanced_sum_up_impl::type; + +template +< + typename Real, + typename Exp +> +struct eps_pow +{ + static constexpr Real value = + std::numeric_limits::epsilon()/2.0 + * eps_pow>::value; +}; + +template struct eps_pow> +{ + static constexpr Real value = 1.0; +}; + +template +< + typename Real, + typename L, + typename S = boost::mp11::mp_size +> +struct eval_eps_polynomial +{ +private: + using last = boost::mp11::mp_back; + using s2last = boost::mp11::mp_back< boost::mp11::mp_pop_back >; +public: + static constexpr Real value = + s2last::value + * eps_pow>::value + + last::value * eps_pow::value; +}; + +template +< + typename Real, + typename L +> +struct eval_eps_polynomial> +{ + static constexpr Real value = + boost::mp11::mp_front::value + * std::numeric_limits::epsilon() / 2.0; +}; + +template +< + typename Real, + typename L +> +struct eval_eps_polynomial> +{ + static constexpr Real value = 0; +}; + +}} // namespace detail::generic_robust_predicates + +}} // namespace boost::geometry + +#endif // BOOST_GEOMETRY_EXTENSIONS_GENERIC_ROBUST_PREDICATES_STRATEGIES_CARTESIAN_DETAIL_ERROR_BOUND_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..4ac3e558b0 --- /dev/null +++ b/include/boost/geometry/extensions/generic_robust_predicates/strategies/cartesian/detail/expression_tree.hpp @@ -0,0 +1,250 @@ +// 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 + +namespace boost { namespace geometry +{ + +namespace detail { namespace generic_robust_predicates +{ + +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; + 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; + +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 > +struct max_argn_impl; + +template using max_argn = typename max_argn_impl::type; + +template +struct max_argn_impl +{ +private: + using children_list = boost::mp11::mp_rename; + using children_max_argn = + boost::mp11::mp_transform; +public: + using type = boost::mp11::mp_max_element + < + children_max_argn, + boost::mp11::mp_less + >; +}; + +template +struct max_argn_impl +{ + using type = boost::mp11::mp_size_t; +}; + +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>; + +}} // 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..6042b1834e --- /dev/null +++ b/include/boost/geometry/extensions/generic_robust_predicates/strategies/cartesian/detail/expressions.hpp @@ -0,0 +1,103 @@ +// 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 + +namespace boost { namespace geometry +{ + +namespace detail { namespace generic_robust_predicates +{ + +template +< + typename A11, typename A12, + typename A21, typename A22 +> +using det2x2 = difference + < + product, + product + >; + +using orient2d = det2x2 + < + difference <_1, _5>, difference<_2, _6>, + difference <_3, _5>, difference<_4, _6> + >; + +template +< + typename A11, typename A12, typename A13, + typename A21, typename A22, typename A23, + typename A31, typename A32, typename A33 +> +struct det3x3_helper +{ +private: + using minor1 = product>; + using minor2 = product>; + using minor3 = product>; +public: + using type = sum>; +}; + +template +< + typename A11, typename A12, typename A13, + typename A21, typename A22, typename A23, + typename A31, typename A32, typename A33 +> +using det3x3 = typename det3x3_helper + < + A11, A12, A13, + A21, A22, A23, + A31, A32, A33 + >::type; + +using orient3d = det3x3 + < + difference<_1, _10>, difference<_2, _11>, difference<_3, _12>, + difference<_4, _10>, difference<_5, _11>, difference<_6, _12>, + difference<_7, _10>, difference<_8, _11>, difference<_9, _12> + >; + +struct incircle_helper +{ +private: + using adx = difference<_1, _7>; + using ady = difference<_2, _8>; + using bdx = difference<_3, _7>; + using bdy = difference<_4, _8>; + using cdx = difference<_5, _7>; + using cdy = difference<_6, _8>; + using alift = sum, product>; + using blift = sum, product>; + using clift = sum, product>; +public: + using type = det3x3 + < + alift, adx, ady, + blift, bdx, bdy, + clift, cdx, cdy + >; +}; + +using incircle = incircle_helper::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..59e4d3f47b --- /dev/null +++ b/include/boost/geometry/extensions/generic_robust_predicates/strategies/cartesian/detail/interval_error_bound.hpp @@ -0,0 +1,243 @@ +// 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 +{ + +template +struct interval_min_impl +{ + using type = Expression; +}; + +template +struct interval_max_impl +{ + using type = Expression; +}; + +template +struct interval_min_impl, max_argn> +{ +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, max_argn> +{ +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, max_argn> +{ +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, max_argn> +{ +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, max_argn> +{ +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, max_argn> +{ +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, max_argn> +{ +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, max_argn> +{ +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, max_argn> +{ +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, max_argn> +{ +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, max_argn> +{ + using type = argument; +}; + +template +struct interval_max_impl, max_argn> +{ + using type = argument; +}; + +template +struct interval_impl +{ + using type = Expression; +}; + +template +struct interval_impl, max_argn> +{ +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, max_argn> +{ +private: + using left = typename interval_impl::type; + using right = typename interval_impl::type; +public: + using type = difference; +}; + +template +struct interval_impl, max_argn> +{ +private: + using left = typename interval_impl::type; + using right = typename interval_impl::type; +public: + using type = sum; +}; + +template +struct interval_impl, max_argn> +{ +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::value>::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..26ca2884fa --- /dev/null +++ b/include/boost/geometry/extensions/generic_robust_predicates/strategies/cartesian/detail/semi_static_filter.hpp @@ -0,0 +1,91 @@ +// 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 +{ + +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: + template + static inline int apply(const Reals&... args) + { + std::array input + {{ static_cast(args)... }}; + std::array::value> results; + approximate_interim(results, input); + const ct error_bound = + get_approx(results, input); + const ct det = + get_approx(results, input); + 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..96725db461 --- /dev/null +++ b/include/boost/geometry/extensions/generic_robust_predicates/strategies/cartesian/detail/stage_a.hpp @@ -0,0 +1,141 @@ +// 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 +#include + +#include +#include +#include + +namespace boost { namespace geometry +{ + +namespace detail { namespace generic_robust_predicates +{ + +template +struct stage_a_error_bound_impl +{ +private: + using ct = CalculationType; + using root = Expression; + using stack = typename boost::mp11::mp_unique>; + using evals = typename boost::mp11::mp_remove_if; + using interim_evals = typename boost::mp11::mp_remove + < + boost::mp11::mp_remove_if, + root + >; + using interim_errors = evals_error; + using final_children = add_children + < + boost::mp11::mp_second + < + boost::mp11::mp_map_find + < + interim_errors, + typename root::left + > + >, + boost::mp11::mp_second + < + boost::mp11::mp_map_find + < + interim_errors, + typename root::right + > + > + >; + using final_children_ltp = error_map_list_to_product; + using final_children_ltp_abs = abs_all; + using final_children_sum_fold = error_map_balanced_sum_up; + using final_coeff = coeff_round + < + div_by_1_m_eps + < + mult_by_1_p_eps + < + boost::mp11::mp_second + > + > + >; + static constexpr ct final_coeff_value = eval_eps_polynomial::value; + + struct final_coeff_constant : public static_constant_interface + { + static constexpr ct value = final_coeff_value; + static constexpr bool non_negative = true; + }; + + using error_expression_variable = boost::mp11::mp_front; +public: + using type = product; +}; + +template +< + typename Expression, + typename CalculationType +> +using stage_a_error_bound = + typename stage_a_error_bound_impl::type; + +template +< + typename Expression, + typename CalculationType +> +using stage_a_semi_static = semi_static_filter + < + Expression, + CalculationType, + stage_a_error_bound + >; + +template +< + typename Expression, + typename CalculationType +> +using stage_a_static = static_filter + < + Expression, + CalculationType, + interval> + >; + +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/static_filter.hpp b/include/boost/geometry/extensions/generic_robust_predicates/strategies/cartesian/detail/static_filter.hpp new file mode 100644 index 0000000000..d93eb6497f --- /dev/null +++ b/include/boost/geometry/extensions/generic_robust_predicates/strategies/cartesian/detail/static_filter.hpp @@ -0,0 +1,98 @@ +// 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 +{ + +template +< + typename Expression, + typename CalculationType, + typename ErrorExpression +> +class static_filter +{ +private: + using ct = CalculationType; + 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 + >; + ct m_error_bound; +public: + inline static_filter() {} + + ct error_bound() const { return m_error_bound; } + + template + inline static_filter(const Reals&... args) + : m_error_bound(approximate_value( + std::array + {static_cast(args)...})) + { + static_assert(sizeof...(Reals) == max_argn::value, + "Number of constructor arguments is incompatible with error expression."); + } + + template + inline int apply(const Reals&... args) const + { + std::array input {static_cast(args)...}; + std::array::value> results; + approximate_interim(results, input); + const ct det = get_approx(results, 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()(const Reals&... 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 diff --git a/include/boost/geometry/extensions/generic_robust_predicates/strategies/cartesian/detail/static_util.hpp b/include/boost/geometry/extensions/generic_robust_predicates/strategies/cartesian/detail/static_util.hpp new file mode 100644 index 0000000000..c23f548645 --- /dev/null +++ b/include/boost/geometry/extensions/generic_robust_predicates/strategies/cartesian/detail/static_util.hpp @@ -0,0 +1,130 @@ +// 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_UTIL_HPP +#define BOOST_GEOMETRY_EXTENSIONS_GENERIC_ROBUST_PREDICATES_STRATEGIES_CARTESIAN_DETAIL_STATIC_UTIL_HPP + +#include +#include +#include + + +namespace boost { namespace geometry +{ + +namespace detail { namespace generic_robust_predicates +{ + +template +< + typename LOV, + typename Void = boost::mp11::mp_same +> +struct empty_or_void +{ + using type = boost::mp11::mp_empty; +}; + +template +struct empty_or_void +{ + using type = boost::mp11::mp_true; +}; + +template +struct log_2_floor_impl +{ + using type = typename boost::mp11::mp_plus + < + typename log_2_floor_impl>::type, + boost::mp11::mp_int<1> + >; +}; + +template <> +struct log_2_floor_impl> +{ + using type = boost::mp11::mp_int<0>; +}; + +template using log_2_floor = typename log_2_floor_impl::type; + +template +struct log_2_ceil_impl +{ +private: + using floor = log_2_floor; +public: + using type = boost::mp11::mp_int + < + (1 << floor::value) == N::value ? floor::value : floor::value + 1 + >; +}; + +template using log_2_ceil = typename log_2_ceil_impl::type; + +template using is_not_zero = boost::mp11::mp_bool; + +template +< + typename M, + typename K, + typename Contained = boost::mp11::mp_map_contains +> +struct mp_map_at_second_or_void +{ + using type = void; +}; + +template +struct mp_map_at_second_or_void +{ + using type = boost::mp11::mp_second>; +}; + +template using inc = boost::mp11::mp_int; + +template using indexed = boost::mp11::mp_transform + < + boost::mp11::mp_list, + boost::mp11::mp_iota>, + L + >; + +template using strip_index = + boost::mp11::mp_transform; + +template +< + typename Map, + typename Key, + typename Contains = boost::mp11::mp_map_contains +> +struct val_or_empty_list +{ + using type = boost::mp11::mp_second>; +}; + +template +< + typename Map, + typename Key +> +struct val_or_empty_list +{ + using type = boost::mp11::mp_list<>; +}; + +}} // namespace detail::generic_robust_predicates + +}} // namespace boost::geometry + +#endif // BOOST_GEOMETRY_EXTENSIONS_GENERIC_ROBUST_PREDICATES_STRATEGIES_CARTESIAN_DETAIL_STATIC_UTIL_HPP From aa2afbc62aa67604c550cc73c3ff1520e0e4e7b5 Mon Sep 17 00:00:00 2001 From: Tinko Bartels Date: Tue, 11 Aug 2020 09:56:07 +0200 Subject: [PATCH 2/6] Style fixes, removal of unused parameters and addition of some explanations in comments. --- .../static_side_2d.cpp | 4 +- .../cartesian/detail/almost_static_filter.hpp | 24 ++++++++-- .../cartesian/detail/approximate.hpp | 47 ++++++++++++------- .../cartesian/detail/coefficient_list.hpp | 32 +++++++++++++ .../cartesian/detail/error_bound.hpp | 18 +++++++ .../cartesian/detail/expression_tree.hpp | 19 ++++++++ .../cartesian/detail/expressions.hpp | 8 ++-- .../cartesian/detail/interval_error_bound.hpp | 7 +++ .../cartesian/detail/semi_static_filter.hpp | 18 ++++++- .../strategies/cartesian/detail/stage_a.hpp | 3 ++ .../cartesian/detail/static_filter.hpp | 16 +++++-- 11 files changed, 163 insertions(+), 33 deletions(-) diff --git a/extensions/example/generic_robust_predicates/static_side_2d.cpp b/extensions/example/generic_robust_predicates/static_side_2d.cpp index c18ef3c4f1..04ba5b0515 100644 --- a/extensions/example/generic_robust_predicates/static_side_2d.cpp +++ b/extensions/example/generic_robust_predicates/static_side_2d.cpp @@ -56,7 +56,7 @@ struct side_robust_with_static_filter bg::get<1>(p2), bg::get<0>(p), bg::get<1>(p)); - if(sign != bg::detail::generic_robust_predicates::sign_uncertain) + if (sign != bg::detail::generic_robust_predicates::sign_uncertain) { return sign; } @@ -68,7 +68,7 @@ struct side_robust_with_static_filter } }; -int main(int argc, char** argv) +int main() { point p1(0.0, 0.0); point p2(1.0, 1.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 index 6dcfd615dc..7727c1c4bc 100644 --- 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 @@ -28,7 +28,7 @@ template struct make_filter_impl { template - static Filter apply(const ExtremaArray& extrema, const Reals&... args) + static Filter apply(ExtremaArray const& extrema, Reals const&... args) { return make_filter_impl ::apply(extrema, args..., extrema[N]); @@ -39,13 +39,27 @@ template struct make_filter_impl { template - static Filter apply(const ExtremaArray& extrema, const Reals&... args) + static Filter apply(ExtremaArray const& extrema, Reals const&... args) { Filter f(args...); return f; } }; +//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, @@ -72,13 +86,13 @@ class almost_static_filter std::numeric_limits::infinity()); } template - int apply(const Reals&... args) const + int apply(Reals const&... args) const { return m_filter.apply(args...); } template - inline void update_extrema(const Reals&... args) + inline void update_extrema(Reals const&... args) { std::array input {{ static_cast(args)... }}; for(int i = 0; i < m_extrema.size() / 2; ++i) @@ -92,7 +106,7 @@ class almost_static_filter } template - inline bool update_extrema_check(const Reals&... args) + inline bool update_extrema_check(Reals const&... args) { bool changed = false; std::array input {{ static_cast(args)... }}; diff --git a/include/boost/geometry/extensions/generic_robust_predicates/strategies/cartesian/detail/approximate.hpp b/include/boost/geometry/extensions/generic_robust_predicates/strategies/cartesian/detail/approximate.hpp index 65d4f100bf..53dabca59e 100644 --- a/include/boost/geometry/extensions/generic_robust_predicates/strategies/cartesian/detail/approximate.hpp +++ b/include/boost/geometry/extensions/generic_robust_predicates/strategies/cartesian/detail/approximate.hpp @@ -28,6 +28,10 @@ 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. +//The most important template in this file is approximate_interim. + template < typename Node, @@ -37,7 +41,7 @@ template > struct get_nth_real_impl { - static inline Real apply(const InputArr& input) + static inline Real apply(InputArr const& input) { return input[N - 1]; } @@ -51,7 +55,7 @@ template > struct get_nth_real_impl { - static inline Real apply(const InputArr&) + static inline Real apply(InputArr const&) { return Node::value; } @@ -64,7 +68,7 @@ template typename Real, typename InputArr > -inline Real get_nth_real(const InputArr& input) +inline Real get_nth_real(InputArr const& input) { return get_nth_real_impl::apply(input); } @@ -80,7 +84,7 @@ template > struct get_approx_impl { - static inline Real apply(Arr& interim_results, const InputArr& input) + static inline Real apply(Arr& interim_results, InputArr const&) { return interim_results[boost::mp11::mp_find::value]; } @@ -96,7 +100,7 @@ template > struct get_approx_impl { - static inline Real apply(Arr& interim_results, const InputArr& input) + static inline Real apply(Arr&, InputArr const& input) { return get_nth_real(input); } @@ -110,7 +114,7 @@ template typename Arr, typename InputArr > -inline Real get_approx(Arr& interim_results, const InputArr& input) +inline Real get_approx(Arr& interim_results, InputArr const& input) { return get_approx_impl < @@ -145,7 +149,7 @@ template > struct approximate_remainder_impl { - static inline void apply(Arr& interim_results, const InputArr& input) + static inline void apply(Arr& interim_results, InputArr const& input) { using node = boost::mp11::mp_front; approximate_interim_impl @@ -178,7 +182,7 @@ struct approximate_remainder_impl InputArr > { - static inline void apply(Arr& interim_results, const InputArr& input) {} + static inline void apply(Arr&, InputArr const&) {} }; template @@ -189,7 +193,7 @@ template typename Arr, typename InputArr > -inline void approximate_remainder(Arr& interim_results, const InputArr& input) +inline void approximate_remainder(Arr& interim_results, InputArr const& input) { approximate_remainder_impl < @@ -220,7 +224,7 @@ struct approximate_interim_impl InputArr > { - static inline void apply(Arr& interim_results, const InputArr& input) + static inline void apply(Arr& interim_results, InputArr const& input) { using node = boost::mp11::mp_front; interim_results[boost::mp11::mp_find::value] = @@ -255,7 +259,7 @@ struct approximate_interim_impl InputArr > { - static inline void apply(Arr& interim_results, const InputArr& input) + static inline void apply(Arr& interim_results, InputArr const& input) { using node = boost::mp11::mp_front; interim_results[boost::mp11::mp_find::value] = std::max( @@ -290,7 +294,7 @@ struct approximate_interim_impl InputArr > { - static inline void apply(Arr& interim_results, const InputArr& input) + static inline void apply(Arr& interim_results, InputArr const& input) { using node = boost::mp11::mp_front; interim_results[boost::mp11::mp_find::value] = std::min( @@ -325,7 +329,7 @@ struct approximate_interim_impl InputArr > { - static inline void apply(Arr& interim_results, const InputArr& input) + static inline void apply(Arr& interim_results, InputArr const& input) { using node = boost::mp11::mp_front; interim_results[boost::mp11::mp_find::value] = @@ -360,7 +364,7 @@ struct approximate_interim_impl InputArr > { - static inline void apply(Arr& interim_results, const InputArr& input) + static inline void apply(Arr& interim_results, InputArr const& input) { using node = boost::mp11::mp_front; interim_results[boost::mp11::mp_find::value] = @@ -395,7 +399,7 @@ struct approximate_interim_impl InputArr > { - static inline void apply(Arr& interim_results, const InputArr& input) + static inline void apply(Arr& interim_results, InputArr const& input) { using node = boost::mp11::mp_front; interim_results[boost::mp11::mp_find::value] = @@ -432,7 +436,7 @@ struct approximate_interim_impl InputArr > { - static inline void apply(Arr& interim_results, const InputArr& input) + static inline void apply(Arr& interim_results, InputArr const& input) { approximate_remainder < @@ -443,6 +447,13 @@ struct approximate_interim_impl } }; +//All expects an boost::mp11::mp_list of all expressions that need to be +//evaluated. Remaining expects an boost::mp11::mp_list of the expressions +//that are left to be evaluated. In the first call it is expected to be +//equal to All and it serves as an anchor for template recursion. Real is a +//floating-point type. The remaining template arguments are deduced from +//parameters. + template < typename All, @@ -451,7 +462,7 @@ template typename Arr, typename InputArr > -inline void approximate_interim(Arr& interim_results, const InputArr& input) +inline void approximate_interim(Arr& interim_results, InputArr const& input) { approximate_remainder < @@ -462,7 +473,7 @@ inline void approximate_interim(Arr& interim_results, const InputArr& input) } template -inline Real approximate_value(const InputArr& input) +inline Real approximate_value(InputArr const& input) { using stack = typename boost::mp11::mp_unique>; using evals = typename boost::mp11::mp_remove_if; diff --git a/include/boost/geometry/extensions/generic_robust_predicates/strategies/cartesian/detail/coefficient_list.hpp b/include/boost/geometry/extensions/generic_robust_predicates/strategies/cartesian/detail/coefficient_list.hpp index 46e39a0acf..331f983e2f 100644 --- a/include/boost/geometry/extensions/generic_robust_predicates/strategies/cartesian/detail/coefficient_list.hpp +++ b/include/boost/geometry/extensions/generic_robust_predicates/strategies/cartesian/detail/coefficient_list.hpp @@ -25,6 +25,18 @@ namespace boost { namespace geometry namespace detail { namespace generic_robust_predicates { +//The templates in this file are helpers for the computation of error maps and +//error expressions (see also error_bound.hpp). A coefficient list is a +//boost::mp11::mp_list that contains integers that represent the coefficients +//of an polynomial. This is used for the polynomials in epsilon that occur in +//forward error analysis, e.g. epsilon + 2 * epsilon^2 + epsilon^3 would be +//represented by the type mp_list, mp_int<2>, mp_int<1>>. +// +//See p. 39 in "Adaptive Precision Floating-Point Arithmetic and Fast Robust +//Geometric Predicates" by Jonathan Richard Shewchuk (can be downloaded from +//https://www.cs.cmu.edu/~quake/robust.html) for the inspiration for the +//methods implemented below. + //only meant to be used in assertions template using is_mp_int = boost::mp11::mp_same>; @@ -68,12 +80,21 @@ struct coeff_truncate_impl static_assert(is_coeff_list::value, "type should be a coefficient list"); }; +//Because of the way floating-point numbers are rounded, for a given epsilon +//polynomial we only care about the first two non-zero coefficients. +//e.g. if a*eps + b*eps^2 + c*eps^3 is considered then c * eps^3 is always +//neglible unless c is so large that it must have resulted from the analysis +//of an expression that is so complex that it cannot be feasibly processed +//by the methods implemented in this feature. template using coeff_truncate = typename coeff_truncate_impl::type; template using app_zero_b = boost::mp11::mp_push_front>; template using app_zero_f = boost::mp11::mp_push_back>; + +//For a coefficient list representing the polynomial p, the following template +//produces the coefficient list representing p * (1 + eps). template using mult_by_1_p_eps = coeff_truncate < boost::mp11::mp_transform @@ -110,6 +131,9 @@ struct mult_by_1_p_eps_pow_impl using type = L; }; + +//For a coefficient list representing the polynomial p, the following template +//produces the cofficient list representing p * (1 + eps)^N template using mult_by_1_p_eps_pow = coeff_truncate < typename mult_by_1_p_eps_pow_impl::type @@ -123,6 +147,9 @@ template using div_by_1_m_eps_helper = boost::mp11::mp_plus >; +//For a coefficient list representing a polynomial p in eps, the following +//template computes a coefficient list representing a polynomial q such that +//q * (1 - p) >= p. template using div_by_1_m_eps = coeff_truncate < boost::mp11::mp_push_back @@ -406,6 +433,11 @@ struct coeff_round_impl { using type = L; }; template struct coeff_round_impl { using type = L; }; +//The following template rounds a coefficient list up such that the +//corresponding polynomial evaluated for epsilon is a representable +//floating-point number. +// +//e.g. 3 * eps + 12 * eps^2 + 1 * eps^3 is rounded to 3 * eps + 16 * eps^2 template using coeff_round = typename coeff_round_impl::type; template using indexed_value_product = diff --git a/include/boost/geometry/extensions/generic_robust_predicates/strategies/cartesian/detail/error_bound.hpp b/include/boost/geometry/extensions/generic_robust_predicates/strategies/cartesian/detail/error_bound.hpp index 6df2bafda6..d42a519aa9 100644 --- a/include/boost/geometry/extensions/generic_robust_predicates/strategies/cartesian/detail/error_bound.hpp +++ b/include/boost/geometry/extensions/generic_robust_predicates/strategies/cartesian/detail/error_bound.hpp @@ -26,6 +26,24 @@ namespace boost { namespace geometry namespace detail { namespace generic_robust_predicates { +//The templates in this file are utilities for the generation of error +//expressions for given arithmetic expressions. This is done using forward +//error analysis methods at compile time. This analysis produces an error map +//that is folded into an error expression. +// +//An Error Map is a template fulfilling the boost::mp11 map concept. The keys K +//are expressions (as in expression trees and the values V are lists of +//integers. Each key-value pair represents the contribution of an error term to +//the total error, i.e. the error expression is something like the sum of +//abs * (V[0] * epsilon + V[1] * epsilon^2 + ...) +// +//over all key-value pairs KV in the error map. +// +//See p. 39 in "Adaptive Precision Floating-Point Arithmetic and Fast Robust +//Geometric Predicates" by Jonathan Richard Shewchuk (can be downloaded from +//https://www.cs.cmu.edu/~quake/robust.html) for the inspiration for the +//methods implemented below. + template using second_is_coeff_list = is_coeff_list< boost::mp11::mp_second >; 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 index 4ac3e558b0..91bbd6d3d8 100644 --- 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 @@ -25,6 +25,19 @@ 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 }; @@ -144,6 +157,12 @@ struct static_constant_interface : public leaf 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, 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 index 6042b1834e..b586605eed 100644 --- 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 @@ -32,10 +32,10 @@ using det2x2 = difference >; using orient2d = det2x2 - < - difference <_1, _5>, difference<_2, _6>, - difference <_3, _5>, difference<_4, _6> - >; + < + difference <_1, _5>, difference<_2, _6>, + difference <_3, _5>, difference<_4, _6> + >; template < 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 index 59e4d3f47b..9a9de3d1ff 100644 --- 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 @@ -20,6 +20,13 @@ 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 { 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 index 26ca2884fa..43095caab6 100644 --- 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 @@ -27,6 +27,22 @@ 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, @@ -55,7 +71,7 @@ struct semi_static_filter using ct = CalculationType; public: template - static inline int apply(const Reals&... args) + static inline int apply(Reals const&... args) { std::array input {{ static_cast(args)... }}; 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 index 96725db461..9fafd9f13b 100644 --- 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 @@ -31,6 +31,9 @@ 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 struct stage_a_error_bound_impl { 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 index d93eb6497f..6f88418a7b 100644 --- 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 @@ -27,6 +27,16 @@ 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, @@ -50,7 +60,7 @@ class static_filter ct error_bound() const { return m_error_bound; } template - inline static_filter(const Reals&... args) + inline static_filter(Reals const&... args) : m_error_bound(approximate_value( std::array {static_cast(args)...})) @@ -60,7 +70,7 @@ class static_filter } template - inline int apply(const Reals&... args) const + inline int apply(Reals const&... args) const { std::array input {static_cast(args)...}; std::array::value> results; @@ -85,7 +95,7 @@ class static_filter } template - inline int operator()(const Reals&... args) const + inline int operator()(Reals const&... args) const { return apply(args...); } From 4da968eed3134feb26368ad051dbc467f2da5ce8 Mon Sep 17 00:00:00 2001 From: Tinko Bartels Date: Thu, 27 Aug 2020 01:17:04 +0200 Subject: [PATCH 3/6] Improvements to type passing. --- .../generic_robust_predicates/approximate.cpp | 11 +- .../cartesian/detail/almost_static_filter.hpp | 1 + .../cartesian/detail/approximate.hpp | 354 ++++++++---------- .../cartesian/detail/expression_tree.hpp | 15 + .../cartesian/detail/result_propagation.hpp | 91 +++++ .../cartesian/detail/semi_static_filter.hpp | 12 +- .../cartesian/detail/static_filter.hpp | 29 +- 7 files changed, 304 insertions(+), 209 deletions(-) create mode 100644 include/boost/geometry/extensions/generic_robust_predicates/strategies/cartesian/detail/result_propagation.hpp diff --git a/extensions/test/generic_robust_predicates/approximate.cpp b/extensions/test/generic_robust_predicates/approximate.cpp index a7a4b76487..4a12d3ecae 100644 --- a/extensions/test/generic_robust_predicates/approximate.cpp +++ b/extensions/test/generic_robust_predicates/approximate.cpp @@ -33,6 +33,7 @@ void test_all() using bg::detail::generic_robust_predicates::approximate_value; using bg::detail::generic_robust_predicates::get_approx; using bg::detail::generic_robust_predicates::approximate_interim; + using bg::detail::generic_robust_predicates::argument_list; using ct = CalculationType; ct r1 = approximate_value, ct>( std::array{1.0, 2.0}); @@ -47,14 +48,16 @@ void test_all() difference<_3, _4> >; using evals = post_order; + using arg_list_input = argument_list<4>; + using arg_list = boost::mp11::mp_list; std::array::value> r; std::array input {5.0, 3.0, 2.0, 8.0}; - approximate_interim(r, input); - ct r3 = get_approx(r, input); + approximate_interim(r, input); + ct r3 = get_approx(r, input); BOOST_CHECK_EQUAL(2.0, r3); - ct r4 = get_approx(r, input); + ct r4 = get_approx(r, input); BOOST_CHECK_EQUAL(-6.0, r4); - ct r5 = get_approx(r, input); + ct r5 = get_approx(r, input); BOOST_CHECK_EQUAL(-12.0, r5); } 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 index 7727c1c4bc..20d7e968e5 100644 --- 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 @@ -75,6 +75,7 @@ class almost_static_filter extrema_array m_extrema; StaticFilter m_filter; public: + static constexpr std::size_t arg_count = StaticFilter::arg_count; const StaticFilter& filter() const { return m_filter; } inline almost_static_filter() { diff --git a/include/boost/geometry/extensions/generic_robust_predicates/strategies/cartesian/detail/approximate.hpp b/include/boost/geometry/extensions/generic_robust_predicates/strategies/cartesian/detail/approximate.hpp index 53dabca59e..807d6a0e6d 100644 --- a/include/boost/geometry/extensions/generic_robust_predicates/strategies/cartesian/detail/approximate.hpp +++ b/include/boost/geometry/extensions/generic_robust_predicates/strategies/cartesian/detail/approximate.hpp @@ -15,12 +15,14 @@ #include #include #include +#include #include #include #include #include +#include namespace boost { namespace geometry { @@ -35,415 +37,381 @@ namespace detail { namespace generic_robust_predicates template < typename Node, - std::size_t N, + typename InputList, typename Real, - typename InputArr -> -struct get_nth_real_impl -{ - static inline Real apply(InputArr const& input) - { - return input[N - 1]; - } -}; - -template -< - typename Node, - typename Real, - typename InputArr -> -struct get_nth_real_impl -{ - static inline Real apply(InputArr const&) - { - return Node::value; - } -}; - -template -< - typename Node, - std::size_t N, - typename Real, - typename InputArr -> -inline Real get_nth_real(InputArr const& input) -{ - return get_nth_real_impl::apply(input); -} - -template -< - typename All, - typename Node, - typename Real, - typename Arr, - bool is_leaf, - typename InputArr + bool IsStaticConstant, + typename ...InputArr > struct get_approx_impl { - static inline Real apply(Arr& interim_results, InputArr const&) + static inline Real apply(const InputArr&... inputs) { - return interim_results[boost::mp11::mp_find::value]; + using indices = index_pair; + return std::get( + std::get( + std::forward_as_tuple(inputs...) + ) + ); } }; template < - typename All, typename Node, + typename InputList, typename Real, - typename Arr, - typename InputArr + typename ...InputArr > -struct get_approx_impl +struct get_approx_impl + < + Node, + InputList, + Real, + true, + InputArr... + > { - static inline Real apply(Arr&, InputArr const& input) + static inline Real apply(const InputArr&...) { - return get_nth_real(input); + return Node::value; } }; template < - typename All, typename Node, + typename InputList, typename Real, - typename Arr, - typename InputArr + typename ...InputArr > -inline Real get_approx(Arr& interim_results, InputArr const& input) +inline Real get_approx(const InputArr&... inputs) { return get_approx_impl < - All, Node, + InputList, Real, - Arr, - is_leaf::value, - InputArr - >::apply(interim_results, input); + is_static_constant::value, + InputArr... + >::apply(inputs...); } template < - typename All, typename Remaining, + typename InputList, typename Real, - typename Arr, operator_types Op, - typename InputArr + typename ...Arr > struct approximate_interim_impl {}; template < - typename All, typename Remaining, + typename InputList, typename Real, - typename Arr, bool Empty, - typename InputArr + typename ...Arr > struct approximate_remainder_impl { - static inline void apply(Arr& interim_results, InputArr const& input) + static inline void apply(Arr&... inputs) { using node = boost::mp11::mp_front; approximate_interim_impl < - All, Remaining, + InputList, Real, - Arr, node::operator_type, - InputArr - >::apply(interim_results, input); + Arr... + >::apply(inputs...); } }; template < - typename All, typename Remaining, + typename InputList, typename Real, - typename Arr, - typename InputArr + typename ...Arr > struct approximate_remainder_impl < - All, Remaining, + InputList, Real, - Arr, true, - InputArr + Arr... > { - static inline void apply(Arr&, InputArr const&) {} + static inline void apply(Arr&...) {} }; template < - typename All, typename Remaining, + typename InputList, typename Real, - typename Arr, - typename InputArr + typename ...Arr > -inline void approximate_remainder(Arr& interim_results, InputArr const& input) +inline void approximate_remainder(Arr&... inputs) { approximate_remainder_impl < - All, Remaining, + InputList, Real, - Arr, boost::mp11::mp_empty::value, - InputArr - >::apply(interim_results, input); + Arr... + >::apply(inputs...); } template < - typename All, typename Remaining, + typename InputList, typename Real, - typename Arr, - typename InputArr + typename ...Arr > struct approximate_interim_impl < - All, Remaining, + InputList, Real, - Arr, operator_types::product, - InputArr + Arr... > { - static inline void apply(Arr& interim_results, InputArr const& input) + static inline void apply(Arr&... inputs) { using node = boost::mp11::mp_front; - interim_results[boost::mp11::mp_find::value] = - get_approx(interim_results, - input) - * get_approx(interim_results, - input); + using indices = index_pair; + Real& out = std::get( + std::get( + std::forward_as_tuple(inputs...) + ) + ); + out = get_approx(inputs...) + * get_approx(inputs...); approximate_remainder < - All, boost::mp11::mp_pop_front, + InputList, Real - >(interim_results, input); + >(inputs...); } }; template < - typename All, typename Remaining, + typename InputList, typename Real, - typename Arr, - typename InputArr + typename ...Arr > struct approximate_interim_impl < - All, Remaining, + InputList, Real, - Arr, operator_types::max, - InputArr + Arr... > { - static inline void apply(Arr& interim_results, InputArr const& input) + static inline void apply(Arr&... inputs) { using node = boost::mp11::mp_front; - interim_results[boost::mp11::mp_find::value] = std::max( - get_approx(interim_results, - input), - get_approx(interim_results, - input)); + using indices = index_pair; + Real& out = std::get( + std::get( + std::forward_as_tuple(inputs...) + ) + ); + out = std::max( + get_approx(inputs...), + get_approx(inputs...) + ); approximate_remainder < - All, boost::mp11::mp_pop_front, + InputList, Real - >(interim_results, input); + >(inputs...); } }; template < - typename All, typename Remaining, + typename InputList, typename Real, - typename Arr, - typename InputArr + typename ...Arr > struct approximate_interim_impl < - All, Remaining, + InputList, Real, - Arr, operator_types::min, - InputArr + Arr... > { - static inline void apply(Arr& interim_results, InputArr const& input) + static inline void apply(Arr&... inputs) { using node = boost::mp11::mp_front; - interim_results[boost::mp11::mp_find::value] = std::min( - get_approx(interim_results, - input), - get_approx(interim_results, - input)); + using indices = index_pair; + Real& out = std::get( + std::get( + std::forward_as_tuple(inputs...) + ) + ); + out = std::min( + get_approx(inputs...), + get_approx(inputs...) + ); approximate_remainder < - All, boost::mp11::mp_pop_front, + InputList, Real - >(interim_results, input); + >(inputs...); } }; template < - typename All, typename Remaining, + typename InputList, typename Real, - typename Arr, - typename InputArr + typename ...Arr > struct approximate_interim_impl < - All, Remaining, + InputList, Real, - Arr, operator_types::sum, - InputArr + Arr... > { - static inline void apply(Arr& interim_results, InputArr const& input) + static inline void apply(Arr&... inputs) { using node = boost::mp11::mp_front; - interim_results[boost::mp11::mp_find::value] = - get_approx(interim_results, - input) - + get_approx(interim_results, - input); + using indices = index_pair; + Real& out = std::get( + std::get( + std::forward_as_tuple(inputs...) + ) + ); + out = get_approx(inputs...) + + get_approx(inputs...); approximate_remainder < - All, boost::mp11::mp_pop_front, + InputList, Real - >(interim_results, input); + >(inputs...); } }; template < - typename All, typename Remaining, + typename InputList, typename Real, - typename Arr, - typename InputArr + typename ...Arr > struct approximate_interim_impl < - All, Remaining, + InputList, Real, - Arr, operator_types::difference, - InputArr + Arr... > { - static inline void apply(Arr& interim_results, InputArr const& input) + static inline void apply(Arr&... inputs) { using node = boost::mp11::mp_front; - interim_results[boost::mp11::mp_find::value] = - get_approx(interim_results, - input) - - get_approx(interim_results, - input); + using indices = index_pair; + Real& out = std::get( + std::get( + std::forward_as_tuple(inputs...) + ) + ); + out = get_approx(inputs...) + - get_approx(inputs...); approximate_remainder < - All, boost::mp11::mp_pop_front, + InputList, Real - >(interim_results, input); + >(inputs...); } }; template < - typename All, typename Remaining, + typename InputList, typename Real, - typename Arr, - typename InputArr + typename ...Arr > struct approximate_interim_impl < - All, Remaining, + InputList, Real, - Arr, operator_types::abs, - InputArr + Arr... > { - static inline void apply(Arr& interim_results, InputArr const& input) + static inline void apply(Arr&... inputs) { using node = boost::mp11::mp_front; - interim_results[boost::mp11::mp_find::value] = - std::abs(get_approx - < - All, - typename node::child, - Real - >(interim_results, input)); + using indices = index_pair; + Real& out = std::get( + std::get( + std::forward_as_tuple(inputs...) + ) + ); + out = std::abs( + get_approx(inputs...)); approximate_remainder < - All, boost::mp11::mp_pop_front, + InputList, Real - >(interim_results, input); + >(inputs...); } }; template < - typename All, typename Remaining, + typename InputList, typename Real, - typename Arr, - typename InputArr + typename ...Arr > struct approximate_interim_impl < - All, Remaining, + InputList, Real, - Arr, operator_types::no_op, - InputArr + Arr... > { - static inline void apply(Arr& interim_results, InputArr const& input) + static inline void apply(Arr&... inputs) { approximate_remainder < - All, boost::mp11::mp_pop_front, + InputList, Real - >(interim_results, input); + >(inputs...); } }; @@ -456,29 +424,33 @@ struct approximate_interim_impl template < - typename All, typename Remaining, + typename InputList, typename Real, - typename Arr, - typename InputArr + typename ...Arr > -inline void approximate_interim(Arr& interim_results, InputArr const& input) +inline void approximate_interim(Arr&... inputs) { approximate_remainder < - All, Remaining, + InputList, Real - >(interim_results, input); + >(inputs...); } template -inline Real approximate_value(InputArr const& input) +inline Real approximate_value(const InputArr& input) { using stack = typename boost::mp11::mp_unique>; using evals = typename boost::mp11::mp_remove_if; + using arg_list = boost::mp11::mp_list + < + evals, + argument_list::value> + >; std::array::value> results; - approximate_interim(results, input); + approximate_interim(results, input); return results.back(); } 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 index 91bbd6d3d8..dbdaff547e 100644 --- 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 @@ -249,6 +249,21 @@ struct max_argn_impl using type = boost::mp11::mp_size_t; }; +template > +struct is_static_constant_impl +{ + using type = boost::mp11::mp_false; +}; + +template +struct is_static_constant_impl +{ + using type = boost::mp11::mp_bool; +}; + +template +using is_static_constant = typename is_static_constant_impl::type; + using _1 = argument<1>; using _2 = argument<2>; using _3 = argument<3>; diff --git a/include/boost/geometry/extensions/generic_robust_predicates/strategies/cartesian/detail/result_propagation.hpp b/include/boost/geometry/extensions/generic_robust_predicates/strategies/cartesian/detail/result_propagation.hpp new file mode 100644 index 0000000000..4e25c0be96 --- /dev/null +++ b/include/boost/geometry/extensions/generic_robust_predicates/strategies/cartesian/detail/result_propagation.hpp @@ -0,0 +1,91 @@ +// 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_RESULT_PROPAGATION_HPP +#define BOOST_GEOMETRY_EXTENSIONS_GENERIC_ROBUST_PREDICATES_STRATEGIES_CARTESIAN_DETAIL_RESULT_PROPAGATION_HPP + +#include +#include +#include + +#include + +namespace boost { namespace geometry +{ + +namespace detail { namespace generic_robust_predicates +{ + +template +< + std::size_t Last, + typename In = boost::mp11::mp_list<>, + std::size_t At = 1 +> +struct argument_list_impl +{ +private: + using next = boost::mp11::mp_push_back>; +public: + using type = typename argument_list_impl::type; +}; + +template +< + std::size_t Last, + typename In +> +struct argument_list_impl + < + Last, + In, + Last + > +{ + using type = boost::mp11::mp_push_back>; +}; + +template using argument_list = + typename argument_list_impl::type; + +template +< + typename Needle, + typename Haystacks +> +struct index_pair_impl +{ +private: + template + using contains = typename boost::mp11::mp_bind_back + < + boost::mp11::mp_contains, + Needle + >::template fn; + using outer = boost::mp11::mp_find_if; + using inner_list = boost::mp11::mp_at; + using inner = boost::mp11::mp_find; +public: + using type = std::pair; +}; + +template +< + typename Needle, + typename Haystacks +> +using index_pair = typename index_pair_impl::type; + +}} // namespace detail::generic_robust_predicates + +}} // namespace boost::geometry + +#endif // BOOST_GEOMETRY_EXTENSIONS_GENERIC_ROBUST_PREDICATES_STRATEGIES_CARTESIAN_DETAIL_RESULT_PROPAGATION_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 index 43095caab6..a763be9e1c 100644 --- 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 @@ -20,6 +20,7 @@ #include #include +#include namespace boost { namespace geometry { @@ -71,16 +72,17 @@ struct semi_static_filter using ct = CalculationType; public: template - static inline int apply(Reals const&... args) + static inline int apply(const Reals&... args) { + using arg_list_input = argument_list; + using arg_list = boost::mp11::mp_list; std::array input {{ static_cast(args)... }}; std::array::value> results; - approximate_interim(results, input); + approximate_interim(results, input); const ct error_bound = - get_approx(results, input); - const ct det = - get_approx(results, input); + get_approx(results, input); + const ct det = get_approx(results, input); if (det > error_bound) { return 1; 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 index 6f88418a7b..5acdbd9a5d 100644 --- 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 @@ -16,10 +16,10 @@ #include #include -#include -#include -#include +#include +#include +#include namespace boost { namespace geometry { @@ -54,28 +54,39 @@ class static_filter post_order >; ct m_error_bound; + static constexpr std::size_t extrema_count = + max_argn::value; public: + static constexpr std::size_t arg_count = max_argn::value; + inline static_filter() {} ct error_bound() const { return m_error_bound; } template - inline static_filter(Reals const&... args) + inline static_filter(const Reals&... args) : m_error_bound(approximate_value( - std::array + std::array {static_cast(args)...})) { - static_assert(sizeof...(Reals) == max_argn::value, + static_assert(sizeof...(Reals) == extrema_count, "Number of constructor arguments is incompatible with error expression."); } + inline static_filter(const std::array& extrema) + : m_error_bound(approximate_value(extrema)) {} + + inline static_filter(const ct& b) : m_error_bound(b) {} + template - inline int apply(Reals const&... args) const + inline int apply(const Reals&... args) const { + using arg_list_input = argument_list; + using arg_list = boost::mp11::mp_list; std::array input {static_cast(args)...}; std::array::value> results; - approximate_interim(results, input); - const ct det = get_approx(results, input); + approximate_interim(results, input); + const ct det = get_approx(results, input); if (det > m_error_bound) { return 1; From 62ca59208ca5f157250971489120179c49760acf Mon Sep 17 00:00:00 2001 From: Tinko Bartels Date: Sun, 30 Aug 2020 19:04:33 +0200 Subject: [PATCH 4/6] Correction of arithmetic mistake and naming consistency improvement. --- .../strategies/cartesian/detail/expressions.hpp | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) 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 index b586605eed..723cd114d9 100644 --- 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 @@ -43,14 +43,14 @@ template typename A21, typename A22, typename A23, typename A31, typename A32, typename A33 > -struct det3x3_helper +struct det3x3_impl { private: using minor1 = product>; using minor2 = product>; using minor3 = product>; public: - using type = sum>; + using type = sum, minor3>; }; template @@ -59,7 +59,7 @@ template typename A21, typename A22, typename A23, typename A31, typename A32, typename A33 > -using det3x3 = typename det3x3_helper +using det3x3 = typename det3x3_impl < A11, A12, A13, A21, A22, A23, @@ -73,7 +73,7 @@ using orient3d = det3x3 difference<_7, _10>, difference<_8, _11>, difference<_9, _12> >; -struct incircle_helper +struct incircle_impl { private: using adx = difference<_1, _7>; @@ -94,7 +94,7 @@ struct incircle_helper >; }; -using incircle = incircle_helper::type; +using incircle = incircle_impl::type; }} // namespace detail::generic_robust_predicates From 543f2d4b7ee181c142a58adc712db8be3fc57738 Mon Sep 17 00:00:00 2001 From: Tinko Bartels Date: Thu, 4 Mar 2021 21:04:07 +0100 Subject: [PATCH 5/6] Adds Stage B filter and Staged Predicate, major refactoring of Stage A and Stage D --- .../cartesian/detail/almost_static_filter.hpp | 73 +- .../cartesian/detail/approximate.hpp | 461 ---- .../cartesian/detail/coefficient_list.hpp | 503 ---- .../cartesian/detail/error_bound.hpp | 623 ----- .../cartesian/detail/expansion_arithmetic.hpp | 2125 +++++++++++++++++ .../cartesian/detail/expansion_eval.hpp | 853 +++++++ .../cartesian/detail/expression_eval.hpp | 281 +++ .../cartesian/detail/expression_tree.hpp | 64 +- .../cartesian/detail/expressions.hpp | 512 +++- .../cartesian/detail/interval_error_bound.hpp | 156 +- .../cartesian/detail/result_propagation.hpp | 91 - .../cartesian/detail/semi_static_filter.hpp | 35 +- .../strategies/cartesian/detail/stage_a.hpp | 115 +- .../cartesian/detail/stage_a_error_bound.hpp | 235 ++ .../strategies/cartesian/detail/stage_b.hpp | 253 ++ .../strategies/cartesian/detail/stage_d.hpp | 116 + .../cartesian/detail/staged_predicate.hpp | 238 ++ .../cartesian/detail/static_filter.hpp | 50 +- .../cartesian/detail/static_util.hpp | 130 - 19 files changed, 4769 insertions(+), 2145 deletions(-) delete mode 100644 include/boost/geometry/extensions/generic_robust_predicates/strategies/cartesian/detail/approximate.hpp delete mode 100644 include/boost/geometry/extensions/generic_robust_predicates/strategies/cartesian/detail/coefficient_list.hpp delete mode 100644 include/boost/geometry/extensions/generic_robust_predicates/strategies/cartesian/detail/error_bound.hpp create mode 100644 include/boost/geometry/extensions/generic_robust_predicates/strategies/cartesian/detail/expansion_arithmetic.hpp create mode 100644 include/boost/geometry/extensions/generic_robust_predicates/strategies/cartesian/detail/expansion_eval.hpp create mode 100644 include/boost/geometry/extensions/generic_robust_predicates/strategies/cartesian/detail/expression_eval.hpp delete mode 100644 include/boost/geometry/extensions/generic_robust_predicates/strategies/cartesian/detail/result_propagation.hpp create mode 100644 include/boost/geometry/extensions/generic_robust_predicates/strategies/cartesian/detail/stage_a_error_bound.hpp create mode 100644 include/boost/geometry/extensions/generic_robust_predicates/strategies/cartesian/detail/stage_b.hpp create mode 100644 include/boost/geometry/extensions/generic_robust_predicates/strategies/cartesian/detail/stage_d.hpp create mode 100644 include/boost/geometry/extensions/generic_robust_predicates/strategies/cartesian/detail/staged_predicate.hpp delete mode 100644 include/boost/geometry/extensions/generic_robust_predicates/strategies/cartesian/detail/static_util.hpp 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 index 20d7e968e5..6c2e9bdd92 100644 --- 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 @@ -24,28 +24,6 @@ namespace boost { namespace geometry namespace detail { namespace generic_robust_predicates { -template -struct make_filter_impl -{ - template - static Filter apply(ExtremaArray const& extrema, Reals const&... args) - { - return make_filter_impl - ::apply(extrema, args..., extrema[N]); - } -}; - -template -struct make_filter_impl -{ - template - static Filter apply(ExtremaArray const& extrema, Reals const&... args) - { - Filter f(args...); - return f; - } -}; - //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: // @@ -69,14 +47,15 @@ template class almost_static_filter { private: - using expression_max_argn = max_argn; using ct = CalculationType; - using extrema_array = std::array; + using extrema_array = std::array>; extrema_array m_extrema; StaticFilter m_filter; public: - static constexpr std::size_t arg_count = StaticFilter::arg_count; - const StaticFilter& filter() const { return m_filter; } + 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(), @@ -86,42 +65,42 @@ class almost_static_filter m_extrema.end(), std::numeric_limits::infinity()); } - template - int apply(Reals const&... args) const + template + int apply(CTs const&... args) const { return m_filter.apply(args...); } - template - inline void update_extrema(Reals const&... args) + template + inline void update_extrema(CTs const&... args) { - std::array input {{ static_cast(args)... }}; - for(int i = 0; i < m_extrema.size() / 2; ++i) + 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(int i = m_extrema.size() / 2; i < m_extrema.size(); ++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(Reals const&... args) + template + inline bool update_extrema_check(CTs const&... args) { bool changed = false; - std::array input {{ static_cast(args)... }}; - for(int i = 0; i < m_extrema.size() / 2; ++i) + std::array input {{ static_cast(args)... }}; + for(unsigned int i = 0; i < m_extrema.size() / 2; ++i) { - if(input[i] > m_extrema[i]) + if (input[i] > m_extrema[i]) { changed = true; m_extrema[i] = input[i]; } } - for(int i = m_extrema.size() / 2; i < m_extrema.size(); ++i) + for(unsigned int i = m_extrema.size() / 2; i < m_extrema.size(); ++i) { - if(input[i] < m_extrema[i]) + if (input[i] < m_extrema[i]) { changed = true; m_extrema[i] = input[i]; @@ -130,14 +109,16 @@ class almost_static_filter return changed; } + template + inline void update(CTs const&... args) + { + update_extrema(args...); + update_filter(); + } + inline void update_filter() { - m_filter = make_filter_impl - < - StaticFilter, - 0, - 2 * expression_max_argn::value - >::apply(m_extrema); + m_filter = StaticFilter(m_extrema); } }; diff --git a/include/boost/geometry/extensions/generic_robust_predicates/strategies/cartesian/detail/approximate.hpp b/include/boost/geometry/extensions/generic_robust_predicates/strategies/cartesian/detail/approximate.hpp deleted file mode 100644 index 807d6a0e6d..0000000000 --- a/include/boost/geometry/extensions/generic_robust_predicates/strategies/cartesian/detail/approximate.hpp +++ /dev/null @@ -1,461 +0,0 @@ -// 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_APPROXIMATE_HPP -#define BOOST_GEOMETRY_EXTENSIONS_GENERIC_ROBUST_PREDICATES_STRATEGIES_CARTESIAN_DETAIL_APPROXIMATE_HPP - -#include -#include -#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. -//The most important template in this file is approximate_interim. - -template -< - typename Node, - typename InputList, - typename Real, - bool IsStaticConstant, - typename ...InputArr -> -struct get_approx_impl -{ - static inline Real apply(const InputArr&... inputs) - { - using indices = index_pair; - return std::get( - std::get( - std::forward_as_tuple(inputs...) - ) - ); - } -}; - -template -< - typename Node, - typename InputList, - typename Real, - typename ...InputArr -> -struct get_approx_impl - < - Node, - InputList, - Real, - true, - InputArr... - > -{ - static inline Real apply(const InputArr&...) - { - return Node::value; - } -}; - -template -< - typename Node, - typename InputList, - typename Real, - typename ...InputArr -> -inline Real get_approx(const InputArr&... inputs) -{ - return get_approx_impl - < - Node, - InputList, - Real, - is_static_constant::value, - InputArr... - >::apply(inputs...); -} - -template -< - typename Remaining, - typename InputList, - typename Real, - operator_types Op, - typename ...Arr -> -struct approximate_interim_impl {}; - -template -< - typename Remaining, - typename InputList, - typename Real, - bool Empty, - typename ...Arr -> -struct approximate_remainder_impl -{ - static inline void apply(Arr&... inputs) - { - using node = boost::mp11::mp_front; - approximate_interim_impl - < - Remaining, - InputList, - Real, - node::operator_type, - Arr... - >::apply(inputs...); - } -}; - -template -< - typename Remaining, - typename InputList, - typename Real, - typename ...Arr -> -struct approximate_remainder_impl - < - Remaining, - InputList, - Real, - true, - Arr... - > -{ - static inline void apply(Arr&...) {} -}; - -template -< - typename Remaining, - typename InputList, - typename Real, - typename ...Arr -> -inline void approximate_remainder(Arr&... inputs) -{ - approximate_remainder_impl - < - Remaining, - InputList, - Real, - boost::mp11::mp_empty::value, - Arr... - >::apply(inputs...); -} - -template -< - typename Remaining, - typename InputList, - typename Real, - typename ...Arr -> -struct approximate_interim_impl - < - Remaining, - InputList, - Real, - operator_types::product, - Arr... - > -{ - static inline void apply(Arr&... inputs) - { - using node = boost::mp11::mp_front; - using indices = index_pair; - Real& out = std::get( - std::get( - std::forward_as_tuple(inputs...) - ) - ); - out = get_approx(inputs...) - * get_approx(inputs...); - approximate_remainder - < - boost::mp11::mp_pop_front, - InputList, - Real - >(inputs...); - } -}; - -template -< - typename Remaining, - typename InputList, - typename Real, - typename ...Arr -> -struct approximate_interim_impl - < - Remaining, - InputList, - Real, - operator_types::max, - Arr... - > -{ - static inline void apply(Arr&... inputs) - { - using node = boost::mp11::mp_front; - using indices = index_pair; - Real& out = std::get( - std::get( - std::forward_as_tuple(inputs...) - ) - ); - out = std::max( - get_approx(inputs...), - get_approx(inputs...) - ); - approximate_remainder - < - boost::mp11::mp_pop_front, - InputList, - Real - >(inputs...); - } -}; - -template -< - typename Remaining, - typename InputList, - typename Real, - typename ...Arr -> -struct approximate_interim_impl - < - Remaining, - InputList, - Real, - operator_types::min, - Arr... - > -{ - static inline void apply(Arr&... inputs) - { - using node = boost::mp11::mp_front; - using indices = index_pair; - Real& out = std::get( - std::get( - std::forward_as_tuple(inputs...) - ) - ); - out = std::min( - get_approx(inputs...), - get_approx(inputs...) - ); - approximate_remainder - < - boost::mp11::mp_pop_front, - InputList, - Real - >(inputs...); - } -}; - -template -< - typename Remaining, - typename InputList, - typename Real, - typename ...Arr -> -struct approximate_interim_impl - < - Remaining, - InputList, - Real, - operator_types::sum, - Arr... - > -{ - static inline void apply(Arr&... inputs) - { - using node = boost::mp11::mp_front; - using indices = index_pair; - Real& out = std::get( - std::get( - std::forward_as_tuple(inputs...) - ) - ); - out = get_approx(inputs...) - + get_approx(inputs...); - approximate_remainder - < - boost::mp11::mp_pop_front, - InputList, - Real - >(inputs...); - } -}; - -template -< - typename Remaining, - typename InputList, - typename Real, - typename ...Arr -> -struct approximate_interim_impl - < - Remaining, - InputList, - Real, - operator_types::difference, - Arr... - > -{ - static inline void apply(Arr&... inputs) - { - using node = boost::mp11::mp_front; - using indices = index_pair; - Real& out = std::get( - std::get( - std::forward_as_tuple(inputs...) - ) - ); - out = get_approx(inputs...) - - get_approx(inputs...); - approximate_remainder - < - boost::mp11::mp_pop_front, - InputList, - Real - >(inputs...); - } -}; - -template -< - typename Remaining, - typename InputList, - typename Real, - typename ...Arr -> -struct approximate_interim_impl - < - Remaining, - InputList, - Real, - operator_types::abs, - Arr... - > -{ - static inline void apply(Arr&... inputs) - { - using node = boost::mp11::mp_front; - using indices = index_pair; - Real& out = std::get( - std::get( - std::forward_as_tuple(inputs...) - ) - ); - out = std::abs( - get_approx(inputs...)); - approximate_remainder - < - boost::mp11::mp_pop_front, - InputList, - Real - >(inputs...); - } -}; - -template -< - typename Remaining, - typename InputList, - typename Real, - typename ...Arr -> -struct approximate_interim_impl - < - Remaining, - InputList, - Real, - operator_types::no_op, - Arr... - > -{ - static inline void apply(Arr&... inputs) - { - approximate_remainder - < - boost::mp11::mp_pop_front, - InputList, - Real - >(inputs...); - } -}; - -//All expects an boost::mp11::mp_list of all expressions that need to be -//evaluated. Remaining expects an boost::mp11::mp_list of the expressions -//that are left to be evaluated. In the first call it is expected to be -//equal to All and it serves as an anchor for template recursion. Real is a -//floating-point type. The remaining template arguments are deduced from -//parameters. - -template -< - typename Remaining, - typename InputList, - typename Real, - typename ...Arr -> -inline void approximate_interim(Arr&... inputs) -{ - approximate_remainder - < - Remaining, - InputList, - Real - >(inputs...); -} - -template -inline Real approximate_value(const InputArr& input) -{ - using stack = typename boost::mp11::mp_unique>; - using evals = typename boost::mp11::mp_remove_if; - using arg_list = boost::mp11::mp_list - < - evals, - argument_list::value> - >; - std::array::value> results; - approximate_interim(results, input); - return results.back(); -} - -}} // namespace detail::generic_robust_predicates - -}} // namespace boost::geometry - -#endif // BOOST_GEOMETRY_EXTENSIONS_GENERIC_ROBUST_PREDICATES_STRATEGIES_CARTESIAN_DETAIL_APPROXIMATE_HPP diff --git a/include/boost/geometry/extensions/generic_robust_predicates/strategies/cartesian/detail/coefficient_list.hpp b/include/boost/geometry/extensions/generic_robust_predicates/strategies/cartesian/detail/coefficient_list.hpp deleted file mode 100644 index 331f983e2f..0000000000 --- a/include/boost/geometry/extensions/generic_robust_predicates/strategies/cartesian/detail/coefficient_list.hpp +++ /dev/null @@ -1,503 +0,0 @@ -// 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_COEFFICIENT_LIST_HPP -#define BOOST_GEOMETRY_EXTENSIONS_GENERIC_ROBUST_PREDICATES_STRATEGIES_CARTESIAN_DETAIL_COEFFICIENT_LIST_HPP - -#include - -#include -#include - -#include - -namespace boost { namespace geometry -{ - -namespace detail { namespace generic_robust_predicates -{ - -//The templates in this file are helpers for the computation of error maps and -//error expressions (see also error_bound.hpp). A coefficient list is a -//boost::mp11::mp_list that contains integers that represent the coefficients -//of an polynomial. This is used for the polynomials in epsilon that occur in -//forward error analysis, e.g. epsilon + 2 * epsilon^2 + epsilon^3 would be -//represented by the type mp_list, mp_int<2>, mp_int<1>>. -// -//See p. 39 in "Adaptive Precision Floating-Point Arithmetic and Fast Robust -//Geometric Predicates" by Jonathan Richard Shewchuk (can be downloaded from -//https://www.cs.cmu.edu/~quake/robust.html) for the inspiration for the -//methods implemented below. - -//only meant to be used in assertions -template using is_mp_int = - boost::mp11::mp_same>; - -template using is_coeff_list = boost::mp11::mp_and - < - boost::mp11::mp_all_of, - boost::mp11::mp_similar> - >; - -template -struct coeff_truncate_impl -{ -private: - static_assert(is_coeff_list::value, "L must be a coefficient list"); - using first_nz = boost::mp11::mp_find_if; - using after_second_nz = - boost::mp11::mp_min< - boost::mp11::mp_size_t, - boost::mp11::mp_size - >; - using tail = boost::mp11::mp_erase_c; - using head = boost::mp11::mp_erase - < - L, - after_second_nz, - boost::mp11::mp_size - >; - using round_up = boost::mp11::mp_less - < - boost::mp11::mp_find_if, - boost::mp11::mp_size - >; -public: - using type = boost::mp11::mp_if - < - round_up, - boost::mp11::mp_push_back>, - head - >; - static_assert(is_coeff_list::value, "type should be a coefficient list"); -}; - -//Because of the way floating-point numbers are rounded, for a given epsilon -//polynomial we only care about the first two non-zero coefficients. -//e.g. if a*eps + b*eps^2 + c*eps^3 is considered then c * eps^3 is always -//neglible unless c is so large that it must have resulted from the analysis -//of an expression that is so complex that it cannot be feasibly processed -//by the methods implemented in this feature. -template using coeff_truncate = typename coeff_truncate_impl::type; - -template using app_zero_b = - boost::mp11::mp_push_front>; -template using app_zero_f = - boost::mp11::mp_push_back>; - -//For a coefficient list representing the polynomial p, the following template -//produces the coefficient list representing p * (1 + eps). -template using mult_by_1_p_eps = coeff_truncate - < - boost::mp11::mp_transform - < - boost::mp11::mp_plus, - app_zero_b, - app_zero_f - > - >; - -template -< - typename L, - typename N, - typename done = boost::mp11::mp_bool -> -struct mult_by_1_p_eps_pow_impl -{ -private: - using next = mult_by_1_p_eps; -public: - using type = - typename mult_by_1_p_eps_pow_impl - < - next, - boost::mp11::mp_int - >::type; -}; - -template -struct mult_by_1_p_eps_pow_impl -{ -public: - using type = L; -}; - - -//For a coefficient list representing the polynomial p, the following template -//produces the cofficient list representing p * (1 + eps)^N -template using mult_by_1_p_eps_pow = coeff_truncate - < - typename mult_by_1_p_eps_pow_impl::type - >; - -template using div_by_1_m_eps_helper = - boost::mp11::mp_partial_sum - < - L, - boost::mp11::mp_int<0>, - boost::mp11::mp_plus - >; - -//For a coefficient list representing a polynomial p in eps, the following -//template computes a coefficient list representing a polynomial q such that -//q * (1 - p) >= p. -template using div_by_1_m_eps = coeff_truncate - < - boost::mp11::mp_push_back - < - boost::mp11::mp_pop_back>, - inc>> - > - >; - -template -< - typename L1, - typename L2, - typename L, - typename L1_Empty = typename empty_or_void::type, - typename L2_Empty = typename empty_or_void::type -> -struct coeff_merge_impl {}; - -template -< - typename L1, - typename L2, - typename L -> -struct coeff_merge_impl -< - L1, - L2, - L, - boost::mp11::mp_false, - boost::mp11::mp_false -> -{ - using type = typename coeff_merge_impl - < - boost::mp11::mp_pop_front, - boost::mp11::mp_pop_front, - boost::mp11::mp_push_back - < - L, - boost::mp11::mp_plus - < - boost::mp11::mp_front, - boost::mp11::mp_front - > - > - >::type; -}; - -template -< - typename L1, - typename L2, - typename L -> -struct coeff_merge_impl - < - L1, - L2, - L, - boost::mp11::mp_true, - boost::mp11::mp_true - > -{ - using type = L; -}; - -template -< - typename L1, - typename L2, - typename L -> -struct coeff_merge_impl - < - L1, - L2, - L, - boost::mp11::mp_true, - boost::mp11::mp_false - > -{ - using type = boost::mp11::mp_append; -}; - -template -< - typename L1, - typename L2, - typename L -> -struct coeff_merge_impl - < - L1, - L2, - L, - boost::mp11::mp_false, - boost::mp11::mp_true - > -{ - using type = boost::mp11::mp_append; -}; - -template using coeff_merge = - coeff_truncate - < - typename coeff_merge_impl - < - L1, - L2, - boost::mp11::mp_list<> - >::type - >; - -template -< - typename L1, - typename L2, - typename L, - typename L1_Empty = typename empty_or_void::type, - typename L2_Empty = typename empty_or_void::type -> -struct coeff_max_impl {}; - -template -< - typename L1, - typename L2, - typename L -> -struct coeff_max_impl -{ - using type = typename coeff_max_impl - < - boost::mp11::mp_if - < - boost::mp11::mp_less - < - boost::mp11::mp_front, - boost::mp11::mp_front - >, - boost::mp11::mp_list<>, - boost::mp11::mp_pop_front - >, - boost::mp11::mp_if - < - boost::mp11::mp_less - < - boost::mp11::mp_front, - boost::mp11::mp_front - >, - boost::mp11::mp_list<>, - boost::mp11::mp_pop_front - >, - boost::mp11::mp_push_back - < - L, - boost::mp11::mp_max - < - boost::mp11::mp_front, - boost::mp11::mp_front - > - > - >::type; -}; - -template -< - typename L1, - typename L2, - typename L -> -struct coeff_max_impl - < - L1, - L2, - L, - boost::mp11::mp_true, - boost::mp11::mp_true - > -{ - using type = L; -}; - -template -< - typename L1, - typename L2, - typename L -> -struct coeff_max_impl - < - L1, - L2, - L, - boost::mp11::mp_true, - boost::mp11::mp_false - > -{ - using type = boost::mp11::mp_append; -}; - -template -< - typename L1, - typename L2, - typename L -> -struct coeff_max_impl - < - L1, - L2, - L, - boost::mp11::mp_false, - boost::mp11::mp_true - > -{ - using type = boost::mp11::mp_append; -}; - -template using coeff_max = - coeff_truncate - < - typename coeff_max_impl - < - L1, - L2, - boost::mp11::mp_list<> - >::type - >; - -template -< - typename L, - std::size_t Tail_Size = - boost::mp11::mp_size::value - - boost::mp11::mp_find_if::value -> -struct coeff_round_impl -{ -private: - using first_nz = boost::mp11::mp_find_if; - using tail = boost::mp11::mp_erase_c; - using zero_tail = boost::mp11::mp_same - < - boost::mp11::mp_find_if, - boost::mp11::mp_size - >; - using head = boost::mp11::mp_erase_c - < - L, - first_nz::value + 2, - boost::mp11::mp_size::value - >; - using minor = boost::mp11::mp_if - < - zero_tail, - boost::mp11::mp_back, - inc> - >; - using major = boost::mp11::mp_at; - using major_rounded = boost::mp11::mp_int< 1 << log_2_ceil::value >; - using minor_rounded = boost::mp11::mp_int - < - (minor::value / major_rounded::value) * major_rounded::value < minor::value ? - (minor::value / major_rounded::value + 1) * major_rounded::value - : (minor::value / major_rounded::value) * major_rounded::value - >; -public: - using type = boost::mp11::mp_push_back - < - boost::mp11::mp_pop_back, - minor_rounded - >; -}; - -template -struct coeff_round_impl { using type = L; }; - -template -struct coeff_round_impl { using type = L; }; - -//The following template rounds a coefficient list up such that the -//corresponding polynomial evaluated for epsilon is a representable -//floating-point number. -// -//e.g. 3 * eps + 12 * eps^2 + 1 * eps^3 is rounded to 3 * eps + 16 * eps^2 -template using coeff_round = typename coeff_round_impl::type; - -template using indexed_value_product = - boost::mp11::mp_list - < - boost::mp11::mp_plus - < - boost::mp11::mp_first, - boost::mp11::mp_first - >, - boost::mp11::mp_int - < - boost::mp11::mp_second::value - * boost::mp11::mp_second::value - > - >; - -template struct add_nested -{ - template using fn = - boost::mp11::mp_plus; -}; - -template using list_product_fold = - boost::mp11::mp_map_update_q - < - M, - IV, - add_nested< boost::mp11::mp_second > - >; - -template using index_less = - boost::mp11::mp_less - < - boost::mp11::mp_first, - boost::mp11::mp_first - >; - -template using list_product = - strip_index - < - boost::mp11::mp_sort - < - boost::mp11::mp_fold - < - boost::mp11::mp_product - < - indexed_value_product, - indexed, - indexed - >, - boost::mp11::mp_list<>, - list_product_fold - >, - index_less - > - >; - -}} // namespace detail::generic_robust_predicates - -}} // namespace boost::geometry - -#endif // BOOST_GEOMETRY_EXTENSIONS_GENERIC_ROBUST_PREDICATES_STRATEGIES_CARTESIAN_DETAIL_COEFFICIENT_LIST_HPP diff --git a/include/boost/geometry/extensions/generic_robust_predicates/strategies/cartesian/detail/error_bound.hpp b/include/boost/geometry/extensions/generic_robust_predicates/strategies/cartesian/detail/error_bound.hpp deleted file mode 100644 index d42a519aa9..0000000000 --- a/include/boost/geometry/extensions/generic_robust_predicates/strategies/cartesian/detail/error_bound.hpp +++ /dev/null @@ -1,623 +0,0 @@ -// 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_ERROR_BOUND_HPP -#define BOOST_GEOMETRY_EXTENSIONS_GENERIC_ROBUST_PREDICATES_STRATEGIES_CARTESIAN_DETAIL_ERROR_BOUND_HPP - -#include - -#include -#include -#include - -#include - -namespace boost { namespace geometry -{ - -namespace detail { namespace generic_robust_predicates -{ - -//The templates in this file are utilities for the generation of error -//expressions for given arithmetic expressions. This is done using forward -//error analysis methods at compile time. This analysis produces an error map -//that is folded into an error expression. -// -//An Error Map is a template fulfilling the boost::mp11 map concept. The keys K -//are expressions (as in expression trees and the values V are lists of -//integers. Each key-value pair represents the contribution of an error term to -//the total error, i.e. the error expression is something like the sum of -//abs * (V[0] * epsilon + V[1] * epsilon^2 + ...) -// -//over all key-value pairs KV in the error map. -// -//See p. 39 in "Adaptive Precision Floating-Point Arithmetic and Fast Robust -//Geometric Predicates" by Jonathan Richard Shewchuk (can be downloaded from -//https://www.cs.cmu.edu/~quake/robust.html) for the inspiration for the -//methods implemented below. - -template using second_is_coeff_list = - is_coeff_list< boost::mp11::mp_second >; - -template using is_error_map = boost::mp11::mp_and - < - boost::mp11::mp_is_map, - boost::mp11::mp_all_of, - boost::mp11::mp_similar - > - //TODO: It would be desirable to also validate that the keys are good ->; - -template -struct error_map_insert_impl -{ -private: - using key = boost::mp11::mp_front; - using value = boost::mp11::mp_second; - using other_value = typename mp_map_at_second_or_void::type; - using merged_value = coeff_merge; - using nkv = boost::mp11::mp_list; -public: - using type = boost::mp11::mp_map_replace; -}; - -template using error_map_insert = - typename error_map_insert_impl::type; - -template -struct add_fold_operator -{ -public: - template using fn = boost::mp11::mp_map_insert - < - M, - boost::mp11::mp_list - < - K, - coeff_merge - < - typename mp_map_at_second_or_void::type, - typename mp_map_at_second_or_void::type - > - > - >; -}; - -template using add_children = boost::mp11::mp_fold - < - boost::mp11::mp_set_union - < - boost::mp11::mp_map_keys, - boost::mp11::mp_map_keys - >, - boost::mp11::mp_list<>, - add_fold_operator::template fn - >; - -template -< - typename Exp, - typename EM, - typename Out, - typename LErr = typename val_or_empty_list::type, - typename RErr = typename val_or_empty_list::type, - typename skip_decompose = - boost::mp11::mp_or - < - boost::mp11::mp_not - < - boost::mp11::mp_same - < - typename Exp::error_type, - sum_error_type - > - >, - boost::mp11::mp_same - < - boost::mp11::mp_list<>, - LErr, - RErr - > - > -> -struct decompose_add_impl -{ -private: - using decomp_left = typename decompose_add_impl - < - typename Exp::left, - EM, - Out - >::type; - using decomp_right = typename decompose_add_impl - < - typename Exp::right, - EM, - decomp_left - >::type; -public: - using type = decomp_right; -}; - -template -< - typename Exp, - typename EM, - typename Out, - typename LErr, - typename RErr -> -struct decompose_add_impl - < - Exp, - EM, - Out, - LErr, - RErr, - boost::mp11::mp_true - > -{ - using type = error_map_insert - < - boost::mp11::mp_list - < - Exp, - boost::mp11::mp_list> - >, - Out - >; -}; - -template -< - typename Exp, - typename EM, - typename Out -> -using decompose_add = typename decompose_add_impl::type; - -template -< - typename Exp, - typename EM, - typename LErr = typename val_or_empty_list::type, - typename RErr = typename val_or_empty_list::type, - typename Children_Empty = - boost::mp11::mp_and - < - typename empty_or_void::type, - typename empty_or_void::type - > -> -struct sum_err_impl -{ -private: - static_assert(is_error_map::value, "LErr needs to be a valid error map."); - static_assert(is_error_map::value, "RErr needs to be a valid error map."); - using children = add_children; - static_assert( - is_error_map::value, "children needs to be a valid error map."); -public: - /* - using type = boost::mp11::mp_map_insert< - children, - boost::mp11::mp_list>> - >;*/ - using type = decompose_add; - static_assert(is_error_map::value, "type needs to be a valid error map."); -}; - -template -< - typename Exp, - typename EM, - typename LErr, - typename RErr -> -struct sum_err_impl - < - Exp, - EM, - LErr, - RErr, - boost::mp11::mp_true - > -{ - static_assert(is_error_map::value, "LErr needs to be a valid error map."); - static_assert(is_error_map::value, "RErr needs to be a valid error map."); - using type = boost::mp11::mp_list - < - boost::mp11::mp_list - < - Exp, - boost::mp11::mp_list> - > - >; - static_assert(is_error_map::value, "type needs to be a valid error map."); -}; - -template using sum_err = - typename sum_err_impl::type; - -template using pad_second = boost::mp11::mp_list - < - boost::mp11::mp_front, - boost::mp11::mp_push_front - < - boost::mp11::mp_second, - boost::mp11::mp_int<0> - > - >; - -template using pop_front_second = boost::mp11::mp_list - < - boost::mp11::mp_front, - boost::mp11::mp_pop_front> - >; - -template using increment_first_of_second = - boost::mp11::mp_transform_front; - -template using prod_entry_merge = - boost::mp11::mp_list - < - boost::mp11::mp_flatten - < - boost::mp11::mp_list - < - boost::mp11::mp_first, - boost::mp11::mp_first - > - >, - list_product - < - boost::mp11::mp_second, - boost::mp11::mp_second - > - >; - -template -< - typename Exp, - typename EM, - typename LErr = typename val_or_empty_list::type, - typename RErr = typename val_or_empty_list::type -> -struct prod_children_impl -{ -private: - static_assert(is_error_map::value, "LErr needs to be a valid error map."); - static_assert(is_error_map::value, "RErr needs to be a valid error map."); - using left = typename Exp::left; - using right = typename Exp::right; - using padded_lerr = boost::mp11::mp_transform; - using added_left_lerr = decompose_add; - using padded_rerr = boost::mp11::mp_transform; - using added_right_rerr = decompose_add; - using prod = boost::mp11::mp_product - < - prod_entry_merge, - added_left_lerr, - added_right_rerr - >; - using stripped_prod = boost::mp11::mp_transform; -public: - using type = stripped_prod; - static_assert(is_error_map::value, "type needs to be a valid error map."); -}; - -template using prod_children = - typename prod_children_impl::type; - -template -< - typename Exp, - typename EM, - typename LErr = typename val_or_empty_list::type, - typename RErr = typename val_or_empty_list::type -> -struct product_err_impl -{ -private: - static_assert(is_error_map::value, "LErr needs to be a valid error map."); - static_assert(is_error_map::value, "RErr needs to be a valid error map."); - using children = prod_children; - static_assert(is_error_map::value, "children needs to be a valid error map."); -public: - using type = boost::mp11::mp_map_insert - < - children, - boost::mp11::mp_list - < - Exp, - boost::mp11::mp_list< boost::mp11::mp_int<1> > - > - >; - static_assert(is_error_map::value, "type needs to be a valid error map."); -}; - -template using product_err = - typename product_err_impl::type; - -template -< - typename Errors, - typename Exp, - typename IsSum = boost::mp11::mp_same - < - typename Exp::error_type, - sum_error_type - > -> -struct sum_or_product_err_impl -{ - using type = sum_err; -}; - -template -< - typename Errors, - typename Exp -> -struct sum_or_product_err_impl -{ - using type = product_err; -}; - -template using sum_or_product_err = - typename sum_or_product_err_impl::type; - -template -< - typename Errors, - typename Exp -> -struct error_fold_impl -{ -private: - using lerr = typename val_or_empty_list::type; - using rerr = typename val_or_empty_list::type; - using err = sum_or_product_err; -public: - using type = boost::mp11::mp_map_insert - < - Errors, - boost::mp11::mp_list - >; -}; - - -template using error_fold = - typename error_fold_impl::type; - -template using evals_error = - boost::mp11::mp_fold - < - Evals, - boost::mp11::mp_list<>, - error_fold - >; - -template using is_mp_list = boost::mp11::mp_similar - < - boost::mp11::mp_list<>, - T - >; - -template -struct list_to_product_impl -{ -private: - using key = boost::mp11::mp_front; - using value = boost::mp11::mp_second; - using multiplications = boost::mp11::mp_int::value - 1>; - using nvalue = mult_by_1_p_eps_pow; - using nkey = boost::mp11::mp_reverse_fold< - boost::mp11::mp_pop_back, - boost::mp11::mp_back, - product - >; -public: - using type = boost::mp11::mp_list; -}; - -template using list_to_product = typename list_to_product_impl::type; - -template -< - typename M, - typename KV, - typename KeyMPList = is_mp_list< boost::mp11::mp_front > -> -struct error_map_list_to_product_fold_impl -{ - using type = error_map_insert, M>; -}; - -template -struct error_map_list_to_product_fold_impl - < - M, - KV, - boost::mp11::mp_false - > -{ - using type = error_map_insert; -}; - -template using error_map_list_to_product_fold = - typename error_map_list_to_product_fold_impl::type; - -template -struct error_map_list_to_product_impl -{ - using type = boost::mp11::mp_fold - < - M, - boost::mp11::mp_list<>, - error_map_list_to_product_fold - >; -}; - -template using error_map_list_to_product = - typename error_map_list_to_product_impl::type; - -template -struct abs_unless_non_negative -{ - using type = boost::mp11::mp_if_c - < - Expression::non_negative, - Expression, - abs - >; -}; - -template using abs_fold = - boost::mp11::mp_push_back - < - EM, - boost::mp11::mp_list - < - typename abs_unless_non_negative> - ::type, - boost::mp11::mp_second - > - >; - -template using abs_all = - boost::mp11::mp_fold - < - EM, - boost::mp11::mp_list<>, - abs_fold - >; - -template -< - typename KV1, - typename KV2 -> -struct abs_sum_error_term_impl -{ -private: - using key1 = boost::mp11::mp_front; - using key2 = boost::mp11::mp_front; - using nkey = sum; - using val1 = boost::mp11::mp_second; - using val2 = boost::mp11::mp_second; - using mval = coeff_max; - static_assert(is_coeff_list::value, "merged coefficient list fails coefficient list check"); - using nval = mult_by_1_p_eps; -public: - using type = boost::mp11::mp_list; -}; - -template using abs_sum_error_term = typename abs_sum_error_term_impl::type; - -//TODO: The following could be probably optimized for potentially produce better error bounds in some cases -// if the error map is treated as a minheap by ordering of epsilon-polynomial. -template using error_map_sum_up = - boost::mp11::mp_fold - < - boost::mp11::mp_pop_front, - boost::mp11::mp_first, - abs_sum_error_term - >; - -template::value> -struct error_map_balanced_sum_up_impl -{ -private: - static constexpr std::size_t mid = MS/2; - using left = typename error_map_balanced_sum_up_impl - < - boost::mp11::mp_erase_c - >::type; - using right = typename error_map_balanced_sum_up_impl - < - boost::mp11::mp_erase_c - >::type; -public: - using type = abs_sum_error_term; -}; - -template -struct error_map_balanced_sum_up_impl -{ - using type = boost::mp11::mp_front; -}; - -template using error_map_balanced_sum_up = - typename error_map_balanced_sum_up_impl::type; - -template -< - typename Real, - typename Exp -> -struct eps_pow -{ - static constexpr Real value = - std::numeric_limits::epsilon()/2.0 - * eps_pow>::value; -}; - -template struct eps_pow> -{ - static constexpr Real value = 1.0; -}; - -template -< - typename Real, - typename L, - typename S = boost::mp11::mp_size -> -struct eval_eps_polynomial -{ -private: - using last = boost::mp11::mp_back; - using s2last = boost::mp11::mp_back< boost::mp11::mp_pop_back >; -public: - static constexpr Real value = - s2last::value - * eps_pow>::value - + last::value * eps_pow::value; -}; - -template -< - typename Real, - typename L -> -struct eval_eps_polynomial> -{ - static constexpr Real value = - boost::mp11::mp_front::value - * std::numeric_limits::epsilon() / 2.0; -}; - -template -< - typename Real, - typename L -> -struct eval_eps_polynomial> -{ - static constexpr Real value = 0; -}; - -}} // namespace detail::generic_robust_predicates - -}} // namespace boost::geometry - -#endif // BOOST_GEOMETRY_EXTENSIONS_GENERIC_ROBUST_PREDICATES_STRATEGIES_CARTESIAN_DETAIL_ERROR_BOUND_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 index dbdaff547e..064604a501 100644 --- 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 @@ -14,6 +14,7 @@ #include #include +#include #include #include @@ -63,7 +64,8 @@ struct internal_binary_node : internal_node { using left = Left; using right = Right; - static constexpr operator_arities operator_arity = operator_arities::binary; + static constexpr operator_arities operator_arity = + operator_arities::binary; }; template @@ -76,8 +78,12 @@ struct internal_unary_node : internal_node template struct sum : public internal_binary_node { - static constexpr bool sign_exact = Left::is_leaf && Right::is_leaf; - static constexpr bool non_negative = Left::non_negative && Right::non_negative; + 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; }; @@ -223,46 +229,21 @@ template class Anchor = is_leaf> using post_order = typename post_order_impl, Anchor>::type; -template > -struct max_argn_impl; +template +constexpr std::size_t max_argn = 0; -template using max_argn = typename max_argn_impl::type; +template +constexpr std::size_t max_argn = + std::max(max_argn, + max_argn); -template -struct max_argn_impl -{ -private: - using children_list = boost::mp11::mp_rename; - using children_max_argn = - boost::mp11::mp_transform; -public: - using type = boost::mp11::mp_max_element - < - children_max_argn, - boost::mp11::mp_less - >; -}; - -template -struct max_argn_impl -{ - using type = boost::mp11::mp_size_t; -}; +template +constexpr std::size_t max_argn = + max_argn; -template > -struct is_static_constant_impl -{ - using type = boost::mp11::mp_false; -}; - -template -struct is_static_constant_impl -{ - using type = boost::mp11::mp_bool; -}; - -template -using is_static_constant = typename is_static_constant_impl::type; +template +constexpr std::size_t max_argn = + Expression::argn; using _1 = argument<1>; using _2 = argument<2>; @@ -276,6 +257,9 @@ 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 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 index 723cd114d9..b4b60cb37e 100644 --- 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 @@ -12,6 +12,11 @@ #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 @@ -20,81 +25,480 @@ namespace boost { namespace geometry namespace detail { namespace generic_robust_predicates { -template -< - typename A11, typename A12, - typename A21, typename A22 -> -using det2x2 = difference +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 < - product, - product + det, + boost::mp11::mp_transform_q + < + orient_entry, + boost::mp11::mp_iota_c< N * N > + > >; -using orient2d = det2x2 +template +using collinear = boost::mp11::mp_apply < - difference <_1, _5>, difference<_2, _6>, - difference <_3, _5>, difference<_4, _6> + det, + boost::mp11::mp_transform + < + collinear_entry, + boost::mp11::mp_iota_c + > >; -template -< - typename A11, typename A12, typename A13, - typename A21, typename A22, typename A23, - typename A31, typename A32, typename A33 -> -struct det3x3_impl +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: - using minor1 = product>; - using minor2 = product>; - using minor3 = product>; + 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: - using type = sum, minor3>; + template + using fn = + boost::mp11::mp_if_c + < + (I::value % (N + 1) != 0), + diff>, + lift + >; }; -template -< - typename A11, typename A12, typename A13, - typename A21, typename A22, typename A23, - typename A31, typename A32, typename A33 -> -using det3x3 = typename det3x3_impl +template +using innsphere = boost::mp11::mp_apply < - A11, A12, A13, - A21, A22, A23, - A31, A32, A33 - >::type; - -using orient3d = det3x3 - < - difference<_1, _10>, difference<_2, _11>, difference<_3, _12>, - difference<_4, _10>, difference<_5, _11>, difference<_6, _12>, - difference<_7, _10>, difference<_8, _11>, difference<_9, _12> + det, + boost::mp11::mp_transform_q + < + incircle_entry, + boost::mp11::mp_iota_c< (N + 1) * (N + 1) > + > >; -struct incircle_impl +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 adx = difference<_1, _7>; - using ady = difference<_2, _8>; - using bdx = difference<_3, _7>; - using bdy = difference<_4, _8>; - using cdx = difference<_5, _7>; - using cdy = difference<_6, _8>; - using alift = sum, product>; - using blift = sum, product>; - using clift = sum, product>; + 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 = det3x3 + using type = sum < - alift, adx, ady, - blift, bdx, bdy, - clift, cdx, cdy + difference, + sum + < + difference, + edet + > >; }; -using incircle = incircle_impl::type; +using insphere_no_translation = insphere_no_translation_impl::type; }} // namespace detail::generic_robust_predicates 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 index 9a9de3d1ff..5716017015 100644 --- 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 @@ -27,66 +27,66 @@ namespace detail { namespace generic_robust_predicates //expressions for static filters that compute error bounds for upper and lower //bounds on input values. -template +template struct interval_min_impl { using type = Expression; }; -template +template struct interval_max_impl { using type = Expression; }; -template -struct interval_min_impl, max_argn> +template +struct interval_min_impl, MaxArgn> { private: - using min_left = typename interval_min_impl::type; - using max_right = typename interval_max_impl::type; + 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, max_argn> +template +struct interval_min_impl, MaxArgn> { private: - using min_left = typename interval_min_impl::type; - using min_right = typename interval_min_impl::type; + 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, max_argn> +template +struct interval_max_impl, MaxArgn> { private: - using max_left = typename interval_max_impl::type; - using min_right = typename interval_min_impl::type; + 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, max_argn> +template +struct interval_max_impl, MaxArgn> { private: - using max_left = typename interval_max_impl::type; - using max_right = typename interval_max_impl::type; + 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, max_argn> +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; + 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 < @@ -103,12 +103,12 @@ struct interval_min_impl, max_argn> >; }; -template -struct interval_min_impl, max_argn> +template +struct interval_min_impl, MaxArgn> { private: - using min_child = typename interval_min_impl::type; - using max_child = typename interval_max_impl::type; + using min_child = typename interval_min_impl::type; + using max_child = typename interval_max_impl::type; public: using type = min < @@ -117,14 +117,14 @@ struct interval_min_impl, max_argn> >; }; -template -struct interval_max_impl, max_argn> +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; + 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 < @@ -141,12 +141,12 @@ struct interval_max_impl, max_argn> >; }; -template -struct interval_max_impl, max_argn> +template +struct interval_max_impl, MaxArgn> { private: - using min_child = typename interval_min_impl::type; - using max_child = typename interval_max_impl::type; + using min_child = typename interval_min_impl::type; + using max_child = typename interval_max_impl::type; public: using type = max < @@ -155,93 +155,103 @@ struct interval_max_impl, max_argn> >; }; -template -struct interval_min_impl, max_argn> +template +struct interval_min_impl, MaxArgn> { private: - using min_child = typename interval_min_impl::type; - using max_child = typename interval_max_impl::type; + 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, max_argn> +template +struct interval_max_impl, MaxArgn> { private: - using min_child = typename interval_min_impl::type; - using max_child = typename interval_max_impl::type; + 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, max_argn> +template +struct interval_min_impl, MaxArgn> { - using type = argument; + using type = argument; }; -template -struct interval_max_impl, max_argn> +template +struct interval_max_impl, MaxArgn> { using type = argument; }; -template +template struct interval_impl { using type = Expression; }; -template -struct interval_impl, max_argn> +template +struct interval_impl, MaxArgn> { private: - using min_child = typename interval_min_impl::type; - using max_child = typename interval_max_impl::type; + 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, max_argn> +template +struct interval_impl, MaxArgn> { private: - using left = typename interval_impl::type; - using right = typename interval_impl::type; + 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, max_argn> +template +struct interval_impl, MaxArgn> { private: - using left = typename interval_impl::type; - using right = typename interval_impl::type; + using left = typename interval_impl::type; + using right = typename interval_impl::type; public: using type = sum; }; -template -struct interval_impl, max_argn> +template +struct interval_impl, MaxArgn> { private: - using left = typename interval_impl::type; - using right = typename interval_impl::type; + 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_min = typename interval_min_impl::type; -template -using interval_max = typename interval_max_impl::type; +template +using interval_max = typename interval_max_impl::type; template using interval = - typename interval_impl::value>::type; + typename interval_impl>::type; }} // namespace detail::generic_robust_predicates diff --git a/include/boost/geometry/extensions/generic_robust_predicates/strategies/cartesian/detail/result_propagation.hpp b/include/boost/geometry/extensions/generic_robust_predicates/strategies/cartesian/detail/result_propagation.hpp deleted file mode 100644 index 4e25c0be96..0000000000 --- a/include/boost/geometry/extensions/generic_robust_predicates/strategies/cartesian/detail/result_propagation.hpp +++ /dev/null @@ -1,91 +0,0 @@ -// 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_RESULT_PROPAGATION_HPP -#define BOOST_GEOMETRY_EXTENSIONS_GENERIC_ROBUST_PREDICATES_STRATEGIES_CARTESIAN_DETAIL_RESULT_PROPAGATION_HPP - -#include -#include -#include - -#include - -namespace boost { namespace geometry -{ - -namespace detail { namespace generic_robust_predicates -{ - -template -< - std::size_t Last, - typename In = boost::mp11::mp_list<>, - std::size_t At = 1 -> -struct argument_list_impl -{ -private: - using next = boost::mp11::mp_push_back>; -public: - using type = typename argument_list_impl::type; -}; - -template -< - std::size_t Last, - typename In -> -struct argument_list_impl - < - Last, - In, - Last - > -{ - using type = boost::mp11::mp_push_back>; -}; - -template using argument_list = - typename argument_list_impl::type; - -template -< - typename Needle, - typename Haystacks -> -struct index_pair_impl -{ -private: - template - using contains = typename boost::mp11::mp_bind_back - < - boost::mp11::mp_contains, - Needle - >::template fn; - using outer = boost::mp11::mp_find_if; - using inner_list = boost::mp11::mp_at; - using inner = boost::mp11::mp_find; -public: - using type = std::pair; -}; - -template -< - typename Needle, - typename Haystacks -> -using index_pair = typename index_pair_impl::type; - -}} // namespace detail::generic_robust_predicates - -}} // namespace boost::geometry - -#endif // BOOST_GEOMETRY_EXTENSIONS_GENERIC_ROBUST_PREDICATES_STRATEGIES_CARTESIAN_DETAIL_RESULT_PROPAGATION_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 index a763be9e1c..e91d3cc477 100644 --- 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 @@ -19,8 +19,7 @@ #include #include -#include -#include +#include namespace boost { namespace geometry { @@ -36,7 +35,7 @@ namespace detail { namespace generic_robust_predicates //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 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 @@ -71,18 +70,30 @@ struct semi_static_filter >; using ct = CalculationType; public: - template - static inline int apply(const Reals&... args) + 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) { - using arg_list_input = argument_list; - using arg_list = boost::mp11::mp_list; - std::array input + std::array input {{ static_cast(args)... }}; std::array::value> results; - approximate_interim(results, input); - const ct error_bound = - get_approx(results, input); - const ct det = get_approx(results, input); + 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; 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 index 9fafd9f13b..e1e50ad501 100644 --- 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 @@ -12,13 +12,8 @@ #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 +#include #include #include @@ -34,72 +29,17 @@ 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 -struct stage_a_error_bound_impl -{ -private: - using ct = CalculationType; - using root = Expression; - using stack = typename boost::mp11::mp_unique>; - using evals = typename boost::mp11::mp_remove_if; - using interim_evals = typename boost::mp11::mp_remove - < - boost::mp11::mp_remove_if, - root - >; - using interim_errors = evals_error; - using final_children = add_children - < - boost::mp11::mp_second - < - boost::mp11::mp_map_find - < - interim_errors, - typename root::left - > - >, - boost::mp11::mp_second - < - boost::mp11::mp_map_find - < - interim_errors, - typename root::right - > - > - >; - using final_children_ltp = error_map_list_to_product; - using final_children_ltp_abs = abs_all; - using final_children_sum_fold = error_map_balanced_sum_up; - using final_coeff = coeff_round - < - div_by_1_m_eps - < - mult_by_1_p_eps - < - boost::mp11::mp_second - > - > - >; - static constexpr ct final_coeff_value = eval_eps_polynomial::value; - - struct final_coeff_constant : public static_constant_interface - { - static constexpr ct value = final_coeff_value; - static constexpr bool non_negative = true; - }; - - using error_expression_variable = boost::mp11::mp_front; -public: - using type = product; -}; - template < typename Expression, typename CalculationType > -using stage_a_error_bound = - typename stage_a_error_bound_impl::type; +using stage_a_error_bound_expression = + typename stage_a_error_bound_expression_impl + < + Expression, + CalculationType + >::type; template < @@ -107,11 +47,15 @@ template typename CalculationType > using stage_a_semi_static = semi_static_filter - < - Expression, - CalculationType, - stage_a_error_bound - >; + < + Expression, + CalculationType, + stage_a_error_bound_expression + < + Expression, + CalculationType + > + >; template < @@ -119,11 +63,18 @@ template typename CalculationType > using stage_a_static = static_filter - < - Expression, - CalculationType, - interval> - >; + < + Expression, + CalculationType, + interval + < + stage_a_error_bound_expression + < + Expression, + CalculationType + > + > + >; template < @@ -131,11 +82,11 @@ template typename CalculationType > using stage_a_almost_static = almost_static_filter - < - Expression, - CalculationType, - stage_a_static - >; + < + Expression, + CalculationType, + stage_a_static + >; }} // namespace detail::generic_robust_predicates 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 index 5acdbd9a5d..dce3f2ba58 100644 --- 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 @@ -12,14 +12,14 @@ #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 -#include +#include namespace boost { namespace geometry { @@ -47,46 +47,36 @@ class static_filter { private: using ct = CalculationType; - 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 - >; ct m_error_bound; - static constexpr std::size_t extrema_count = - max_argn::value; + static constexpr std::size_t extrema_count = max_argn; public: - static constexpr std::size_t arg_count = max_argn::value; + static constexpr bool stateful = true; + static constexpr bool updates = false; - inline static_filter() {} + inline static_filter() : m_error_bound(std::numeric_limits::max()) {} - ct error_bound() const { return m_error_bound; } + inline ct error_bound() const { return m_error_bound; } - template - inline static_filter(const Reals&... args) - : m_error_bound(approximate_value( + template + inline static_filter(CTs const&... args) + : m_error_bound(evaluate_expression( std::array {static_cast(args)...})) { - static_assert(sizeof...(Reals) == extrema_count, + static_assert(sizeof...(CTs) == extrema_count, "Number of constructor arguments is incompatible with error expression."); } - inline static_filter(const std::array& extrema) - : m_error_bound(approximate_value(extrema)) {} + inline static_filter(std::array const& extrema) + : m_error_bound(evaluate_expression(extrema)) {} - inline static_filter(const ct& b) : m_error_bound(b) {} + inline static_filter(ct const& b) : m_error_bound(b) {} - template - inline int apply(const Reals&... args) const + template + inline int apply(CTs const&... args) const { - using arg_list_input = argument_list; - using arg_list = boost::mp11::mp_list; - std::array input {static_cast(args)...}; - std::array::value> results; - approximate_interim(results, input); - const ct det = get_approx(results, input); + std::array input {static_cast(args)...}; + const ct det = evaluate_expression(input); if (det > m_error_bound) { return 1; @@ -105,8 +95,8 @@ class static_filter } } - template - inline int operator()(Reals const&... args) const + template + inline int operator()(CTs const&... args) const { return apply(args...); } diff --git a/include/boost/geometry/extensions/generic_robust_predicates/strategies/cartesian/detail/static_util.hpp b/include/boost/geometry/extensions/generic_robust_predicates/strategies/cartesian/detail/static_util.hpp deleted file mode 100644 index c23f548645..0000000000 --- a/include/boost/geometry/extensions/generic_robust_predicates/strategies/cartesian/detail/static_util.hpp +++ /dev/null @@ -1,130 +0,0 @@ -// 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_UTIL_HPP -#define BOOST_GEOMETRY_EXTENSIONS_GENERIC_ROBUST_PREDICATES_STRATEGIES_CARTESIAN_DETAIL_STATIC_UTIL_HPP - -#include -#include -#include - - -namespace boost { namespace geometry -{ - -namespace detail { namespace generic_robust_predicates -{ - -template -< - typename LOV, - typename Void = boost::mp11::mp_same -> -struct empty_or_void -{ - using type = boost::mp11::mp_empty; -}; - -template -struct empty_or_void -{ - using type = boost::mp11::mp_true; -}; - -template -struct log_2_floor_impl -{ - using type = typename boost::mp11::mp_plus - < - typename log_2_floor_impl>::type, - boost::mp11::mp_int<1> - >; -}; - -template <> -struct log_2_floor_impl> -{ - using type = boost::mp11::mp_int<0>; -}; - -template using log_2_floor = typename log_2_floor_impl::type; - -template -struct log_2_ceil_impl -{ -private: - using floor = log_2_floor; -public: - using type = boost::mp11::mp_int - < - (1 << floor::value) == N::value ? floor::value : floor::value + 1 - >; -}; - -template using log_2_ceil = typename log_2_ceil_impl::type; - -template using is_not_zero = boost::mp11::mp_bool; - -template -< - typename M, - typename K, - typename Contained = boost::mp11::mp_map_contains -> -struct mp_map_at_second_or_void -{ - using type = void; -}; - -template -struct mp_map_at_second_or_void -{ - using type = boost::mp11::mp_second>; -}; - -template using inc = boost::mp11::mp_int; - -template using indexed = boost::mp11::mp_transform - < - boost::mp11::mp_list, - boost::mp11::mp_iota>, - L - >; - -template using strip_index = - boost::mp11::mp_transform; - -template -< - typename Map, - typename Key, - typename Contains = boost::mp11::mp_map_contains -> -struct val_or_empty_list -{ - using type = boost::mp11::mp_second>; -}; - -template -< - typename Map, - typename Key -> -struct val_or_empty_list -{ - using type = boost::mp11::mp_list<>; -}; - -}} // namespace detail::generic_robust_predicates - -}} // namespace boost::geometry - -#endif // BOOST_GEOMETRY_EXTENSIONS_GENERIC_ROBUST_PREDICATES_STRATEGIES_CARTESIAN_DETAIL_STATIC_UTIL_HPP From da5b0ffbddd2e1eb5af316073f565eed19ff10a1 Mon Sep 17 00:00:00 2001 From: Tinko Bartels Date: Thu, 4 Mar 2021 21:05:23 +0100 Subject: [PATCH 6/6] Updates tests to reflect interface changes, adds test for staged predicate --- .../test/generic_robust_predicates/Jamfile | 5 ++- .../{approximate.cpp => expression_eval.cpp} | 25 ++++------- .../test/generic_robust_predicates/staged.cpp | 43 +++++++++++++++++++ 3 files changed, 55 insertions(+), 18 deletions(-) rename extensions/test/generic_robust_predicates/{approximate.cpp => expression_eval.cpp} (66%) create mode 100644 extensions/test/generic_robust_predicates/staged.cpp diff --git a/extensions/test/generic_robust_predicates/Jamfile b/extensions/test/generic_robust_predicates/Jamfile index 0595681dce..ef072ff343 100644 --- a/extensions/test/generic_robust_predicates/Jamfile +++ b/extensions/test/generic_robust_predicates/Jamfile @@ -6,7 +6,8 @@ test-suite boost-geometry-extensions-generic_robust_predicates : - [ run approximate.cpp ] - [ run side3d.cpp : : : off ] + [ run expression_eval.cpp ] + [ run side3d.cpp ] + [ run staged.cpp ] ; diff --git a/extensions/test/generic_robust_predicates/approximate.cpp b/extensions/test/generic_robust_predicates/expression_eval.cpp similarity index 66% rename from extensions/test/generic_robust_predicates/approximate.cpp rename to extensions/test/generic_robust_predicates/expression_eval.cpp index 4a12d3ecae..e6a2c5972f 100644 --- a/extensions/test/generic_robust_predicates/approximate.cpp +++ b/extensions/test/generic_robust_predicates/expression_eval.cpp @@ -14,8 +14,8 @@ #include -#include #include +#include template void test_all() @@ -30,15 +30,13 @@ void test_all() 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::approximate_value; - using bg::detail::generic_robust_predicates::get_approx; - using bg::detail::generic_robust_predicates::approximate_interim; - using bg::detail::generic_robust_predicates::argument_list; + using bg::detail::generic_robust_predicates::evaluate_expression; + using bg::detail::generic_robust_predicates::evaluate_expressions; using ct = CalculationType; - ct r1 = approximate_value, ct>( + ct r1 = evaluate_expression>( std::array{1.0, 2.0}); BOOST_CHECK_EQUAL(3.0, r1); - ct r2 = approximate_value, abs<_2>>, ct>( + ct r2 = evaluate_expression, abs<_2>>>( std::array{-10.0, 2.0}); BOOST_CHECK_EQUAL(10.0, r2); @@ -48,17 +46,12 @@ void test_all() difference<_3, _4> >; using evals = post_order; - using arg_list_input = argument_list<4>; - using arg_list = boost::mp11::mp_list; std::array::value> r; std::array input {5.0, 3.0, 2.0, 8.0}; - approximate_interim(r, input); - ct r3 = get_approx(r, input); - BOOST_CHECK_EQUAL(2.0, r3); - ct r4 = get_approx(r, input); - BOOST_CHECK_EQUAL(-6.0, r4); - ct r5 = get_approx(r, input); - BOOST_CHECK_EQUAL(-12.0, r5); + 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]); } 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; +}