Automatically generated, fast floating-point predicates

extensions/example/Jamfile

@@ -13,3 +13,4 @@ project boost-geometry-examples-extensions
     ;
 
 build-project gis ;
+build-project generic_robust_predicates ;

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 ;

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 <iostream>
+
+#include <boost/geometry/core/cs.hpp>
+#include <boost/geometry/geometries/point.hpp>
+
+#include <boost/geometry/extensions/triangulation/strategies/cartesian/side_robust.hpp>
+#include <boost/geometry/extensions/generic_robust_predicates/strategies/cartesian/detail/expressions.hpp>
+#include <boost/geometry/extensions/generic_robust_predicates/strategies/cartesian/detail/stage_a.hpp>
+
+
+namespace bg = boost::geometry;
+using point = bg::model::point<double, 2, bg::cs::cartesian>;
+
+template <typename CalculationType>
+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<double>::apply(p1, p2, p);
+        }
+    }
+};
+
+int main()
+{
+    point p1(0.0, 0.0);
+    point p2(1.0, 1.0);
+    point p (0.0, 1.0);
+    side_robust_with_static_filter<double> static_strategy(2.0, 2.0, 1.0, 1.0);
+    std::cout << "Side value: " << static_strategy.apply(p1, p2, p) << "\n"; 
+    return 0;
+}

extensions/test/Jamfile

@@ -28,3 +28,4 @@ build-project gis ;
 build-project iterators ;
 build-project nsphere ;
 build-project triangulation ;
+build-project generic_robust_predicates ;

extensions/test/generic_robust_predicates/Jamfile

@@ -0,0 +1,13 @@
+# Boost.Geometry (aka GGL, Generic Geometry Library)
+#
+# Use, modification and distribution is subject to the Boost Software License,
+# Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at
+# http://www.boost.org/LICENSE_1_0.txt)
+
+test-suite boost-geometry-extensions-generic_robust_predicates
+    :
+    [ run expression_eval.cpp ]
+    [ run side3d.cpp ]
+    [ run staged.cpp ]
+    ;
+

extensions/test/generic_robust_predicates/expression_eval.cpp

@@ -0,0 +1,62 @@
+// Boost.Geometry (aka GGL, Generic Geometry Library)
+// Unit Test
+
+// Copyright (c) 2020 Tinko Bartels, Berlin, Germany.
+
+// Contributed and/or modified by Tinko Bartels,
+//   as part of Google Summer of Code 2020 program.
+
+// Use, modification and distribution is subject to the Boost Software License,
+// Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at
+// http://www.boost.org/LICENSE_1_0.txt)
+
+#include <array>
+
+#include <geometry_test_common.hpp>
+
+#include <boost/geometry/extensions/generic_robust_predicates/strategies/cartesian/detail/expression_tree.hpp>
+#include <boost/geometry/extensions/generic_robust_predicates/strategies/cartesian/detail/expression_eval.hpp>
+
+template <typename CalculationType>
+void test_all()
+{
+    using bg::detail::generic_robust_predicates::_1;
+    using bg::detail::generic_robust_predicates::_2;
+    using bg::detail::generic_robust_predicates::_3;
+    using bg::detail::generic_robust_predicates::_4;
+    using bg::detail::generic_robust_predicates::sum;
+    using bg::detail::generic_robust_predicates::difference;
+    using bg::detail::generic_robust_predicates::product;
+    using bg::detail::generic_robust_predicates::max;
+    using bg::detail::generic_robust_predicates::abs;
+    using bg::detail::generic_robust_predicates::post_order;
+    using bg::detail::generic_robust_predicates::evaluate_expression;
+    using bg::detail::generic_robust_predicates::evaluate_expressions;
+    using ct = CalculationType;
+    ct r1 = evaluate_expression<sum<_1, _2>>(
+            std::array<ct, 2>{1.0, 2.0});
+    BOOST_CHECK_EQUAL(3.0, r1);
+    ct r2 = evaluate_expression<max<abs<_1>, abs<_2>>>(
+            std::array<ct, 2>{-10.0, 2.0});
+    BOOST_CHECK_EQUAL(10.0, r2);
+
+    using expression = product
+        <
+            difference<_1, _2>,
+            difference<_3, _4>
+        >;
+    using evals = post_order<expression>;
+    std::array<ct, boost::mp11::mp_size<evals>::value> r;
+    std::array<ct, 4> input {5.0, 3.0, 2.0, 8.0};
+    evaluate_expressions(input, r, evals{});
+    BOOST_CHECK_EQUAL(2.0, r[0]);
+    BOOST_CHECK_EQUAL(-6.0, r[1]);
+    BOOST_CHECK_EQUAL(-12.0, r[2]);
+}
+
+
+int test_main(int, char* [])
+{
+    test_all<double>();
+    return 0;
+}

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 <array>
+
+#include <geometry_test_common.hpp>
+
+#include <boost/geometry/extensions/generic_robust_predicates/strategies/cartesian/detail/expressions.hpp>
+#include <boost/geometry/extensions/generic_robust_predicates/strategies/cartesian/detail/stage_a.hpp>
+
+template <typename CalculationType>
+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<orient3d, ct>;
+    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<orient3d, ct> 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<orient3d, ct> 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<double>();
+    return 0;
+}

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 <geometry_test_common.hpp>
+
+#include <boost/geometry/extensions/generic_robust_predicates/strategies/cartesian/detail/expressions.hpp>
+#include <boost/geometry/extensions/generic_robust_predicates/strategies/cartesian/detail/stage_a.hpp>
+#include <boost/geometry/extensions/generic_robust_predicates/strategies/cartesian/detail/stage_d.hpp>
+#include <boost/geometry/extensions/generic_robust_predicates/strategies/cartesian/detail/stage_b.hpp>
+#include <boost/geometry/extensions/generic_robust_predicates/strategies/cartesian/detail/staged_predicate.hpp>
+
+using namespace boost::geometry::detail::generic_robust_predicates;
+
+template <typename CalculationType>
+void test_all()
+{
+    using ct = CalculationType;
+    using expression = orient2d;
+    using filter1 = stage_a_static<expression, ct>;
+    using filter2 = stage_a_almost_static<expression, ct>;
+    using filter3 = stage_a_semi_static<expression, ct>;
+    using filter4 = stage_b<expression, ct>;
+    using filter5 = stage_d<expression, ct>;
+    using staged = staged_predicate<ct, filter1, filter2, filter3, filter4, filter5>;
+    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<double>();
+    return 0;
+}

include/boost/geometry/extensions/generic_robust_predicates/strategies/cartesian/detail/almost_static_filter.hpp

@@ -0,0 +1,129 @@
+// Boost.Geometry (aka GGL, Generic Geometry Library)
+
+// Copyright (c) 2020 Tinko Bartels, Berlin, Germany.
+
+// Contributed and/or modified by Tinko Bartels,
+//   as part of Google Summer of Code 2020 program.
+
+// Use, modification and distribution is subject to the Boost Software License,
+// Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at
+// http://www.boost.org/LICENSE_1_0.txt)
+
+#ifndef BOOST_GEOMETRY_EXTENSIONS_GENERIC_ROBUST_PREDICATES_STRATEGIES_CARTESIAN_DETAIL_ALMOST_STATIC_FILTER_HPP
+#define BOOST_GEOMETRY_EXTENSIONS_GENERIC_ROBUST_PREDICATES_STRATEGIES_CARTESIAN_DETAIL_ALMOST_STATIC_FILTER_HPP
+
+#include <array>
+#include <algorithm>
+#include <limits>
+
+#include <boost/geometry/extensions/generic_robust_predicates/strategies/cartesian/detail/expression_tree.hpp>
+
+namespace boost { namespace geometry
+{
+
+namespace detail { namespace generic_robust_predicates
+{
+
+//The almost static filter holds an instance of a static filter and applies it
+//when it is called. Its static filter can be updated and applied like this:
+//
+//almost_static_filter<...> f;
+//
+//f.update(max1, max2, ..., min1, min2, ...);
+//
+//f.apply(arg1, arg2, ...);
+//
+//Unlike with the static filter, global bounds do not need to be known at
+//construction time and with incremental algorithms where inputs with higher
+//magnitude are added later, the earlier applications of the filter may benefit
+//from more optimistic error bounds.
+
+template
+<
+    typename Expression,
+    typename CalculationType,
+    typename StaticFilter
+>
+class almost_static_filter
+{
+private:
+    using ct = CalculationType;
+    using extrema_array = std::array<ct, 2 * max_argn<Expression>>;
+    extrema_array m_extrema;
+    StaticFilter m_filter;
+public:
+    static constexpr bool stateful = true;
+    static constexpr bool updates = true;
+
+    StaticFilter const& filter() const { return m_filter; }
+    inline almost_static_filter()
+    {
+        std::fill(m_extrema.begin(),
+                  m_extrema.begin() + m_extrema.size() / 2,
+                  -std::numeric_limits<ct>::infinity());
+        std::fill(m_extrema.begin() + m_extrema.size() / 2,
+                  m_extrema.end(),
+                  std::numeric_limits<ct>::infinity());
+    }
+    template <typename ...CTs>
+    int apply(CTs const&... args) const
+    {
+        return m_filter.apply(args...);
+    }
+
+    template <typename ...CTs>
+    inline void update_extrema(CTs const&... args)
+    {
+        std::array<ct, sizeof...(CTs)> input {{ static_cast<ct>(args)... }};
+        for(unsigned int i = 0; i < m_extrema.size() / 2; ++i)
+        {
+            m_extrema[i] = std::max(m_extrema[i], input[i]);
+        }
+        for(unsigned int i = m_extrema.size() / 2; i < m_extrema.size(); ++i)
+        {
+            m_extrema[i] = std::min(m_extrema[i], input[i]);
+        }
+    }
+
+    template <typename ...CTs>
+    inline bool update_extrema_check(CTs const&... args)
+    {
+        bool changed = false;
+        std::array<ct, sizeof...(CTs)> input {{ static_cast<ct>(args)... }};
+        for(unsigned int i = 0; i < m_extrema.size() / 2; ++i)
+        {
+            if (input[i] > m_extrema[i])
+            {
+                changed = true;
+                m_extrema[i] = input[i];
+            }
+        }
+        for(unsigned int i = m_extrema.size() / 2; i < m_extrema.size(); ++i)
+        {
+            if (input[i] < m_extrema[i])
+            {
+                changed = true;
+                m_extrema[i] = input[i];
+            }
+        }
+        return changed;
+    }
+
+    template <typename ...CTs>
+    inline void update(CTs const&... args)
+    {
+        update_extrema(args...);
+        update_filter();
+    }
+
+    inline void update_filter()
+    {
+        m_filter = StaticFilter(m_extrema);
+    }
+};
+
+}} // namespace detail::generic_robust_predicates
+
+}} // namespace boost::geometry
+
+#endif // BOOST_GEOMETRY_EXTENSIONS_GENERIC_ROBUST_PREDICATES_STRATEGIES_CARTESIAN_DETAIL_ALMOST_STATIC_FILTER_HPP

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 <cstddef>
+#include <cmath>
+#include <array>
+
+#include <boost/mp11/integral.hpp>
+#include <boost/mp11/list.hpp>
+#include <boost/mp11/algorithm.hpp>
+
+#include <boost/geometry/extensions/generic_robust_predicates/strategies/cartesian/detail/expression_tree.hpp>
+
+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 <operator_types Op>
+struct evaluate_expression_unary_impl {};
+
+template <>
+struct evaluate_expression_unary_impl<operator_types::abs>
+{
+    template <typename CT>
+    static constexpr void apply(CT const& c, CT& out)
+    {
+        out = std::abs(c);
+    }
+};
+
+template <operator_types Op>
+struct evaluate_expression_binary_impl {};
+
+template <>
+struct evaluate_expression_binary_impl<operator_types::sum>
+{
+    template <typename CT>
+    static constexpr void apply(CT const& l, CT const& r, CT& out)
+    {
+        out = l + r;
+    }
+};
+
+template <>
+struct evaluate_expression_binary_impl<operator_types::product>
+{
+    template <typename CT>
+    static constexpr void apply(CT const& l, CT const& r, CT& out)
+    {
+        out = l * r;
+    }
+};
+
+template <>
+struct evaluate_expression_binary_impl<operator_types::difference>
+{
+    template <typename CT>
+    static constexpr void apply(CT const& l, CT const& r, CT& out)
+    {
+        out = l - r;
+    }
+};
+
+template <>
+struct evaluate_expression_binary_impl<operator_types::min>
+{
+    template <typename CT>
+    static constexpr void apply(CT const& l, CT const& r, CT& out)
+    {
+        out = std::min(l, r);
+    }
+};
+
+template <>
+struct evaluate_expression_binary_impl<operator_types::max>
+{
+    template <typename CT>
+    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 <std::size_t I, typename Expression>
+struct get_arg_or_out_impl<I, Expression, false>
+{
+    template <typename In, typename Out>
+    static constexpr auto apply(In const&, Out const& out)
+    {
+        return out[I];
+    }
+};
+
+template <typename Expression, std::size_t Argn = Expression::argn>
+struct get_arg_or_const
+{
+    template <typename In>
+    static constexpr auto apply(In const& in)
+    {
+        return in[Argn - 1];
+    }
+};
+
+template <typename Expression>
+struct get_arg_or_const<Expression, 0>
+{
+    template <typename In>
+    static constexpr auto apply(In const&)
+    {
+        return Expression::value;
+    }
+};
+
+template <std::size_t I, typename Expression>
+struct get_arg_or_out_impl<I, Expression, true>
+{
+    template <typename In, typename Out>
+    static constexpr auto apply(In const& in, Out const&)
+    {
+        return get_arg_or_const<Expression>::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<Expressions, Expression>::value;
+    static constexpr std::size_t i_left =
+        boost::mp11::mp_find<Expressions, left>::value;
+    static constexpr std::size_t i_right =
+        boost::mp11::mp_find<Expressions, right>::value;
+    static constexpr void apply(
+            InputArr const& input,
+            OutputArr& output
+        )
+    {
+        auto& out = output[i_out];
+        auto const& l =
+            get_arg_or_out_impl<i_left, left>::apply(input, output);
+        auto const& r =
+            get_arg_or_out_impl<i_right, right>::apply(input, output);
+        evaluate_expression_binary_impl<Expression::operator_type>
+            ::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<Expressions, Expression>::value;
+    static constexpr std::size_t i_child =
+        boost::mp11::mp_find<Expressions, child>::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<Expression::operator_type>
+            ::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<Expressions...>)
+{
+    using exp_list = boost::mp11::mp_list<Expressions...>;
+    //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<post_order<Expression>>;
+    using evals = typename boost::mp11::mp_remove_if<stack, is_leaf>;
+    std::array
+        <
+            typename InputArr::value_type,
+            boost::mp11::mp_size<evals>::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

include/boost/geometry/extensions/generic_robust_predicates/strategies/cartesian/detail/expression_tree.hpp

@@ -0,0 +1,268 @@
+// Boost.Geometry (aka GGL, Generic Geometry Library)
+
+// Copyright (c) 2020 Tinko Bartels, Berlin, Germany.
+
+// Contributed and/or modified by Tinko Bartels,
+//   as part of Google Summer of Code 2020 program.
+
+// Use, modification and distribution is subject to the Boost Software License,
+// Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at
+// http://www.boost.org/LICENSE_1_0.txt)
+
+#ifndef BOOST_GEOMETRY_EXTENSIONS_GENERIC_ROBUST_PREDICATES_STRATEGIES_CARTESIAN_DETAIL_EXPRESSION_TREE_HPP
+#define BOOST_GEOMETRY_EXTENSIONS_GENERIC_ROBUST_PREDICATES_STRATEGIES_CARTESIAN_DETAIL_EXPRESSION_TREE_HPP
+
+#include <cstddef>
+#include <type_traits>
+#include <algorithm>
+
+#include <boost/mp11/algorithm.hpp>
+#include <boost/mp11/integral.hpp>
+#include <boost/mp11/list.hpp>
+
+namespace boost { namespace geometry
+{
+
+namespace detail { namespace generic_robust_predicates
+{
+
+//The templates in this file are designed to represent arithmetic expressions
+//In the type system so that they can be processed at compile-time for the
+//automatic creation of floating-point filters.
+//
+//E.g. the arithmetic expression a * b + c can be represented by the type
+//
+//sum< product < argument<1>, argument<2> >, argument<3> >
+//
+//The arguments represent the input variables and are the leaves of the
+//expression tree (other possible leaves are constants). Because they represent
+//placeholders their indexing is one-based rather than zero-based, similar to
+//the placeholders for std::bind.
+
+enum class operator_types {
+    sum, difference, product, abs, no_op, max, min
+};
+
+enum class operator_arities { nullary, unary, binary };
+
+constexpr int sign_uncertain = -2;
+
+struct sum_error_type {};
+struct product_error_type {};
+struct no_error_type {};
+
+template <typename... Children>
+struct internal_node
+{
+    static constexpr bool is_leaf = false;
+    static constexpr bool non_negative = false;
+    using all_children = boost::mp11::mp_list<Children...>; //for convenience
+};
+
+template <typename Left, typename Right>
+struct internal_binary_node : internal_node<Left, Right>
+{
+    using left  = Left;
+    using right = Right;
+    static constexpr operator_arities operator_arity =
+        operator_arities::binary;
+};
+
+template <typename Child>
+struct internal_unary_node : internal_node<Child>
+{
+    using child = Child;
+    static constexpr operator_arities operator_arity = operator_arities::unary;
+};
+
+template <typename Left, typename Right>
+struct sum : public internal_binary_node<Left, Right>
+{
+    static constexpr bool sign_exact =
+           (Left::is_leaf && Right::is_leaf)
+        || (   Left::sign_exact && Right::sign_exact
+            && Left::non_negative && Right::non_negative);
+    static constexpr bool non_negative =    Left::non_negative
+                                         && Right::non_negative;
+    static constexpr operator_types operator_type = operator_types::sum;
+    using error_type = sum_error_type;
+};
+
+template <typename Left, typename Right>
+struct difference : public internal_binary_node<Left, Right>
+{
+    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 <typename Left, typename Right>
+struct product : public internal_binary_node<Left, Right>
+{
+    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<Left, Right>::value;
+    static constexpr operator_types operator_type = operator_types::product;
+    using error_type = product_error_type;
+};
+
+template <typename Left, typename Right>
+struct max : public internal_binary_node<Left, Right>
+{
+    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 <typename Left, typename Right>
+struct min : public internal_binary_node<Left, Right>
+{
+    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 <typename Child>
+struct abs : public internal_unary_node<Child>
+{
+    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 <std::size_t Argn>
+struct argument : public leaf
+{
+    static constexpr std::size_t argn = Argn;
+};
+
+template <typename NumberType>
+struct static_constant_interface : public leaf
+{
+    using value_type = NumberType;
+    static constexpr NumberType value = 0; //override
+    static constexpr std::size_t argn = 0;
+};
+
+template <typename Node>
+using is_leaf = boost::mp11::mp_bool<Node::is_leaf>;
+
+//post_order_impl and post_order are templates for compile-time traversion of
+//expression trees. The post_order traversal was chosen because it guarantees
+//that children (subexpressions) are evaluated before there parents as is
+//necessary for the evaluation of the arithmetic expressions that these
+//expression trees represent.
+
+template
+<
+    typename In,
+    typename Out,
+    template <typename> class Anchor = is_leaf,
+    bool IsBinary = In::operator_arity == operator_arities::binary,
+    bool AtAnchor = Anchor<In>::value
+>
+struct post_order_impl;
+
+template
+<
+    typename In,
+    typename Out,
+    template <typename> class Anchor,
+    bool IsBinary
+>
+struct post_order_impl<In, Out, Anchor, IsBinary, true>
+{
+    using type = Out;
+};
+
+template <typename In, typename Out, template <typename> class Anchor>
+struct post_order_impl<In, Out, Anchor, true, false>
+{
+    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<Out, leftl, rightl>;
+    using type   =
+        boost::mp11::mp_unique<boost::mp11::mp_push_back<merged, In>>;
+};
+
+template <typename In, typename Out, template <typename> class Anchor>
+struct post_order_impl<In, Out, Anchor, false, false>
+{
+    using childl = typename post_order_impl
+            <
+                typename In::child,
+                boost::mp11::mp_list<>,
+                Anchor
+            >::type;
+    using merged = boost::mp11::mp_append<Out, childl>;
+    using type   =
+        boost::mp11::mp_unique<boost::mp11::mp_push_back<merged, In>>;
+};
+
+template <typename In, template <typename> class Anchor = is_leaf>
+using post_order =
+    typename post_order_impl<In, boost::mp11::mp_list<>, Anchor>::type;
+
+template <typename Expression, operator_arities = Expression::operator_arity>
+constexpr std::size_t max_argn = 0;
+
+template <typename Expression>
+constexpr std::size_t max_argn<Expression, operator_arities::binary> =
+    std::max(max_argn<typename Expression::left>,
+             max_argn<typename Expression::right>);
+
+template <typename Expression>
+constexpr std::size_t max_argn<Expression, operator_arities::unary> =
+    max_argn<typename Expression::child>;
+
+template <typename Expression>
+constexpr std::size_t max_argn<Expression, operator_arities::nullary> =
+    Expression::argn;
+
+using  _1 = argument<1>;
+using  _2 = argument<2>;
+using  _3 = argument<3>;
+using  _4 = argument<4>;
+using  _5 = argument<5>;
+using  _6 = argument<6>;
+using  _7 = argument<7>;
+using  _8 = argument<8>;
+using  _9 = argument<9>;
+using _10 = argument<10>;
+using _11 = argument<11>;
+using _12 = argument<12>;
+using _13 = argument<13>;
+using _14 = argument<14>;
+using _15 = argument<15>;
+
+}} // namespace detail::generic_robust_predicates
+
+}} // namespace boost::geometry
+
+#endif // BOOST_GEOMETRY_EXTENSIONS_GENERIC_ROBUST_PREDICATES_STRATEGIES_CARTESIAN_DETAIL_EXPRESSION_TREE_HPP

include/boost/geometry/extensions/generic_robust_predicates/strategies/cartesian/detail/interval_error_bound.hpp

@@ -0,0 +1,260 @@
+// Boost.Geometry (aka GGL, Generic Geometry Library)
+
+// Copyright (c) 2020 Tinko Bartels, Berlin, Germany.
+
+// Contributed and/or modified by Tinko Bartels,
+//   as part of Google Summer of Code 2020 program.
+
+// Use, modification and distribution is subject to the Boost Software License,
+// Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at
+// http://www.boost.org/LICENSE_1_0.txt)
+
+#ifndef BOOST_GEOMETRY_EXTENSIONS_GENERIC_ROBUST_PREDICATES_STRATEGIES_CARTESIAN_DETAIL_INTERVAL_ERROR_BOUND_HPP
+#define BOOST_GEOMETRY_EXTENSIONS_GENERIC_ROBUST_PREDICATES_STRATEGIES_CARTESIAN_DETAIL_INTERVAL_ERROR_BOUND_HPP
+
+#include <boost/geometry/extensions/generic_robust_predicates/strategies/cartesian/detail/expression_tree.hpp>
+
+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 <typename Expression, std::size_t MaxArgn>
+struct interval_min_impl
+{
+    using type = Expression;
+};
+
+template <typename Expression, std::size_t MaxArgn>
+struct interval_max_impl
+{
+    using type = Expression;
+};
+
+template <typename Left, typename Right, std::size_t MaxArgn>
+struct interval_min_impl<difference<Left, Right>, MaxArgn>
+{
+private:
+    using min_left = typename interval_min_impl<Left, MaxArgn>::type;
+    using max_right = typename interval_max_impl<Right, MaxArgn>::type;
+public:
+    using type = difference<min_left, max_right>;
+};
+
+template <typename Left, typename Right, std::size_t MaxArgn>
+struct interval_min_impl<sum<Left, Right>, MaxArgn>
+{
+private:
+    using min_left = typename interval_min_impl<Left, MaxArgn>::type;
+    using min_right = typename interval_min_impl<Right, MaxArgn>::type;
+public:
+    using type = sum<min_left, min_right>;
+};
+
+template <typename Left, typename Right, std::size_t MaxArgn>
+struct interval_max_impl<difference<Left, Right>, MaxArgn>
+{
+private:
+    using max_left = typename interval_max_impl<Left, MaxArgn>::type;
+    using min_right = typename interval_min_impl<Right, MaxArgn>::type;
+public:
+    using type = difference<max_left, min_right>;
+};
+
+template <typename Left, typename Right, std::size_t MaxArgn>
+struct interval_max_impl<sum<Left, Right>, MaxArgn>
+{
+private:
+    using max_left = typename interval_max_impl<Left, MaxArgn>::type;
+    using max_right = typename interval_max_impl<Right, MaxArgn>::type;
+public:
+    using type = sum<max_left, max_right>;
+};
+
+template <typename Left, typename Right, std::size_t MaxArgn>
+struct interval_min_impl<product<Left, Right>, MaxArgn>
+{
+private:
+    using min_left = typename interval_min_impl<Left, MaxArgn>::type;
+    using max_left = typename interval_max_impl<Left, MaxArgn>::type;
+    using min_right = typename interval_min_impl<Right, MaxArgn>::type;
+    using max_right = typename interval_max_impl<Right, MaxArgn>::type;
+public:
+    using type = min
+        <
+            min
+                <
+                    product<min_left, min_right>,
+                    product<min_left, max_right>
+                >,
+            min
+                <
+                    product<max_left, min_right>,
+                    product<max_left, max_right>
+                >
+        >;
+};
+
+template <typename Child, std::size_t MaxArgn>
+struct interval_min_impl<product<Child, Child>, MaxArgn>
+{
+private:
+    using min_child = typename interval_min_impl<Child, MaxArgn>::type;
+    using max_child = typename interval_max_impl<Child, MaxArgn>::type;
+public:
+    using type = min
+        <
+            product<min_child, min_child>,
+            product<max_child, max_child>
+        >;
+};
+
+template <typename Left, typename Right, std::size_t MaxArgn>
+struct interval_max_impl<product<Left, Right>, MaxArgn>
+{
+private:
+    using min_left = typename interval_min_impl<Left, MaxArgn>::type;
+    using max_left = typename interval_max_impl<Left, MaxArgn>::type;
+    using min_right = typename interval_min_impl<Right, MaxArgn>::type;
+    using max_right = typename interval_max_impl<Right, MaxArgn>::type;
+public:
+    using type = max
+        <
+            max
+                <
+                    product<min_left, min_right>,
+                    product<min_left, max_right>
+                >,
+            max
+                <
+                    product<max_left, min_right>,
+                    product<max_left, max_right>
+                >
+        >;
+};
+
+template <typename Child, std::size_t MaxArgn>
+struct interval_max_impl<product<Child, Child>, MaxArgn>
+{
+private:
+    using min_child = typename interval_min_impl<Child, MaxArgn>::type;
+    using max_child = typename interval_max_impl<Child, MaxArgn>::type;
+public:
+    using type = max
+        <
+            product<min_child, min_child>,
+            product<max_child, max_child>
+        >;
+};
+
+template <typename Child, std::size_t MaxArgn>
+struct interval_min_impl<abs<Child>, MaxArgn>
+{
+private:
+    using min_child = typename interval_min_impl<Child, MaxArgn>::type;
+    using max_child = typename interval_max_impl<Child, MaxArgn>::type;
+public:
+    using type = min<abs<min_child>, abs<max_child>>;
+};
+
+template <typename Child, std::size_t MaxArgn>
+struct interval_max_impl<abs<Child>, MaxArgn>
+{
+private:
+    using min_child = typename interval_min_impl<Child, MaxArgn>::type;
+    using max_child = typename interval_max_impl<Child, MaxArgn>::type;
+public:
+    using type = max<abs<min_child>, abs<max_child>>;
+};
+
+template <std::size_t argn, std::size_t MaxArgn>
+struct interval_min_impl<argument<argn>, MaxArgn>
+{
+    using type = argument<argn + MaxArgn>;
+};
+
+template <std::size_t argn, std::size_t MaxArgn>
+struct interval_max_impl<argument<argn>, MaxArgn>
+{
+    using type = argument<argn>;
+};
+
+template <typename Expression, std::size_t MaxArgn>
+struct interval_impl
+{
+    using type = Expression;
+};
+
+template <typename Child, std::size_t MaxArgn>
+struct interval_impl<abs<Child>, MaxArgn>
+{
+private:
+    using min_child = typename interval_min_impl<Child, MaxArgn>::type;
+    using max_child = typename interval_max_impl<Child, MaxArgn>::type;
+public:
+    using type = max<abs<min_child>, abs<max_child>>;
+};
+
+template <typename Left, typename Right, std::size_t MaxArgn>
+struct interval_impl<max<Left, Right>, MaxArgn>
+{
+private:
+    using left = typename interval_impl<Left, MaxArgn>::type;
+    using right = typename interval_impl<Right, MaxArgn>::type;
+public:
+    using type = max<left, right>;
+};
+
+template <typename Left, typename Right, std::size_t MaxArgn>
+struct interval_impl<difference<Left, Right>, MaxArgn>
+{
+private:
+    using left = typename interval_impl<Left, MaxArgn>::type;
+    using right = typename interval_impl<Right, MaxArgn>::type;
+public:
+    using type = difference<left, right>;
+};
+
+template <typename Left, typename Right, std::size_t MaxArgn>
+struct interval_impl<sum<Left, Right>, MaxArgn>
+{
+private:
+    using left = typename interval_impl<Left, MaxArgn>::type;
+    using right = typename interval_impl<Right, MaxArgn>::type;
+public:
+    using type = sum<left, right>;
+};
+
+template <typename Left, typename Right, std::size_t MaxArgn>
+struct interval_impl<product<Left, Right>, MaxArgn>
+{
+private:
+    using left = typename interval_impl<Left, MaxArgn>::type;
+    using right = typename interval_impl<Right, MaxArgn>::type;
+public:
+    using type = product<left, right>;
+};
+
+template <typename Expression, std::size_t MaxArgn>
+using interval_min = typename interval_min_impl<Expression, MaxArgn>::type;
+
+template <typename Expression, std::size_t MaxArgn>
+using interval_max = typename interval_max_impl<Expression, MaxArgn>::type;
+
+template <typename Expression>
+using interval =
+    typename interval_impl<Expression, max_argn<Expression>>::type;
+
+}} // namespace detail::generic_robust_predicates
+
+}} // namespace boost::geometry
+
+#endif // BOOST_GEOMETRY_EXTENSIONS_GENERIC_ROBUST_PREDICATES_STRATEGIES_CARTESIAN_DETAIL_INTERVAL_ERROR_BOUND_HPP

include/boost/geometry/extensions/generic_robust_predicates/strategies/cartesian/detail/semi_static_filter.hpp

@@ -0,0 +1,120 @@
+// Boost.Geometry (aka GGL, Generic Geometry Library)
+
+// Copyright (c) 2020 Tinko Bartels, Berlin, Germany.
+
+// Contributed and/or modified by Tinko Bartels,
+//   as part of Google Summer of Code 2020 program.
+
+// Use, modification and distribution is subject to the Boost Software License,
+// Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at
+// http://www.boost.org/LICENSE_1_0.txt)
+
+#ifndef BOOST_GEOMETRY_EXTENSIONS_GENERIC_ROBUST_PREDICATES_STRATEGIES_CARTESIAN_DETAIL_SEMI_STATIC_FILTER_HPP
+#define BOOST_GEOMETRY_EXTENSIONS_GENERIC_ROBUST_PREDICATES_STRATEGIES_CARTESIAN_DETAIL_SEMI_STATIC_FILTER_HPP
+
+#include <array>
+
+#include <boost/mp11/list.hpp>
+#include <boost/mp11/algorithm.hpp>
+#include <boost/mp11/set.hpp>
+
+#include <boost/geometry/extensions/generic_robust_predicates/strategies/cartesian/detail/expression_tree.hpp>
+#include <boost/geometry/extensions/generic_robust_predicates/strategies/cartesian/detail/expression_eval.hpp>
+
+namespace boost { namespace geometry
+{
+
+namespace detail { namespace generic_robust_predicates
+{
+
+//The semi static filter evaluates an Expression and an ErrorExpression. On
+//application it evaluates both expressions and verifies the sign of the
+//Expression by the following rules:
+//If the absolute value of the number resulting from the evaluation of the
+//expression is larger than or equal to the result for the error expression,
+//then the sign can be returned. If not, then a constant is returned to
+//represent an uncertain sign.
+//
+//The filters that are build from this template are meant to be semi-static in
+//the sense that the error expression including constants that depend only the
+//epsilon of the calculation type are known statically at compile-time but the
+//final value of the error bound depends on the specific inputs for each call
+//to the filter and parts of it need to be computed at each call.
+//
+//This filter is stateless.
+
+template
+<
+    typename Expression,
+    typename CalculationType,
+    typename ErrorExpression
+>
+struct semi_static_filter
+{
+private:
+    using stack = typename boost::mp11::mp_unique<post_order<Expression>>;
+    using evals = typename boost::mp11::mp_remove_if<stack, is_leaf>;
+    using error_eval_stack = boost::mp11::mp_unique
+        <
+            post_order<ErrorExpression>
+        >;
+    using error_eval_stack_remainder = boost::mp11::mp_set_difference
+        <
+            error_eval_stack,
+            evals
+        >;
+    using all_evals = boost::mp11::mp_append
+        <
+            evals,
+            error_eval_stack_remainder
+        >;
+    using ct = CalculationType;
+public:
+    static constexpr bool stateful = false;
+    static constexpr bool updates = false;
+
+    template <typename ...CTs>
+    static inline ct error_bound(CTs const&... args)
+    {
+        std::array<ct, sizeof...(CTs)> input
+            {{ static_cast<ct>(args)... }};
+        return evaluate_expression<ErrorExpression>(input);
+    }
+
+    template <typename ...CTs>
+    static inline int apply(CTs const&... args)
+    {
+        std::array<ct, sizeof...(CTs)> input
+            {{ static_cast<ct>(args)... }};
+        std::array<ct, boost::mp11::mp_size<all_evals>::value> results;
+        evaluate_expressions(input, results, all_evals{});
+        constexpr std::size_t i_eb =
+            boost::mp11::mp_find<all_evals, ErrorExpression>::value;
+        const ct error_bound = results[i_eb];
+        constexpr std::size_t i_e =
+            boost::mp11::mp_find<all_evals, Expression>::value;
+        const ct det = results[i_e];
+        if (det > error_bound)
+        {
+            return 1;
+        }
+        else if (det < -error_bound)
+        {
+            return -1;
+        }
+        else if (error_bound == 0 && det == 0)
+        {
+            return 0;
+        }
+        else
+        {
+            return sign_uncertain;
+        }
+    }
+};
+
+}} // namespace detail::generic_robust_predicates
+
+}} // namespace boost::geometry
+
+#endif // BOOST_GEOMETRY_EXTENSIONS_GENERIC_ROBUST_PREDICATES_STRATEGIES_CARTESIAN_DETAIL_SEMI_STATIC_FILTER_HPP

include/boost/geometry/extensions/generic_robust_predicates/strategies/cartesian/detail/stage_a.hpp

@@ -0,0 +1,95 @@
+// Boost.Geometry (aka GGL, Generic Geometry Library)
+
+// Copyright (c) 2020 Tinko Bartels, Berlin, Germany.
+
+// Contributed and/or modified by Tinko Bartels,
+//   as part of Google Summer of Code 2020 program.
+
+// Use, modification and distribution is subject to the Boost Software License,
+// Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at
+// http://www.boost.org/LICENSE_1_0.txt)
+
+#ifndef BOOST_GEOMETRY_EXTENSIONS_GENERIC_ROBUST_PREDICATES_STRATEGIES_CARTESIAN_DETAIL_STAGE_A_HPP
+#define BOOST_GEOMETRY_EXTENSIONS_GENERIC_ROBUST_PREDICATES_STRATEGIES_CARTESIAN_DETAIL_STAGE_A_HPP
+
+#include <boost/geometry/extensions/generic_robust_predicates/strategies/cartesian/detail/expression_tree.hpp>
+#include <boost/geometry/extensions/generic_robust_predicates/strategies/cartesian/detail/stage_a_error_bound.hpp>
+#include <boost/geometry/extensions/generic_robust_predicates/strategies/cartesian/detail/interval_error_bound.hpp>
+
+#include <boost/geometry/extensions/generic_robust_predicates/strategies/cartesian/detail/semi_static_filter.hpp>
+#include <boost/geometry/extensions/generic_robust_predicates/strategies/cartesian/detail/almost_static_filter.hpp>
+#include <boost/geometry/extensions/generic_robust_predicates/strategies/cartesian/detail/static_filter.hpp>
+
+namespace boost { namespace geometry
+{
+
+namespace detail { namespace generic_robust_predicates
+{
+
+//The following templates apply the ideas implemented in error_bound.hpp and
+//interval_error_bound.hpp to generate error expressions from expressions.
+
+template
+<
+    typename Expression,
+    typename CalculationType
+>
+using stage_a_error_bound_expression =
+    typename stage_a_error_bound_expression_impl
+        <
+            Expression,
+            CalculationType
+        >::type;
+
+template
+<
+    typename Expression,
+    typename CalculationType
+>
+using stage_a_semi_static = semi_static_filter
+    <
+        Expression,
+        CalculationType,
+        stage_a_error_bound_expression
+            <
+                Expression,
+                CalculationType
+            >
+    >;
+
+template
+<
+    typename Expression,
+    typename CalculationType
+>
+using stage_a_static = static_filter
+    <
+        Expression,
+        CalculationType,
+        interval
+            <
+                stage_a_error_bound_expression
+                    <
+                        Expression,
+                        CalculationType
+                    >
+            >
+    >;
+
+template
+<
+    typename Expression,
+    typename CalculationType
+>
+using stage_a_almost_static = almost_static_filter
+    <
+        Expression,
+        CalculationType,
+        stage_a_static<Expression, CalculationType>
+    >;
+
+}} // namespace detail::generic_robust_predicates
+
+}} // namespace boost::geometry
+
+#endif // BOOST_GEOMETRY_EXTENSIONS_GENERIC_ROBUST_PREDICATES_STRATEGIES_CARTESIAN_DETAIL_STAGE_A_HPP

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 <array>
+#include <limits>
+
+#include <boost/mp11/utility.hpp>
+
+#include <boost/geometry/extensions/generic_robust_predicates/strategies/cartesian/detail/expression_tree.hpp>
+
+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 <typename Expression>
+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<Expression>
+>
+struct stage_a_error_bound {};
+
+template <typename Expression>
+struct stage_a_error_bound<Expression, stage_a_error_propagation_cases::exact>
+{
+    using magnitude = abs<Expression>;
+    static constexpr std::array<int, 3> a {0, 0, 0};
+};
+
+template <typename Expression>
+struct stage_a_error_bound
+    <
+        Expression,
+        stage_a_error_propagation_cases::op_on_exacts
+    >
+{
+    using magnitude = abs<Expression>;
+    static constexpr std::array<int, 3> a {1, 0, 0};
+};
+
+constexpr std::array<int, 3> coeff_max(const std::array<int, 3> a,
+                                       const std::array<int, 3> 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<int, 3> coeff_product(const std::array<int, 3> a,
+                                           const std::array<int, 3> b)
+{
+    return std::array<int, 3> {
+            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<int, 3>
+coeff_mult_by_1_plus_eps(const std::array<int, 3> a)
+{
+    return std::array<int, 3> {
+            a[0],
+            a[1] + a[0],
+            a[2] + a[1]
+        };
+}
+
+constexpr std::array<int, 3> coeff_inc_first(const std::array<int, 3> a)
+{
+    return std::array<int, 3> { a[0] + 1, a[1], a[2] };
+}
+
+constexpr std::array<int, 3>
+coeff_div_by_1_minus_eps(const std::array<int, 3> a)
+{
+    return std::array<int, 3> {
+        a[0],
+        a[0] + a[1],
+        a[0] + a[1] + a[2] + 1
+    };
+}
+
+template <int N>
+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 <typename Expression>
+struct stage_a_error_bound
+    <
+        Expression,
+        stage_a_error_propagation_cases::sum_or_diff
+    >
+{
+private:
+    using leb  = stage_a_error_bound<typename Expression::left>;
+    using reb = stage_a_error_bound<typename Expression::right>;
+    static constexpr auto max_a = coeff_max(leb::a, reb::a);
+public:
+    using magnitude = sum<typename leb::magnitude, typename reb::magnitude>;
+    static constexpr std::array<int, 3> a =
+        coeff_inc_first(coeff_mult_by_1_plus_eps(max_a));
+};
+
+template <typename Expression>
+struct stage_a_error_bound
+    <
+        Expression,
+        stage_a_error_propagation_cases::product
+    >
+{
+private:
+    using leb  = stage_a_error_bound<typename Expression::left>;
+    using reb = stage_a_error_bound<typename Expression::right>;
+    static constexpr std::array<int, 3> a_prod = coeff_product(leb::a, reb::a);
+public:
+    using magnitude =
+        product<typename leb::magnitude, typename reb::magnitude>;
+    static constexpr std::array<int, 3> 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<Expression>
+>
+struct stage_a_condition {};
+
+template <typename Expression>
+struct stage_a_condition
+    <
+        Expression,
+        stage_a_error_propagation_cases::sum_or_diff
+    >
+{
+private:
+    using leb = stage_a_error_bound<typename Expression::left>;
+    using reb = stage_a_error_bound<typename Expression::right>;
+    static constexpr auto max_a = coeff_max(leb::a, reb::a);
+    static constexpr std::array<int, 3> 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<a[0]>();
+    static constexpr int eps_square_coeff =
+        a[2] > 0 ? c * ((a[1] + 1) / c + 1) : c * (a[1] / c + 1);
+public:
+    using magnitude = sum<typename leb::magnitude, typename reb::magnitude>;
+    static constexpr std::array<int, 2> coefficients {a[0], eps_square_coeff};
+};
+
+template <typename Expression>
+struct stage_a_condition
+    <
+        Expression,
+        stage_a_error_propagation_cases::product
+    >
+{
+private:
+    using leb = stage_a_error_bound<typename Expression::left>;
+    using reb = stage_a_error_bound<typename Expression::right>;
+    static constexpr auto a_prod = coeff_product(leb::a, reb::a);
+    static constexpr std::array<int, 3> 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<a[0]>();
+    static constexpr int eps_square_coeff =
+        a[2] > 0 ? c * ((a[1] + 1) / c + 1) : c * (a[1] / c + 1);
+public:
+    using magnitude = product<typename leb::magnitude, typename reb::magnitude>;
+    static constexpr std::array<int, 2> 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<Expression>;
+    static constexpr ct eps = std::numeric_limits<ct>::epsilon() / 2.0;
+    struct constant : public static_constant_interface<ct>
+    {
+        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<constant, typename stage_a_cond::magnitude>;
+};
+
+}} // namespace detail::generic_robust_predicates
+
+}} // namespace boost::geometry
+
+#endif // BOOST_GEOMETRY_EXTENSIONS_GENERIC_ROBUST_PREDICATES_STRATEGIES_CARTESIAN_DETAIL_STAGE_A_ERROR_BOUND_HPP

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 <array>
+#include <iterator>
+
+#include <boost/mp11/list.hpp>
+#include <boost/mp11/algorithm.hpp>
+
+#include <boost/geometry/extensions/generic_robust_predicates/strategies/cartesian/detail/stage_d.hpp>
+
+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 <Expression, true>
+{
+    using type = boost::mp11::mp_bool
+        <
+            Expression::left::is_leaf && Expression::right::is_leaf
+        >;
+};
+
+template <typename Expression>
+using is_leaf_difference = typename is_leaf_difference_impl<Expression>::type;
+
+template
+<
+    typename Evals,
+    typename LeafDifferences,
+    typename AccumulatedSizes,
+    typename Iter,
+    typename CT,
+    std::size_t RemainingDifferences =
+        boost::mp11::mp_size<LeafDifferences>::value
+>
+struct all_differences_zero_tail
+{
+    template <typename ...CTs>
+    static bool apply(Iter begin, Iter end, CTs const&... args)
+    {
+        using eval = boost::mp11::mp_front<LeafDifferences>;
+        using left = typename eval::left;
+        using right = typename eval::right;
+        std::array<CT, sizeof...(CTs)> input
+            {{ static_cast<CT>(args)... }};
+        CT left_val = input[left::argn - 1];
+        CT right_val = input[right::argn - 1];
+        using eval_index = boost::mp11::mp_find<Evals, eval>;
+        constexpr std::size_t start =
+            boost::mp11::mp_at<AccumulatedSizes, eval_index>::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<LeafDifferences>,
+                    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 <typename ...CTs>
+    static bool apply(Iter, Iter, CTs const&...)
+    {
+        return true;
+    }
+};
+
+template
+<
+    typename Expression,
+    typename CalculationType,
+    template <int> class ZEPolicy = default_zero_elimination_policy,
+    template <int, int> 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 <typename ...CTs>
+    static inline int apply(CTs const&... args)
+    {
+        using stack = typename boost::mp11::mp_unique<post_order<Expression>>;
+        using evals = typename boost::mp11::mp_remove_if<stack, is_leaf>;
+        using sizes_pre = boost::mp11::mp_transform
+            <
+                expansion_size_stage_b,
+                boost::mp11::mp_pop_back<evals>
+            >;
+        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<accumulated_sizes>::value
+                >;
+        result_array results;
+
+        using leaf_differences =
+            boost::mp11::mp_copy_if<evals, is_leaf_difference>;
+
+        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<ZEPolicy, true>
+        >;
+        std::array
+            <
+                typename result_array::iterator,
+                boost::mp11::mp_size<ze_evals>::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

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 <cstddef>
+#include <array>
+#include <algorithm>
+
+#include <boost/mp11/integral.hpp>
+#include <boost/mp11/list.hpp>
+#include <boost/mp11/function.hpp>
+#include <boost/mp11/algorithm.hpp>
+
+#include <boost/geometry/extensions/generic_robust_predicates/strategies/cartesian/detail/expression_tree.hpp>
+#include <boost/geometry/extensions/generic_robust_predicates/strategies/cartesian/detail/expansion_eval.hpp>
+
+namespace boost { namespace geometry
+{
+
+namespace detail { namespace generic_robust_predicates
+{
+
+template
+<
+    typename Expression,
+    typename CalculationType,
+    template <int> class ZEPolicy = default_zero_elimination_policy,
+    template <int, int> 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 <typename ...CTs>
+    static inline int apply(CTs const&... args)
+    {
+        using stack = typename boost::mp11::mp_unique<post_order<Expression>>;
+        using evals = typename boost::mp11::mp_remove_if<stack, is_leaf>;
+        using sizes_pre = boost::mp11::mp_transform
+            <
+                expansion_size,
+                boost::mp11::mp_pop_back<evals>
+            >;
+        using sizes = boost::mp11::mp_push_back
+            <
+                sizes_pre,
+                boost::mp11::mp_size_t
+                    <
+                        final_expansion_size
+                            <
+                                Expression,
+                                expansion_size<typename Expression::left>::value,
+                                expansion_size<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<accumulated_sizes>::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

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 <type_traits>
+#include <array>
+#include <tuple>
+
+#include <boost/mp11/list.hpp>
+#include <boost/mp11/algorithm.hpp>
+#include <boost/mp11/set.hpp>
+
+#include <boost/geometry/extensions/generic_robust_predicates/strategies/cartesian/detail/expression_tree.hpp>
+
+namespace boost { namespace geometry
+{
+
+namespace detail { namespace generic_robust_predicates
+{
+
+
+template <typename RemainingStages>
+constexpr bool has_next_and_is_stateful =
+    boost::mp11::mp_front<RemainingStages>::stateful;
+
+template <>
+constexpr bool has_next_and_is_stateful<boost::mp11::mp_list<>> = false;
+
+template
+<
+    typename StatefulStages,
+    typename RemainingStages,
+    bool next_is_stateful = has_next_and_is_stateful<RemainingStages>
+>
+struct next_stage
+{
+    template <typename ...CTs>
+    static inline int apply(StatefulStages const& stages,
+                            CTs const&... args)
+    {
+        using stage = boost::mp11::mp_front<RemainingStages>;
+        int sign = stage::template apply<>(args...);
+        if (sign == sign_uncertain)
+        {
+            return next_stage
+                <
+                    StatefulStages,
+                    boost::mp11::mp_pop_front<RemainingStages>
+                >::apply(stages, args...);
+        }
+        else
+        {
+            return sign;
+        }
+    }
+};
+
+template
+<
+    typename StatefulStages,
+    typename RemainingStages
+>
+struct next_stage
+    <
+        StatefulStages,
+        RemainingStages,
+        true
+    >
+{
+    template <typename ...CTs>
+    static inline int apply(StatefulStages const& stages,
+                            CTs const&... args)
+    {
+        using stage = boost::mp11::mp_front<RemainingStages>;
+        int sign = std::get<stage>(stages).apply(args...);
+        if (sign == sign_uncertain)
+        {
+            return next_stage
+                <
+                    StatefulStages,
+                    boost::mp11::mp_pop_front<RemainingStages>
+                >::apply(stages, args...);
+        }
+        else
+        {
+            return sign;
+        }
+    }
+};
+
+template
+<
+    typename StatefulStages
+>
+struct next_stage
+    <
+        StatefulStages,
+        boost::mp11::mp_list<>,
+        false
+    >
+{
+    template <typename ...CTs>
+    static inline int apply(StatefulStages const&,
+                            CTs const&...)
+    {
+        return sign_uncertain;
+    }
+};
+
+template <typename Stage> using is_stateful = boost::mp11::mp_bool<Stage::stateful>;
+
+template <typename Stage, bool updates = Stage::updates>
+struct construct_stage_impl
+{
+    template <typename Array> static inline Stage apply(Array const& a)
+    {
+        Stage out(a);
+        return out;
+    }
+};
+
+template <typename Stage>
+struct construct_stage_impl<Stage, true>
+{
+    template <typename Array> static inline Stage apply(Array const&)
+    {
+        Stage out;
+        return out;
+    }
+};
+
+template <typename Stages>
+struct construct_stages_impl
+{
+    template <typename Array>
+    static inline boost::mp11::mp_rename<Stages, std::tuple>
+    apply(Array const& a)
+    {
+        using stage = boost::mp11::mp_front<Stages>;
+        std::tuple<stage> first(construct_stage_impl<stage>::apply(a));
+        return std::tuple_cat(first,
+                              construct_stages_impl
+                                    <
+                                        boost::mp11::mp_pop_front<Stages>
+                                    >::apply(a));
+    }
+};
+
+template <>
+struct construct_stages_impl<std::tuple<>>
+{
+    template <typename Array>
+    static inline std::tuple<> apply(Array const&)
+    {
+        return std::tuple<>{};
+    }
+};
+
+template <typename Stage>
+using is_updatable = boost::mp11::mp_bool<Stage::updates>;
+
+template <typename UpdatableStages>
+struct update_all_impl
+{
+    template <typename Stages, typename ...CTs>
+    static void apply(Stages& stages, CTs const&... args)
+    {
+        using current_stage = boost::mp11::mp_front<UpdatableStages>;
+        std::get<current_stage>(stages).update(args...);
+        update_all_impl<boost::mp11::mp_pop_front<UpdatableStages>>
+            ::apply(stages, args...);
+    }
+};
+
+template <>
+struct update_all_impl<std::tuple<>>
+{
+    template <typename Stages, typename ...CTs>
+    static void apply(Stages&, CTs const&...) {}
+};
+
+template
+<
+    typename CalculationType,
+    typename ...Stages
+>
+struct staged_predicate
+{
+private:
+    using ct = CalculationType;
+    using stages = std::tuple<Stages...>;
+    using stateful_stages = boost::mp11::mp_copy_if<stages, is_stateful>;
+    using updatable_stages = boost::mp11::mp_copy_if<stages, is_updatable>;
+    stateful_stages m_stages;
+    using all_stages = boost::mp11::mp_list<Stages...>;
+public:
+    static constexpr bool stateful =
+        boost::mp11::mp_any_of<stages, is_stateful>::value;
+    static constexpr bool updates =
+        boost::mp11::mp_any_of<stages, is_updatable>::value;
+    template <typename ...CTs>
+    inline staged_predicate(CTs const&... args) : m_stages(
+            construct_stages_impl<stateful_stages>::apply(
+                std::array<ct, sizeof...(args)>{ static_cast<ct>(args)... }
+            )) {}
+
+    template <typename ...CTs>
+    inline void update(CTs const&... args)
+    {
+        update_all_impl<updatable_stages>::apply(m_stages, args...);
+    }
+
+    template <typename ...CTs>
+    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

include/boost/geometry/extensions/generic_robust_predicates/strategies/cartesian/detail/static_filter.hpp

@@ -0,0 +1,109 @@
+// Boost.Geometry (aka GGL, Generic Geometry Library)
+
+// Copyright (c) 2020 Tinko Bartels, Berlin, Germany.
+
+// Contributed and/or modified by Tinko Bartels,
+//   as part of Google Summer of Code 2020 program.
+
+// Use, modification and distribution is subject to the Boost Software License,
+// Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at
+// http://www.boost.org/LICENSE_1_0.txt)
+
+#ifndef BOOST_GEOMETRY_EXTENSIONS_GENERIC_ROBUST_PREDICATES_STRATEGIES_CARTESIAN_DETAIL_STATIC_FILTER_HPP
+#define BOOST_GEOMETRY_EXTENSIONS_GENERIC_ROBUST_PREDICATES_STRATEGIES_CARTESIAN_DETAIL_STATIC_FILTER_HPP
+
+#include <limits>
+#include <array>
+
+#include <boost/mp11/list.hpp>
+#include <boost/mp11/algorithm.hpp>
+
+#include <boost/geometry/extensions/generic_robust_predicates/strategies/cartesian/detail/expression_tree.hpp>
+#include <boost/geometry/extensions/generic_robust_predicates/strategies/cartesian/detail/expression_eval.hpp>
+
+namespace boost { namespace geometry
+{
+
+namespace detail { namespace generic_robust_predicates
+{
+
+//The static filter works similar to the semi static filter with the exception
+//that the error bound is only computed once at construction time. For this
+//purpose, the filter is stateful and it requires upper and lower bounds on the
+//inputs at compile time like this:
+//
+//static_filter<...> f(max1, max2, ..., min1, min2, ...);
+//
+//for each call to the filter, i.e. f.apply(arg1, arg2, ...); it must hold that
+//min1 <= arg1 <= max1 and so on.
+
+template
+<
+    typename Expression,
+    typename CalculationType,
+    typename ErrorExpression
+>
+class static_filter
+{
+private:
+    using ct = CalculationType;
+    ct m_error_bound;
+    static constexpr std::size_t extrema_count = max_argn<ErrorExpression>;
+public:
+    static constexpr bool stateful = true;
+    static constexpr bool updates = false;
+
+    inline static_filter() : m_error_bound(std::numeric_limits<ct>::max()) {}
+
+    inline ct error_bound() const { return m_error_bound; }
+
+    template <typename ...CTs>
+    inline static_filter(CTs const&... args)
+        : m_error_bound(evaluate_expression<ErrorExpression>(
+                std::array<ct, extrema_count>
+                    {static_cast<ct>(args)...}))
+    {
+        static_assert(sizeof...(CTs) == extrema_count,
+                      "Number of constructor arguments is incompatible with error expression.");
+    }
+
+    inline static_filter(std::array<ct, extrema_count> const& extrema)
+        : m_error_bound(evaluate_expression<ErrorExpression>(extrema)) {}
+
+    inline static_filter(ct const& b) : m_error_bound(b) {}
+
+    template <typename ...CTs>
+    inline int apply(CTs const&... args) const
+    {
+        std::array<ct, sizeof...(CTs)> input {static_cast<ct>(args)...};
+        const ct det = evaluate_expression<Expression>(input);
+        if (det > m_error_bound)
+        {
+            return 1;
+        }
+        else if (det < -m_error_bound)
+        {
+            return -1;
+        }
+        else if (m_error_bound == 0 && det == 0)
+        {
+            return 0;
+        }
+        else
+        {
+            return sign_uncertain;
+        }
+    }
+
+    template <typename ...CTs>
+    inline int operator()(CTs const&... args) const
+    {
+        return apply(args...);
+    }
+};
+
+}} // namespace detail::generic_robust_predicates
+
+}} // namespace boost::geometry
+
+#endif // BOOST_GEOMETRY_EXTENSIONS_GENERIC_ROBUST_PREDICATES_STRATEGIES_CARTESIAN_DETAIL_STATIC_FILTER_HPP