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 <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(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<double> 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 : : : <debug-symbols>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 <array>
+
+#include <geometry_test_common.hpp>
+
+#include <boost/geometry/extensions/generic_robust_predicates/strategies/cartesian/detail/approximate.hpp>
+#include <boost/geometry/extensions/generic_robust_predicates/strategies/cartesian/detail/expression_tree.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::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<sum<_1, _2>, ct>(
+            std::array<ct, 2>{1.0, 2.0});
+    BOOST_CHECK_EQUAL(3.0, r1);
+    ct r2 = approximate_value<max<abs<_1>, abs<_2>>, ct>(
+            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};
+    approximate_interim<evals, evals, ct>(r, input);
+    ct r3 = get_approx<evals, typename expression::left, ct>(r, input);
+    BOOST_CHECK_EQUAL(2.0, r3);
+    ct r4 = get_approx<evals, typename expression::right, ct>(r, input);
+    BOOST_CHECK_EQUAL(-6.0, r4);
+    ct r5 = get_approx<evals, expression, ct>(r, input);
+    BOOST_CHECK_EQUAL(-12.0, r5);
+}
+
+
+int test_main(int, char* [])
+{
+    test_all<double>();
+    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 <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;
+}
diff --git a/include/boost/geometry/algorithms/detail/buffer/buffer_inserter.hpp b/include/boost/geometry/algorithms/detail/buffer/buffer_inserter.hpp
index 9d1faee770..976f011e26 100644
--- a/include/boost/geometry/algorithms/detail/buffer/buffer_inserter.hpp
+++ b/include/boost/geometry/algorithms/detail/buffer/buffer_inserter.hpp
@@ -1,6 +1,6 @@
 // Boost.Geometry (aka GGL, Generic Geometry Library)
 
-// Copyright (c) 2012-2014 Barend Gehrels, Amsterdam, the Netherlands.
+// Copyright (c) 2012-2020 Barend Gehrels, Amsterdam, the Netherlands.
 
 // This file was modified by Oracle on 2017.
 // Modifications copyright (c) 2017 Oracle and/or its affiliates.
@@ -31,16 +31,13 @@
 
 #include <boost/geometry/strategies/buffer.hpp>
 #include <boost/geometry/strategies/side.hpp>
-#include <boost/geometry/algorithms/detail/make/make.hpp>
+#include <boost/geometry/algorithms/detail/direction_code.hpp>
 #include <boost/geometry/algorithms/detail/buffer/buffered_piece_collection.hpp>
 #include <boost/geometry/algorithms/detail/buffer/line_line_intersection.hpp>
 
-#include <boost/geometry/algorithms/assign.hpp>
 #include <boost/geometry/algorithms/num_interior_rings.hpp>
 #include <boost/geometry/algorithms/simplify.hpp>
 
-#include <boost/geometry/arithmetic/infinite_line_functions.hpp>
-
 #include <boost/geometry/views/detail/normalized_view.hpp>
 
 
@@ -113,23 +110,11 @@ struct buffer_range
             JoinStrategy const& join_strategy,
             EndStrategy const& end_strategy,
             RobustPolicy const& ,
-            SideStrategy const& side_strategy) // side strategy
+            SideStrategy const& side_strategy)
     {
-        output_point_type intersection_point;
-        geometry::assign_zero(intersection_point);
-
-        geometry::strategy::buffer::join_selector join
-                = get_join_type(penultimate_input, previous_input, input, side_strategy);
-        if (join == geometry::strategy::buffer::join_convex)
-        {
-            // Calculate the intersection-point formed by the two sides.
-            // It might be that the two sides are not convex, but continue
-            // or spikey, we then change the join-type
-            join = line_line_intersection::apply(
-                        perp1, perp2, prev_perp1, prev_perp2,
-                        intersection_point);
-
-        }
+        geometry::strategy::buffer::join_selector const join
+                = get_join_type(penultimate_input, previous_input, input,
+                                side_strategy);
 
         switch(join)
         {
@@ -163,6 +148,11 @@ struct buffer_range
                 {
                     // The corner is convex, we create a join
                     // TODO (future) - avoid a separate vector, add the piece directly
+                    output_point_type const
+                        intersection_point
+                            = line_line_intersection::apply(perp1, perp2,
+                                    prev_perp1, prev_perp2);
+
                     std::vector<output_point_type> range_out;
                     if (join_strategy.apply(intersection_point,
                                 previous_input, prev_perp2, perp1,
@@ -177,14 +167,14 @@ struct buffer_range
         }
     }
 
-    static inline bool similar_direction(output_point_type const& p0,
+    // Returns true if collinear point p2 continues after p0 and p1.
+    // If it turns back (spike), it returns false.
+    static inline bool same_direction(output_point_type const& p0,
             output_point_type const& p1,
             output_point_type const& p2)
     {
-        typedef model::infinite_line<coordinate_type> line_type;
-        line_type const p = detail::make::make_infinite_line<coordinate_type>(p0, p1);
-        line_type const q = detail::make::make_infinite_line<coordinate_type>(p1, p2);
-        return arithmetic::similar_direction(p, q);
+        typedef typename cs_tag<output_point_type>::type cs_tag;
+        return direction_code<cs_tag>(p0, p1, p2) == 1;
     }
 
     template <typename SideStrategy>
@@ -197,8 +187,7 @@ struct buffer_range
         int const side = side_strategy.apply(p0, p1, p2);
         return side == -1 ? geometry::strategy::buffer::join_convex
             :  side == 1  ? geometry::strategy::buffer::join_concave
-            :  similar_direction(p0, p1, p2)
-                          ? geometry::strategy::buffer::join_continue
+            :  same_direction(p0, p1, p2) ? geometry::strategy::buffer::join_continue
             : geometry::strategy::buffer::join_spike;
     }
 
diff --git a/include/boost/geometry/algorithms/detail/buffer/buffer_policies.hpp b/include/boost/geometry/algorithms/detail/buffer/buffer_policies.hpp
index 4d071dc1d3..6e10bc3bd3 100644
--- a/include/boost/geometry/algorithms/detail/buffer/buffer_policies.hpp
+++ b/include/boost/geometry/algorithms/detail/buffer/buffer_policies.hpp
@@ -139,8 +139,8 @@ struct buffer_turn_operation
     {}
 };
 
-// Version for buffer including type of location, is_opposite, and helper variables
-template <typename Point, typename RobustPoint, typename SegmentRatio>
+// Version of turn_info for buffer with its turn index and other helper variables
+template <typename Point, typename SegmentRatio>
 struct buffer_turn_info
     : public detail::overlay::turn_info
         <
@@ -153,9 +153,6 @@ struct buffer_turn_info
 
     std::size_t turn_index;
 
-    // Still necessary for turn-in-original calculations
-    RobustPoint robust_point;
-
     // Information if turn can be used. It is not traversable if it is within
     // another piece, or within the original (depending on deflation),
     // or (for deflate) if there are not enough points to traverse it.
diff --git a/include/boost/geometry/algorithms/detail/buffer/buffered_piece_collection.hpp b/include/boost/geometry/algorithms/detail/buffer/buffered_piece_collection.hpp
index 8bc7e370ff..438921be52 100644
--- a/include/boost/geometry/algorithms/detail/buffer/buffered_piece_collection.hpp
+++ b/include/boost/geometry/algorithms/detail/buffer/buffered_piece_collection.hpp
@@ -123,16 +123,11 @@ struct buffered_piece_collection
 {
     typedef typename geometry::point_type<Ring>::type point_type;
     typedef typename geometry::coordinate_type<Ring>::type coordinate_type;
-    typedef typename geometry::robust_point_type
-    <
-        point_type,
-        RobustPolicy
-    >::type robust_point_type;
 
     // Robust ring/polygon type, always clockwise
-    typedef geometry::model::ring<robust_point_type> robust_ring_type;
+    typedef geometry::model::ring<point_type> clockwise_ring_type;
 
-    typedef geometry::model::box<point_type> normal_box_type;
+    typedef geometry::model::box<point_type> box_type;
 
     typedef typename IntersectionStrategy::side_strategy_type side_strategy_type;
     typedef typename IntersectionStrategy::envelope_strategy_type envelope_strategy_type;
@@ -150,32 +145,21 @@ struct buffered_piece_collection
 
     typedef typename IntersectionStrategy::template point_in_geometry_strategy
         <
-            robust_point_type,
-            robust_ring_type
+            point_type,
+            clockwise_ring_type
         >::type point_in_geometry_strategy_type;
 
-    typedef typename geometry::rescale_policy_type
-        <
-            typename geometry::point_type<Ring>::type,
-            typename IntersectionStrategy::cs_tag
-        >::type rescale_policy_type;
-
-    typedef geometry::segment_ratio
-    <
-        typename geometry::coordinate_type<robust_point_type>::type
-    > ratio_type;
 
     typedef buffer_turn_info
     <
         point_type,
-        robust_point_type,
-        ratio_type
+        typename segment_ratio_type<point_type, RobustPolicy>::type
     > buffer_turn_info_type;
 
     typedef buffer_turn_operation
     <
         point_type,
-        ratio_type
+        typename segment_ratio_type<point_type, RobustPolicy>::type
     > buffer_turn_operation_type;
 
     typedef std::vector<buffer_turn_info_type> turn_vector_type;
@@ -237,8 +221,7 @@ struct buffered_piece_collection
 
     struct original_ring
     {
-        typedef geometry::model::box<robust_point_type> robust_box_type;
-        typedef geometry::sections<robust_box_type, 1> sections_type;
+        typedef geometry::sections<box_type, 1> sections_type;
 
         // Creates an empty instance
         inline original_ring()
@@ -246,7 +229,7 @@ struct buffered_piece_collection
             , m_has_interiors(false)
         {}
 
-        inline original_ring(robust_ring_type const& ring,
+        inline original_ring(clockwise_ring_type const& ring,
                 bool is_interior, bool has_interiors,
                 envelope_strategy_type const& envelope_strategy,
                 expand_strategy_type const& expand_strategy)
@@ -266,8 +249,8 @@ struct buffered_piece_collection
                     envelope_strategy, expand_strategy);
         }
 
-        robust_ring_type m_ring;
-        robust_box_type m_box;
+        clockwise_ring_type m_ring;
+        box_type m_box;
         sections_type m_sections;
 
         bool m_is_interior;
@@ -293,7 +276,7 @@ struct buffered_piece_collection
 
     // Specificly for offsetted rings around points
     // but also for large joins with many points
-    typedef geometry::sections<normal_box_type, 2> sections_type;
+    typedef geometry::sections<box_type, 2> sections_type;
     sections_type monotonic_sections;
 
     // Define the clusters, mapping cluster_id -> turns
@@ -329,8 +312,7 @@ struct buffered_piece_collection
         , m_envelope_strategy(intersection_strategy.get_envelope_strategy())
         , m_expand_strategy(intersection_strategy.get_expand_strategy())
         , m_point_in_geometry_strategy(intersection_strategy
-            .template get_point_in_geometry_strategy<robust_point_type,
-                        robust_ring_type>())
+            .template get_point_in_geometry_strategy<point_type, clockwise_ring_type>())
         , m_robust_policy(robust_policy)
     {}
 
@@ -496,7 +478,7 @@ struct buffered_piece_collection
 
         geometry::partition
             <
-                normal_box_type,
+                box_type,
                 include_turn_policy,
                 detail::partition::include_all_policy
             >::apply(m_turns, original_rings, visitor,
@@ -534,9 +516,6 @@ struct buffered_piece_collection
         {
             buffer_turn_info_type& turn = *it;
 
-            // Calculate rescaled turn points
-            geometry::recalculate(turn.robust_point, turn.point, m_robust_policy);
-
             // Update member used
             turn.turn_index = index;
 
@@ -612,7 +591,7 @@ struct buffered_piece_collection
                                                    m_envelope_strategy);
             geometry::partition
                 <
-                    normal_box_type
+                    box_type
                 >::apply(monotonic_sections, visitor,
                          get_section_box_type(),
                          overlaps_section_box_type());
@@ -639,7 +618,7 @@ struct buffered_piece_collection
 
             geometry::partition
                 <
-                    normal_box_type
+                    box_type
                 >::apply(m_turns, m_pieces, visitor,
                          turn_get_box(), turn_ovelaps_box_type(),
                          piece_get_box(), piece_ovelaps_box_type());
@@ -771,7 +750,7 @@ struct buffered_piece_collection
             // For rescaling, it is recalculated. Without rescaling, it
             // is just assigning (note that this Ring type is the
             // GeometryOut type, which might differ from the input ring type)
-            geometry::model::ring<robust_point_type> adapted_ring;
+            clockwise_ring_type clockwise_ring;
 
             typedef detail::normalized_view<InputRing const> view_type;
             view_type const view(input_ring);
@@ -779,13 +758,11 @@ struct buffered_piece_collection
             for (typename boost::range_iterator<view_type const>::type it =
                 boost::begin(view); it != boost::end(view); ++it)
             {
-                robust_point_type adapted_point;
-                geometry::recalculate(adapted_point, *it, m_robust_policy);
-                adapted_ring.push_back(adapted_point);
+                clockwise_ring.push_back(*it);
             }
 
             original_rings.back()
-                = original_ring(adapted_ring,
+                = original_ring(clockwise_ring,
                     is_interior, has_interiors,
                     m_envelope_strategy, m_expand_strategy);
         }
@@ -1085,9 +1062,6 @@ struct buffered_piece_collection
     {
         typedef typename IntersectionStrategy::disjoint_point_box_strategy_type d_pb_strategy_type;
 
-        robust_point_type any_point;
-        geometry::recalculate(any_point, point, m_robust_policy);
-
         signed_size_type count_in_original = 0;
 
         // Check of the robust point of this outputted ring is in
@@ -1104,7 +1078,7 @@ struct buffered_piece_collection
             {
                 continue;
             }
-            if (detail::disjoint::disjoint_point_box(any_point,
+            if (detail::disjoint::disjoint_point_box(point,
                                                      original.m_box,
                                                      d_pb_strategy_type()))
             {
@@ -1112,7 +1086,7 @@ struct buffered_piece_collection
             }
 
             int const geometry_code
-                = detail::within::point_in_geometry(any_point,
+                = detail::within::point_in_geometry(point,
                     original.m_ring, m_point_in_geometry_strategy);
 
             if (geometry_code == -1)
diff --git a/include/boost/geometry/algorithms/detail/buffer/line_line_intersection.hpp b/include/boost/geometry/algorithms/detail/buffer/line_line_intersection.hpp
index ec50d91a8f..45adb5abd2 100644
--- a/include/boost/geometry/algorithms/detail/buffer/line_line_intersection.hpp
+++ b/include/boost/geometry/algorithms/detail/buffer/line_line_intersection.hpp
@@ -1,6 +1,6 @@
 // Boost.Geometry (aka GGL, Generic Geometry Library)
 
-// Copyright (c) 2012-2019 Barend Gehrels, Amsterdam, the Netherlands.
+// Copyright (c) 2012-2020 Barend Gehrels, Amsterdam, the Netherlands.
 
 // Use, modification and distribution is subject to the Boost Software License,
 // Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at
@@ -9,10 +9,12 @@
 #ifndef BOOST_GEOMETRY_ALGORITHMS_DETAIL_BUFFER_LINE_LINE_INTERSECTION_HPP
 #define BOOST_GEOMETRY_ALGORITHMS_DETAIL_BUFFER_LINE_LINE_INTERSECTION_HPP
 
-
+#include <boost/geometry/core/assert.hpp>
 #include <boost/geometry/arithmetic/infinite_line_functions.hpp>
 #include <boost/geometry/algorithms/detail/make/make.hpp>
-#include <boost/geometry/strategies/buffer.hpp>
+
+#include <boost/core/ignore_unused.hpp>
+
 
 namespace boost { namespace geometry
 {
@@ -26,10 +28,8 @@ namespace detail { namespace buffer
 struct line_line_intersection
 {
     template <typename Point>
-    static inline strategy::buffer::join_selector
-    apply(Point const& pi, Point const& pj,
-          Point const& qi, Point const& qj,
-          Point& ip)
+    static inline Point
+    apply(Point const& pi, Point const& pj, Point const& qi, Point const& qj)
     {
         typedef typename coordinate_type<Point>::type ct;
         typedef model::infinite_line<ct> line_type;
@@ -37,18 +37,14 @@ struct line_line_intersection
         line_type const p = detail::make::make_infinite_line<ct>(pi, pj);
         line_type const q = detail::make::make_infinite_line<ct>(qi, qj);
 
-        if (arithmetic::intersection_point(p, q, ip))
-        {
-            return strategy::buffer::join_convex;
-        }
-
-        // The lines do not intersect.
-        // Distinguish between continuing lines (having a similar direction)
-        // and spikes (having the opposite direction).
-        return arithmetic::similar_direction(p, q)
-            ? strategy::buffer::join_continue
-            : strategy::buffer::join_spike
-            ;
+        // The input lines are not parallel, they intersect, because
+        // their join type is checked before.
+        Point ip;
+        bool const intersecting = arithmetic::intersection_point(p, q, ip);
+        BOOST_GEOMETRY_ASSERT(intersecting);
+        boost::ignore_unused(intersecting);
+
+        return ip;
     }
 };
 
diff --git a/include/boost/geometry/algorithms/detail/buffer/turn_in_original_visitor.hpp b/include/boost/geometry/algorithms/detail/buffer/turn_in_original_visitor.hpp
index 4e20f9662a..ab9ab08daa 100644
--- a/include/boost/geometry/algorithms/detail/buffer/turn_in_original_visitor.hpp
+++ b/include/boost/geometry/algorithms/detail/buffer/turn_in_original_visitor.hpp
@@ -35,7 +35,12 @@ struct original_get_box
     template <typename Box, typename Original>
     static inline void apply(Box& total, Original const& original)
     {
-        geometry::expand(total, original.m_box);
+        typedef typename strategy::expand::services::default_strategy
+            <
+                box_tag, typename cs_tag<Box>::type
+            >::type expand_strategy_type;
+
+        geometry::expand(total, original.m_box, expand_strategy_type());
     }
 };
 
@@ -72,7 +77,7 @@ struct turn_in_original_ovelaps_box
         }
 
         return ! geometry::detail::disjoint::disjoint_point_box(
-                    turn.robust_point, box, DisjointPointBoxStrategy());
+                    turn.point, box, DisjointPointBoxStrategy());
     }
 };
 
@@ -226,13 +231,13 @@ class turn_in_original_visitor
             return true;
         }
 
-        if (geometry::disjoint(turn.robust_point, original.m_box))
+        if (geometry::disjoint(turn.point, original.m_box))
         {
             // Skip all disjoint
             return true;
         }
 
-        int const code = point_in_original(turn.robust_point, original, m_point_in_geometry_strategy);
+        int const code = point_in_original(turn.point, original, m_point_in_geometry_strategy);
 
         if (code == -1)
         {
diff --git a/include/boost/geometry/algorithms/detail/direction_code.hpp b/include/boost/geometry/algorithms/detail/direction_code.hpp
index 38c65afa2e..0e81a202a0 100644
--- a/include/boost/geometry/algorithms/detail/direction_code.hpp
+++ b/include/boost/geometry/algorithms/detail/direction_code.hpp
@@ -1,6 +1,6 @@
 // Boost.Geometry (aka GGL, Generic Geometry Library)
 
-// Copyright (c) 2015 Barend Gehrels, Amsterdam, the Netherlands.
+// Copyright (c) 2015-2020 Barend Gehrels, Amsterdam, the Netherlands.
 
 // This file was modified by Oracle on 2015, 2017, 2019.
 // Modifications copyright (c) 2015-2019 Oracle and/or its affiliates.
@@ -54,21 +54,29 @@ struct direction_code_impl<cartesian_tag>
 
         typedef model::infinite_line<calc_t> line_type;
 
-        // point b is often equal to the specified point, check that first
-        line_type const q = detail::make::make_infinite_line<calc_t>(segment_b, point);
-        if (arithmetic::is_degenerate(q))
+        // Situation and construction of perpendicular line
+        //
+        //     P1     a--------------->b   P2
+        //                             |
+        //                             |
+        //                             v
+        //
+        // P1 is located right of the (directional) perpendicular line
+        // and therefore gets a negative side_value, and returns -1.
+        // P2 is to the left of the perpendicular line and returns 1.
+        // If the specified point is located on top of b, it returns 0.
+
+        line_type const line
+            = detail::make::make_perpendicular_line<calc_t>(segment_a,
+                segment_b, segment_b);
+
+        if (arithmetic::is_degenerate(line))
         {
             return 0;
         }
 
-        line_type const p = detail::make::make_infinite_line<calc_t>(segment_a, segment_b);
-        if (arithmetic::is_degenerate(p))
-        {
-            return 0;
-        }
-
-        // p extends a-b if direction is similar
-        return arithmetic::similar_direction(p, q) ? 1 : -1;
+        calc_t const sv = arithmetic::side_value(line, point);
+        return sv == 0 ? 0 : sv > 0 ? 1 : -1;
     }
 };
 
diff --git a/include/boost/geometry/algorithms/detail/intersection/areal_areal.hpp b/include/boost/geometry/algorithms/detail/intersection/areal_areal.hpp
index a9bf98fabe..994c3c7daa 100644
--- a/include/boost/geometry/algorithms/detail/intersection/areal_areal.hpp
+++ b/include/boost/geometry/algorithms/detail/intersection/areal_areal.hpp
@@ -26,7 +26,7 @@ namespace detail { namespace intersection
 template
 <
     typename GeometryOut,
-    typename OutTag = typename detail::intersection::tag
+    typename OutTag = typename geometry::detail::setop_insert_output_tag
                         <
                             typename geometry::detail::output_geometry_value
                                 <
@@ -85,9 +85,10 @@ struct intersection_areal_areal_<TupledOut, tupled_output_tag>
 
         boost::ignore_unused
             <
-                detail::intersection::expect_output_pla
+                geometry::detail::expect_output
                     <
-                        Areal1, Areal2, single_out
+                        Areal1, Areal2, single_out,
+                        point_tag, linestring_tag, polygon_tag
                     >
             >();
 
diff --git a/include/boost/geometry/algorithms/detail/intersection/multi.hpp b/include/boost/geometry/algorithms/detail/intersection/multi.hpp
index 6a8c1f6c62..3b2b1461ca 100644
--- a/include/boost/geometry/algorithms/detail/intersection/multi.hpp
+++ b/include/boost/geometry/algorithms/detail/intersection/multi.hpp
@@ -432,8 +432,8 @@ struct intersection_insert
         TupledOut,
         OverlayType,
         ReverseMultiLinestring, ReverseRing,
-        multi_linestring_tag, ring_tag, detail::intersection::tupled_output_tag,
-        linear_tag, areal_tag, detail::intersection::tupled_output_tag
+        multi_linestring_tag, ring_tag, detail::tupled_output_tag,
+        linear_tag, areal_tag, detail::tupled_output_tag
     > : detail::intersection::intersection_of_multi_linestring_with_areal
             <
                 ReverseRing,
@@ -441,7 +441,19 @@ struct intersection_insert
                 OverlayType,
                 true
             >
-      , detail::intersection::expect_output_pl<MultiLinestring, Ring, TupledOut>
+      , detail::expect_output
+            <
+                MultiLinestring, Ring, TupledOut,
+                // NOTE: points can be the result only in case of intersection.
+                // TODO: union should require L and A
+                typename boost::mpl::if_c
+                    <
+                        (OverlayType == overlay_intersection),
+                        point_tag,
+                        void
+                    >::type,
+                linestring_tag
+            >
 {};
 
 
@@ -458,8 +470,8 @@ struct intersection_insert
         TupledOut,
         OverlayType,
         ReverseMultiLinestring, ReversePolygon,
-        multi_linestring_tag, polygon_tag, detail::intersection::tupled_output_tag,
-        linear_tag, areal_tag, detail::intersection::tupled_output_tag
+        multi_linestring_tag, polygon_tag, detail::tupled_output_tag,
+        linear_tag, areal_tag, detail::tupled_output_tag
     > : detail::intersection::intersection_of_multi_linestring_with_areal
             <
                 ReversePolygon,
@@ -467,7 +479,19 @@ struct intersection_insert
                 OverlayType,
                 true
             >
-    , detail::intersection::expect_output_pl<MultiLinestring, Polygon, TupledOut>
+      , detail::expect_output
+            <
+                MultiLinestring, Polygon, TupledOut,
+                // NOTE: points can be the result only in case of intersection.
+                // TODO: union should require L and A
+                typename boost::mpl::if_c
+                    <
+                        (OverlayType == overlay_intersection),
+                        point_tag,
+                        void
+                    >::type,
+                linestring_tag
+            >
 {};
 
 template
@@ -483,8 +507,8 @@ struct intersection_insert
         TupledOut,
         OverlayType,
         ReversePolygon, ReverseMultiLinestring,
-        polygon_tag, multi_linestring_tag, detail::intersection::tupled_output_tag,
-        areal_tag, linear_tag, detail::intersection::tupled_output_tag
+        polygon_tag, multi_linestring_tag, detail::tupled_output_tag,
+        areal_tag, linear_tag, detail::tupled_output_tag
     > : detail::intersection::intersection_of_areal_with_multi_linestring
             <
                 ReversePolygon,
@@ -492,7 +516,22 @@ struct intersection_insert
                 OverlayType,
                 true
             >
-    , detail::intersection::expect_output_pl<Polygon, MultiLinestring, TupledOut>
+      , detail::expect_output
+            <
+                Polygon, MultiLinestring, TupledOut,
+                // NOTE: points can be the result only in case of intersection.
+                // TODO: union should require L and A
+                // TODO: in general the result of difference should depend on the first argument
+                //       but this specialization calls L/A in reality so the first argument is linear.
+                //       So expect only L for difference?
+                typename boost::mpl::if_c
+                    <
+                        (OverlayType == overlay_intersection),
+                        point_tag,
+                        void
+                    >::type,
+                linestring_tag
+            >
 {};
 
 template
@@ -508,8 +547,8 @@ struct intersection_insert
         TupledOut,
         OverlayType,
         ReverseMultiLinestring, ReverseMultiPolygon,
-        multi_linestring_tag, multi_polygon_tag, detail::intersection::tupled_output_tag,
-        linear_tag, areal_tag, detail::intersection::tupled_output_tag
+        multi_linestring_tag, multi_polygon_tag, detail::tupled_output_tag,
+        linear_tag, areal_tag, detail::tupled_output_tag
     > : detail::intersection::intersection_of_multi_linestring_with_areal
             <
                 ReverseMultiPolygon,
@@ -517,7 +556,19 @@ struct intersection_insert
                 OverlayType,
                 true
             >
-    , detail::intersection::expect_output_pl<MultiLinestring, MultiPolygon, TupledOut>
+      , detail::expect_output
+            <
+                MultiLinestring, MultiPolygon, TupledOut,
+                // NOTE: points can be the result only in case of intersection.
+                // TODO: union should require L and A
+                typename boost::mpl::if_c
+                    <
+                        (OverlayType == overlay_intersection),
+                        point_tag,
+                        void
+                    >::type,
+                linestring_tag
+            >
 {};
 
 
diff --git a/include/boost/geometry/algorithms/detail/overlay/handle_colocations.hpp b/include/boost/geometry/algorithms/detail/overlay/handle_colocations.hpp
index 951e91706d..c12ddc6481 100644
--- a/include/boost/geometry/algorithms/detail/overlay/handle_colocations.hpp
+++ b/include/boost/geometry/algorithms/detail/overlay/handle_colocations.hpp
@@ -752,6 +752,46 @@ struct is_turn_index
     signed_size_type m_index;
 };
 
+template
+<
+    typename Sbs,
+    typename Point,
+    typename Turns,
+    typename Geometry1,
+    typename Geometry2
+>
+inline bool fill_sbs(Sbs& sbs, Point& turn_point,
+                     cluster_info const& cinfo,
+                     Turns const& turns,
+                     Geometry1 const& geometry1, Geometry2 const& geometry2)
+{
+    typedef typename boost::range_value<Turns>::type turn_type;
+
+    std::set<signed_size_type> const& ids = cinfo.turn_indices;
+
+    if (ids.empty())
+    {
+        return false;
+    }
+
+    bool first = true;
+    for (std::set<signed_size_type>::const_iterator sit = ids.begin();
+         sit != ids.end(); ++sit)
+    {
+        signed_size_type turn_index = *sit;
+        turn_type const& turn = turns[turn_index];
+        if (first )
+        {
+            turn_point = turn.point;
+        }
+        for (int i = 0; i < 2; i++)
+        {
+            sbs.add(turn.operations[i], turn_index, i, geometry1, geometry2, first);
+            first = false;
+        }
+    }
+    return true;
+}
 
 template
 <
@@ -783,32 +823,14 @@ inline void gather_cluster_properties(Clusters& clusters, Turns& turns,
          mit != clusters.end(); ++mit)
     {
         cluster_info& cinfo = mit->second;
-        std::set<signed_size_type> const& ids = cinfo.turn_indices;
-        if (ids.empty())
-        {
-            continue;
-        }
 
         sbs_type sbs(strategy);
         point_type turn_point; // should be all the same for all turns in cluster
-
-        bool first = true;
-        for (std::set<signed_size_type>::const_iterator sit = ids.begin();
-             sit != ids.end(); ++sit)
+        if (! fill_sbs(sbs, turn_point, cinfo, turns, geometry1, geometry2))
         {
-            signed_size_type turn_index = *sit;
-            turn_type const& turn = turns[turn_index];
-            if (first)
-            {
-                turn_point = turn.point;
-            }
-            for (int i = 0; i < 2; i++)
-            {
-                turn_operation_type const& op = turn.operations[i];
-                sbs.add(op, turn_index, i, geometry1, geometry2, first);
-                first = false;
-            }
+            continue;
         }
+
         sbs.apply(turn_point);
 
         sbs.find_open();
@@ -823,7 +845,7 @@ inline void gather_cluster_properties(Clusters& clusters, Turns& turns,
         // polygons
         for (std::size_t i = 0; i < sbs.m_ranked_points.size(); i++)
         {
-            const typename sbs_type::rp& ranked = sbs.m_ranked_points[i];
+            typename sbs_type::rp const& ranked = sbs.m_ranked_points[i];
             turn_type& turn = turns[ranked.turn_index];
             turn_operation_type& op = turn.operations[ranked.operation_index];
 
diff --git a/include/boost/geometry/algorithms/detail/overlay/intersection_insert.hpp b/include/boost/geometry/algorithms/detail/overlay/intersection_insert.hpp
index afc1adb9d6..568cc46fdd 100644
--- a/include/boost/geometry/algorithms/detail/overlay/intersection_insert.hpp
+++ b/include/boost/geometry/algorithms/detail/overlay/intersection_insert.hpp
@@ -525,73 +525,6 @@ struct intersection_areal_linear_point
 };
 
 
-struct tupled_output_tag {};
-
-
-template
-<
-    typename GeometryOut,
-    bool IsTupled = geometry::detail::is_tupled_range_values<GeometryOut>::value
->
-struct tag
-    : geometry::tag<GeometryOut>
-{};
-
-template <typename GeometryOut>
-struct tag<GeometryOut, true>
-{
-    typedef tupled_output_tag type;
-};
-
-
-template <typename Geometry1, typename Geometry2, typename TupledOut>
-struct expect_output_p
-{
-    static const bool is_point_found = geometry::tuples::exists_if
-        <
-            TupledOut, geometry::detail::is_tag_same_as_pred<point_tag>::template pred
-        >::value;
-
-    BOOST_MPL_ASSERT_MSG
-        (
-            is_point_found, POINTLIKE_GEOMETRY_EXPECTED_IN_TUPLED_OUTPUT,
-            (types<Geometry1, Geometry2, TupledOut>)
-        );
-};
-
-template <typename Geometry1, typename Geometry2, typename TupledOut>
-struct expect_output_pl
-    : expect_output_p<Geometry1, Geometry2, TupledOut>
-{
-    static const bool is_linestring_found = geometry::tuples::exists_if
-        <
-            TupledOut, geometry::detail::is_tag_same_as_pred<linestring_tag>::template pred
-        >::value;
-
-    BOOST_MPL_ASSERT_MSG
-        (
-            is_linestring_found, LINEAR_GEOMETRY_EXPECTED_IN_TUPLED_OUTPUT,
-            (types<Geometry1, Geometry2, TupledOut>)
-        );
-};
-
-template <typename Geometry1, typename Geometry2, typename TupledOut>
-struct expect_output_pla
-    : expect_output_pl<Geometry1, Geometry2, TupledOut>
-{
-    static const bool is_polygon_found = geometry::tuples::exists_if
-        <
-            TupledOut, geometry::detail::is_tag_same_as_pred<polygon_tag>::template pred
-        >::value;
-
-    BOOST_MPL_ASSERT_MSG
-    (
-        is_polygon_found, AREAL_GEOMETRY_EXPECTED_IN_TUPLED_OUTPUT,
-        (types<Geometry1, Geometry2, TupledOut>)
-    );
-};
-
-
 }} // namespace detail::intersection
 #endif // DOXYGEN_NO_DETAIL
 
@@ -614,7 +547,7 @@ template
     // tag dispatching:
     typename TagIn1 = typename geometry::tag<Geometry1>::type,
     typename TagIn2 = typename geometry::tag<Geometry2>::type,
-    typename TagOut = typename detail::intersection::tag<GeometryOut>::type,
+    typename TagOut = typename detail::setop_insert_output_tag<GeometryOut>::type,
     // metafunction finetuning helpers:
     typename CastedTagIn1 = typename geometry::tag_cast<TagIn1, areal_tag, linear_tag, pointlike_tag>::type,
     typename CastedTagIn2 = typename geometry::tag_cast<TagIn2, areal_tag, linear_tag, pointlike_tag>::type,
@@ -996,13 +929,21 @@ struct intersection_insert
     <
         Linear1, Linear2, TupledOut, OverlayType,
         Reverse1, Reverse2,
-        TagIn1, TagIn2, detail::intersection::tupled_output_tag,
-        linear_tag, linear_tag, detail::intersection::tupled_output_tag
+        TagIn1, TagIn2, detail::tupled_output_tag,
+        linear_tag, linear_tag, detail::tupled_output_tag
     >
-    // NOTE: This is not fully correct because points can be the result only in
-    // case of intersection but intersection_insert is called also by difference.
-    // So this requirement could be relaxed in the future.
-    : detail::intersection::expect_output_pl<Linear1, Linear2, TupledOut>
+    : detail::expect_output
+        <
+            Linear1, Linear2, TupledOut,
+            // NOTE: points can be the result only in case of intersection.
+            typename boost::mpl::if_c
+                <
+                    (OverlayType == overlay_intersection),
+                    point_tag,
+                    void
+                >::type,
+            linestring_tag
+        >
 {
     // NOTE: The order of geometries in TupledOut tuple/pair must correspond to the order
     // iterators in OutputIterators tuple/pair.
@@ -1113,10 +1054,10 @@ struct intersection_insert
     <
         PointLike1, PointLike2, TupledOut, OverlayType,
         Reverse1, Reverse2,
-        TagIn1, TagIn2, detail::intersection::tupled_output_tag,
-        pointlike_tag, pointlike_tag, detail::intersection::tupled_output_tag
+        TagIn1, TagIn2, detail::tupled_output_tag,
+        pointlike_tag, pointlike_tag, detail::tupled_output_tag
     >
-    : detail::intersection::expect_output_p<PointLike1, PointLike2, TupledOut>
+    : detail::expect_output<PointLike1, PointLike2, TupledOut, point_tag>
 {
     // NOTE: The order of geometries in TupledOut tuple/pair must correspond to the order
     // of iterators in OutputIterators tuple/pair.
@@ -1240,16 +1181,16 @@ struct intersection_insert
     <
         PointLike, Linear, TupledOut, OverlayType,
         Reverse1, Reverse2,
-        TagIn1, TagIn2, detail::intersection::tupled_output_tag,
-        pointlike_tag, linear_tag, detail::intersection::tupled_output_tag
+        TagIn1, TagIn2, detail::tupled_output_tag,
+        pointlike_tag, linear_tag, detail::tupled_output_tag
     >
     // Reuse the implementation for PointLike/PointLike.
     : intersection_insert
         <
             PointLike, Linear, TupledOut, OverlayType,
             Reverse1, Reverse2,
-            TagIn1, TagIn2, detail::intersection::tupled_output_tag,
-            pointlike_tag, pointlike_tag, detail::intersection::tupled_output_tag
+            TagIn1, TagIn2, detail::tupled_output_tag,
+            pointlike_tag, pointlike_tag, detail::tupled_output_tag
         >
 {};
 
@@ -1265,8 +1206,8 @@ struct intersection_insert
     <
         Linestring, MultiPoint, TupledOut, overlay_intersection,
         Reverse1, Reverse2,
-        linestring_tag, multi_point_tag, detail::intersection::tupled_output_tag,
-        linear_tag, pointlike_tag, detail::intersection::tupled_output_tag
+        linestring_tag, multi_point_tag, detail::tupled_output_tag,
+        linear_tag, pointlike_tag, detail::tupled_output_tag
     >
 {
     template <typename RobustPolicy, typename OutputIterators, typename Strategy>
@@ -1368,16 +1309,16 @@ struct intersection_insert
     <
         PointLike, Areal, TupledOut, OverlayType,
         Reverse1, Reverse2,
-        TagIn1, TagIn2, detail::intersection::tupled_output_tag,
-        pointlike_tag, areal_tag, detail::intersection::tupled_output_tag
+        TagIn1, TagIn2, detail::tupled_output_tag,
+        pointlike_tag, areal_tag, detail::tupled_output_tag
     >
     // Reuse the implementation for PointLike/PointLike.
     : intersection_insert
         <
             PointLike, Areal, TupledOut, OverlayType,
             Reverse1, Reverse2,
-            TagIn1, TagIn2, detail::intersection::tupled_output_tag,
-            pointlike_tag, pointlike_tag, detail::intersection::tupled_output_tag
+            TagIn1, TagIn2, detail::tupled_output_tag,
+            pointlike_tag, pointlike_tag, detail::tupled_output_tag
         >
 {};
 
@@ -1394,8 +1335,8 @@ struct intersection_insert
     <
         Areal, MultiPoint, TupledOut, overlay_intersection,
         Reverse1, Reverse2,
-        TagIn1, multi_point_tag, detail::intersection::tupled_output_tag,
-        areal_tag, pointlike_tag, detail::intersection::tupled_output_tag
+        TagIn1, multi_point_tag, detail::tupled_output_tag,
+        areal_tag, pointlike_tag, detail::tupled_output_tag
     >
 {
     template <typename RobustPolicy, typename OutputIterators, typename Strategy>
@@ -1426,8 +1367,8 @@ struct intersection_insert
         TupledOut,
         OverlayType,
         ReverseLinestring, ReversePolygon,
-        linestring_tag, polygon_tag, detail::intersection::tupled_output_tag,
-        linear_tag, areal_tag, detail::intersection::tupled_output_tag
+        linestring_tag, polygon_tag, detail::tupled_output_tag,
+        linear_tag, areal_tag, detail::tupled_output_tag
     > : detail::intersection::intersection_of_linestring_with_areal
             <
                 ReversePolygon,
@@ -1450,8 +1391,8 @@ struct intersection_insert
         TupledOut,
         OverlayType,
         ReverseLinestring, ReverseRing,
-        linestring_tag, ring_tag, detail::intersection::tupled_output_tag,
-        linear_tag, areal_tag, detail::intersection::tupled_output_tag
+        linestring_tag, ring_tag, detail::tupled_output_tag,
+        linear_tag, areal_tag, detail::tupled_output_tag
     > : detail::intersection::intersection_of_linestring_with_areal
             <
                 ReverseRing,
diff --git a/include/boost/geometry/algorithms/detail/overlay/sort_by_side.hpp b/include/boost/geometry/algorithms/detail/overlay/sort_by_side.hpp
index 903055813f..484a439e0a 100644
--- a/include/boost/geometry/algorithms/detail/overlay/sort_by_side.hpp
+++ b/include/boost/geometry/algorithms/detail/overlay/sort_by_side.hpp
@@ -186,7 +186,14 @@ struct less_by_side
             return on_same(first, second);
         }
 
-        int const side_first_wrt_second = -side_second_wrt_first;
+        int const side_first_wrt_second = m_strategy.apply(m_turn_point, second.point, first.point);
+        if (side_second_wrt_first != -side_first_wrt_second)
+        {
+            // (FP) accuracy error in side calculation, the sides are not opposite.
+            // In that case they can be handled as collinear.
+            // If not, then the sort-order might not be stable.
+            return on_same(first, second);
+        }
 
         // Both are on same side, and not collinear
         // Union: return true if second is right w.r.t. first, so -1,
@@ -367,9 +374,11 @@ public :
         }
     }
 
-    template <signed_size_type segment_identifier::*Member, typename Map>
-    void find_open_generic(Map& handled, bool check)
+    void find_open_by_piece_index()
     {
+        // For buffers, use piece index
+        std::set<signed_size_type> handled;
+
         for (std::size_t i = 0; i < m_ranked_points.size(); i++)
         {
             const rp& ranked = m_ranked_points[i];
@@ -378,17 +387,35 @@ public :
                 continue;
             }
 
-            signed_size_type const& index = ranked.seg_id.*Member;
-            if (check && (index < 0 || index > 1))
+            signed_size_type const& index = ranked.seg_id.piece_index;
+            if (handled.count(index) > 0)
             {
-                // Should not occur
                 continue;
             }
-            if (! handled[index])
+            find_polygons_for_source<&segment_identifier::piece_index>(index, i);
+            handled.insert(index);
+        }
+    }
+
+    void find_open_by_source_index()
+    {
+        // Check for source index 0 and 1
+        bool handled[2] = {false, false};
+        for (std::size_t i = 0; i < m_ranked_points.size(); i++)
+        {
+            const rp& ranked = m_ranked_points[i];
+            if (ranked.direction != dir_from)
+            {
+                continue;
+            }
+
+            signed_size_type const& index = ranked.seg_id.source_index;
+            if (index < 0 || index > 1 || handled[index])
             {
-                find_polygons_for_source<Member>(index, i);
-                handled[index] = true;
+                continue;
             }
+            find_polygons_for_source<&segment_identifier::source_index>(index, i);
+            handled[index] = true;
         }
     }
 
@@ -396,21 +423,11 @@ public :
     {
         if (BOOST_GEOMETRY_CONDITION(OverlayType == overlay_buffer))
         {
-            // For buffers, use piece index
-            std::map<signed_size_type, bool> handled;
-            find_open_generic
-                <
-                    &segment_identifier::piece_index
-                >(handled, false);
+            find_open_by_piece_index();
         }
         else
         {
-            // For other operations, by source (there should only source 0,1)
-            bool handled[2] = {false, false};
-            find_open_generic
-                <
-                    &segment_identifier::source_index
-                >(handled, true);
+            find_open_by_source_index();
         }
     }
 
diff --git a/include/boost/geometry/algorithms/detail/tupled_output.hpp b/include/boost/geometry/algorithms/detail/tupled_output.hpp
index 4301716939..b66ce54797 100644
--- a/include/boost/geometry/algorithms/detail/tupled_output.hpp
+++ b/include/boost/geometry/algorithms/detail/tupled_output.hpp
@@ -9,13 +9,17 @@
 #ifndef BOOST_GEOMETRY_ALGORITHMS_DETAIL_TUPLED_OUTPUT_HPP
 #define BOOST_GEOMETRY_ALGORITHMS_DETAIL_TUPLED_OUTPUT_HPP
 
+#include <boost/geometry/algorithms/convert.hpp>
 #include <boost/geometry/core/config.hpp>
 #include <boost/geometry/core/tag.hpp>
+#include <boost/geometry/core/tag_cast.hpp>
 #include <boost/geometry/core/tags.hpp>
+#include <boost/geometry/geometries/concepts/check.hpp>
 #include <boost/geometry/util/range.hpp>
 #include <boost/geometry/util/tuples.hpp>
 
 #include <boost/mpl/and.hpp>
+#include <boost/mpl/if.hpp>
 #include <boost/range/value_type.hpp>
 #include <boost/type_traits/detail/yes_no_type.hpp>
 #include <boost/type_traits/integral_constant.hpp>
@@ -53,7 +57,7 @@ struct is_multi_geometry
 
 // true for point, linestring or polygon
 template <typename T>
-struct is_multi_geometry_value
+struct is_multi_geometry_element
     : boost::integral_constant
         <
             bool,
@@ -86,23 +90,6 @@ struct is_range
 {};
 
 
-// geometry tag of Rng value_type
-template <typename Rng>
-struct range_value_tag
-    : geometry::tag<typename boost::range_value<Rng>::type>
-{};
-
-// true if geometry tag of Rng is the same as Tag
-template <typename Rng, typename Tag>
-struct is_range_value_tag_same_as
-    : boost::is_same
-        <
-            typename range_value_tag<Rng>::type,
-            Tag
-        >
-{};
-
-
 template <typename T, bool IsRange = is_range<T>::value>
 struct is_tupled_output_element_base
     : boost::integral_constant<bool, false>
@@ -114,12 +101,13 @@ struct is_tupled_output_element_base<T, true>
         <
             bool,
             (is_multi_geometry<T>::value
-            ||
-            ((! is_geometry<T>::value)
-                &&
-                ((is_range_value_tag_same_as<T, point_tag>::value)
-                || (is_range_value_tag_same_as<T, linestring_tag>::value)
-                || (is_range_value_tag_same_as<T, polygon_tag>::value))))
+                ||
+                ((! is_geometry<T>::value)
+                    &&
+                    is_multi_geometry_element
+                        <
+                            typename boost::range_value<T>::type
+                        >::value))
         >
 {};
 
@@ -134,7 +122,7 @@ struct is_tupled_output_element
 
 // true if Output is not a geometry (so e.g. tuple was not adapted to any
 // concept) and at least one of the tuple elements is a multi-geometry or
-// a range of points linestrings or polygons
+// a range of points, linestrings or polygons
 template <typename Output>
 struct is_tupled_output_check
     : boost::mpl::and_
@@ -146,28 +134,15 @@ struct is_tupled_output_check
 {};
 
 
-// true if T is a point, linestring or polygon
-template <typename T>
-struct is_tupled_range_values_element
-    : boost::integral_constant
-        <
-            bool,
-            ((boost::is_same<typename geometry::tag<T>::type, point_tag>::value)
-                || (boost::is_same<typename geometry::tag<T>::type, linestring_tag>::value)
-                || (boost::is_same<typename geometry::tag<T>::type, polygon_tag>::value))
-        >
-{};
-
-
 // true if T is not a geometry (so e.g. tuple was not adapted to any
 // concept) and at least one of the tuple elements is a point, linesting
 // or polygon
 template <typename T>
-struct is_tupled_range_values_check
+struct is_tupled_single_output_check
     : boost::mpl::and_
         <
             boost::is_same<typename geometry::tag<T>::type, void>,
-            geometry::tuples::exists_if<T, is_tupled_range_values_element>
+            geometry::tuples::exists_if<T, is_multi_geometry_element>
         >
 {};
 
@@ -222,15 +197,15 @@ struct is_tupled_output<Output, true>
 
 
 // true if T is boost::tuple, boost::tuples::cons, std::pair or std::tuple
-// and is_tupled_range_values_check defiend above passes
+// and is_tupled_single_output_check defiend above passes
 template <typename T, bool IsTupled = is_tupled<T>::value>
-struct is_tupled_range_values
+struct is_tupled_single_output
     : boost::integral_constant<bool, false>
 {};
 
 template <typename T>
-struct is_tupled_range_values<T, true>
-    : is_tupled_range_values_check<T>
+struct is_tupled_single_output<T, true>
+    : is_tupled_single_output_check<T>
 {};
 
 
@@ -527,6 +502,234 @@ struct output_geometry_access<GeometryOut, Tag, DefaultTag, DefaultTag>
 };
 
 
+template <typename Geometry>
+struct output_geometry_concept_check
+{
+    static void apply()
+    {
+        concepts::check<Geometry>();
+    }
+};
+
+template <typename First, typename Second>
+struct output_geometry_concept_check<std::pair<First, Second> >
+{
+    static void apply()
+    {
+        concepts::check<First>();
+        concepts::check<Second>();
+    }
+};
+
+template <typename Tuple,
+          size_t I = 0,
+          size_t N = geometry::tuples::size<Tuple>::value>
+struct output_geometry_concept_check_t
+{
+    static void apply()
+    {
+        concepts::check<typename geometry::tuples::element<I, Tuple>::type>();
+        output_geometry_concept_check_t<Tuple, I + 1, N>::apply();
+    }
+};
+
+template <typename Tuple, size_t N>
+struct output_geometry_concept_check_t<Tuple, N, N>
+{
+    static void apply()
+    {}
+};
+
+template
+<
+    class T0, class T1, class T2, class T3, class T4,
+    class T5, class T6, class T7, class T8, class T9
+>
+struct output_geometry_concept_check<boost::tuple<T0, T1, T2, T3, T4, T5, T6, T7, T8, T9> >
+    : output_geometry_concept_check_t<boost::tuple<T0, T1, T2, T3, T4, T5, T6, T7, T8, T9> >
+{};
+
+template <typename HT, typename TT>
+struct output_geometry_concept_check<boost::tuples::cons<HT, TT> >
+    : output_geometry_concept_check_t<boost::tuples::cons<HT, TT> >
+{};
+
+#ifdef BOOST_GEOMETRY_CXX11_TUPLE
+
+template <typename ...Ts>
+struct output_geometry_concept_check<std::tuple<Ts...> >
+    : output_geometry_concept_check_t<std::tuple<Ts...> >
+{};
+
+#endif // BOOST_GEOMETRY_CXX11_TUPLE
+
+
+struct tupled_output_tag {};
+
+
+template <typename GeometryOut>
+struct setop_insert_output_tag
+    : boost::mpl::if_c
+        <
+            geometry::detail::is_tupled_single_output<GeometryOut>::value,
+            tupled_output_tag,
+            typename geometry::tag<GeometryOut>::type
+        >
+{};
+
+
+template <typename Geometry1, typename Geometry2, typename TupledOut, bool IsFound, typename Tag>
+struct expect_output_assert_base;
+
+template <typename Geometry1, typename Geometry2, typename TupledOut, bool IsFound>
+struct expect_output_assert_base<Geometry1, Geometry2, TupledOut, IsFound, pointlike_tag>
+{
+    BOOST_MPL_ASSERT_MSG
+        (
+            IsFound, POINTLIKE_GEOMETRY_EXPECTED_IN_TUPLED_OUTPUT,
+            (types<Geometry1, Geometry2, TupledOut>)
+        );
+};
+
+template <typename Geometry1, typename Geometry2, typename TupledOut, bool IsFound>
+struct expect_output_assert_base<Geometry1, Geometry2, TupledOut, IsFound, linear_tag>
+{
+    BOOST_MPL_ASSERT_MSG
+    (
+        IsFound, LINEAR_GEOMETRY_EXPECTED_IN_TUPLED_OUTPUT,
+        (types<Geometry1, Geometry2, TupledOut>)
+    );
+};
+
+template <typename Geometry1, typename Geometry2, typename TupledOut, bool IsFound>
+struct expect_output_assert_base<Geometry1, Geometry2, TupledOut, IsFound, areal_tag>
+{
+    BOOST_MPL_ASSERT_MSG
+    (
+        IsFound, AREAL_GEOMETRY_EXPECTED_IN_TUPLED_OUTPUT,
+        (types<Geometry1, Geometry2, TupledOut>)
+    );
+};
+
+
+template <typename Geometry1, typename Geometry2, typename TupledOut, typename Tag>
+struct expect_output_assert
+    : expect_output_assert_base
+        <
+            Geometry1, Geometry2, TupledOut,
+            geometry::tuples::exists_if
+                <
+                    TupledOut,
+                    is_tag_same_as_pred<Tag>::template pred
+                >::value,
+            typename geometry::tag_cast
+                <
+                    Tag, pointlike_tag, linear_tag, areal_tag
+                >::type
+        >
+{};
+
+template <typename Geometry1, typename Geometry2, typename TupledOut>
+struct expect_output_assert<Geometry1, Geometry2, TupledOut, void>
+{};
+
+template
+<
+    typename Geometry1, typename Geometry2, typename TupledOut,
+    typename Tag1,
+    typename Tag2 = void,
+    typename Tag3 = void
+>
+struct expect_output
+    : expect_output_assert<Geometry1, Geometry2, TupledOut, Tag1>
+    , expect_output_assert<Geometry1, Geometry2, TupledOut, Tag2>
+    , expect_output_assert<Geometry1, Geometry2, TupledOut, Tag3>
+{};
+
+template
+<
+    typename Geometry1, typename Geometry2, typename TupledOut,
+    typename Tag1, typename Tag2
+>
+struct expect_output<Geometry1, Geometry2, TupledOut, Tag1, Tag2, void>
+    : expect_output_assert<Geometry1, Geometry2, TupledOut, Tag1>
+    , expect_output_assert<Geometry1, Geometry2, TupledOut, Tag2>
+{};
+
+template
+<
+    typename Geometry1, typename Geometry2, typename TupledOut,
+    typename Tag1
+>
+struct expect_output<Geometry1, Geometry2, TupledOut, Tag1, void, void>
+    : expect_output_assert<Geometry1, Geometry2, TupledOut, Tag1>
+{};
+
+
+template <typename CastedTag>
+struct single_tag_from_base_tag;
+
+template <>
+struct single_tag_from_base_tag<pointlike_tag>
+{
+    typedef point_tag type;
+};
+
+template <>
+struct single_tag_from_base_tag<linear_tag>
+{
+    typedef linestring_tag type;
+};
+
+template <>
+struct single_tag_from_base_tag<areal_tag>
+{
+    typedef polygon_tag type;
+};
+
+
+template
+<
+    typename Geometry,
+    typename SingleOut,
+    bool IsMulti = geometry::detail::is_multi_geometry<Geometry>::value
+>
+struct convert_to_output
+{
+    template <typename OutputIterator>
+    static OutputIterator apply(Geometry const& geometry,
+                                OutputIterator oit)
+    {
+        SingleOut single_out;
+        geometry::convert(geometry, single_out);
+        *oit++ = single_out;
+        return oit;
+    }
+};
+
+template
+<
+    typename Geometry,
+    typename SingleOut
+>
+struct convert_to_output<Geometry, SingleOut, true>
+{
+    template <typename OutputIterator>
+    static OutputIterator apply(Geometry const& geometry,
+                                OutputIterator oit)
+    {
+        typedef typename boost::range_iterator<Geometry const>::type iterator;
+        for (iterator it = boost::begin(geometry); it != boost::end(geometry); ++it)
+        {
+            SingleOut single_out;
+            geometry::convert(*it, single_out);
+            *oit++ = single_out;
+        }
+        return oit;
+    }
+};
+
+
 } // namespace detail
 #endif // DOXYGEN_NO_DETAIL
 
diff --git a/include/boost/geometry/algorithms/difference.hpp b/include/boost/geometry/algorithms/difference.hpp
index 7a96e130dc..be96eea5cb 100644
--- a/include/boost/geometry/algorithms/difference.hpp
+++ b/include/boost/geometry/algorithms/difference.hpp
@@ -2,8 +2,8 @@
 
 // Copyright (c) 2007-2012 Barend Gehrels, Amsterdam, the Netherlands.
 
-// This file was modified by Oracle on 2017, 2019.
-// Modifications copyright (c) 2017, 2019, Oracle and/or its affiliates.
+// This file was modified by Oracle on 2017, 2019, 2020.
+// Modifications copyright (c) 2017-2020, Oracle and/or its affiliates.
 
 // Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle
 
@@ -19,6 +19,7 @@
 #include <boost/variant/static_visitor.hpp>
 #include <boost/variant/variant_fwd.hpp>
 
+#include <boost/geometry/algorithms/detail/intersection/multi.hpp>
 #include <boost/geometry/algorithms/detail/overlay/intersection_insert.hpp>
 #include <boost/geometry/policies/robustness/get_rescale_policy.hpp>
 #include <boost/geometry/strategies/default_strategy.hpp>
@@ -32,6 +33,146 @@ namespace boost { namespace geometry
 namespace detail { namespace difference
 {
 
+template
+<
+    typename Geometry1,
+    typename Geometry2,
+    typename SingleOut,
+    typename OutTag = typename detail::setop_insert_output_tag<SingleOut>::type,
+    bool ReturnGeometry1 = (topological_dimension<Geometry1>::value
+                            > topological_dimension<Geometry2>::value)
+>
+struct call_intersection_insert
+{
+    template
+    <
+        typename OutputIterator,
+        typename RobustPolicy,
+        typename Strategy
+    >
+    static inline OutputIterator apply(Geometry1 const& geometry1,
+                                       Geometry2 const& geometry2,
+                                       RobustPolicy const& robust_policy,
+                                       OutputIterator out,
+                                       Strategy const& strategy)
+    {
+        return geometry::dispatch::intersection_insert
+            <
+                Geometry1, Geometry2,
+                SingleOut,
+                overlay_difference,
+                geometry::detail::overlay::do_reverse<geometry::point_order<Geometry1>::value>::value,
+                geometry::detail::overlay::do_reverse<geometry::point_order<Geometry2>::value, true>::value
+            >::apply(geometry1, geometry2, robust_policy, out, strategy);
+    }
+};
+
+template
+<
+    typename Geometry1,
+    typename Geometry2,
+    typename SingleOut
+>
+struct call_intersection_insert_tupled_base
+{
+    typedef typename geometry::detail::single_tag_from_base_tag
+        <
+            typename geometry::tag_cast
+                <
+                    typename geometry::tag<Geometry1>::type,
+                    pointlike_tag, linear_tag, areal_tag
+                >::type
+        >::type single_tag;
+
+    typedef detail::expect_output
+        <
+            Geometry1, Geometry2, SingleOut, single_tag
+        > expect_check;
+
+    typedef typename geometry::detail::output_geometry_access
+        <
+            SingleOut, single_tag, single_tag
+        > access;
+};
+
+template
+<
+    typename Geometry1,
+    typename Geometry2,
+    typename SingleOut
+>
+struct call_intersection_insert
+    <
+        Geometry1, Geometry2, SingleOut,
+        detail::tupled_output_tag,
+        false
+    >
+    : call_intersection_insert_tupled_base<Geometry1, Geometry2, SingleOut>
+{
+    typedef call_intersection_insert_tupled_base<Geometry1, Geometry2, SingleOut> base_t;
+
+    template
+    <
+        typename OutputIterator,
+        typename RobustPolicy,
+        typename Strategy
+    >
+    static inline OutputIterator apply(Geometry1 const& geometry1,
+                                       Geometry2 const& geometry2,
+                                       RobustPolicy const& robust_policy,
+                                       OutputIterator out,
+                                       Strategy const& strategy)
+    {
+        base_t::access::get(out) = call_intersection_insert
+            <
+                Geometry1, Geometry2,
+                typename base_t::access::type
+            >::apply(geometry1, geometry2, robust_policy,
+                     base_t::access::get(out), strategy);
+
+        return out;
+    }
+};
+
+template
+<
+    typename Geometry1,
+    typename Geometry2,
+    typename SingleOut
+>
+struct call_intersection_insert
+    <
+        Geometry1, Geometry2, SingleOut,
+        detail::tupled_output_tag,
+        true
+    >
+    : call_intersection_insert_tupled_base<Geometry1, Geometry2, SingleOut>
+{
+    typedef call_intersection_insert_tupled_base<Geometry1, Geometry2, SingleOut> base_t;
+
+    template
+    <
+        typename OutputIterator,
+        typename RobustPolicy,
+        typename Strategy
+    >
+    static inline OutputIterator apply(Geometry1 const& geometry1,
+                                       Geometry2 const& geometry2,
+                                       RobustPolicy const& robust_policy,
+                                       OutputIterator out,
+                                       Strategy const& strategy)
+    {
+        base_t::access::get(out) = geometry::detail::convert_to_output
+            <
+                Geometry1,
+                typename base_t::access::type
+            >::apply(geometry1, base_t::access::get(out));
+
+        return out;
+    }
+};
+
+
 /*!
 \brief_calc2{difference} \brief_strategy
 \ingroup difference
@@ -65,7 +206,8 @@ inline OutputIterator difference_insert(Geometry1 const& geometry1,
 {
     concepts::check<Geometry1 const>();
     concepts::check<Geometry2 const>();
-    concepts::check<GeometryOut>();
+    //concepts::check<GeometryOut>();
+    geometry::detail::output_geometry_concept_check<GeometryOut>::apply();
 
     typedef typename geometry::rescale_overlay_policy_type
         <
@@ -75,16 +217,12 @@ inline OutputIterator difference_insert(Geometry1 const& geometry1,
         >::type rescale_policy_type;
 
     rescale_policy_type robust_policy
-            = geometry::get_rescale_policy<rescale_policy_type>(
-                geometry1, geometry2, strategy);
+        = geometry::get_rescale_policy<rescale_policy_type>(
+            geometry1, geometry2, strategy);
 
-    return geometry::dispatch::intersection_insert
+    return geometry::detail::difference::call_intersection_insert
         <
-            Geometry1, Geometry2,
-            GeometryOut,
-            overlay_difference,
-            geometry::detail::overlay::do_reverse<geometry::point_order<Geometry1>::value>::value,
-            geometry::detail::overlay::do_reverse<geometry::point_order<Geometry2>::value, true>::value
+            Geometry1, Geometry2, GeometryOut
         >::apply(geometry1, geometry2, robust_policy, out, strategy);
 }
 
@@ -146,11 +284,14 @@ struct difference
                              Collection & output_collection,
                              Strategy const& strategy)
     {
-        typedef typename boost::range_value<Collection>::type geometry_out;
+        typedef typename geometry::detail::output_geometry_value
+            <
+                Collection
+            >::type single_out;
 
-        detail::difference::difference_insert<geometry_out>(
+        detail::difference::difference_insert<single_out>(
             geometry1, geometry2,
-            range::back_inserter(output_collection),
+            geometry::detail::output_geometry_back_inserter(output_collection),
             strategy);
     }
 
@@ -165,11 +306,13 @@ struct difference
                              Collection & output_collection,
                              default_strategy)
     {
-        typedef typename boost::range_value<Collection>::type geometry_out;
+        typedef typename strategy::relate::services::default_strategy
+            <
+                Geometry1,
+                Geometry2
+            >::type strategy_type;
         
-        detail::difference::difference_insert<geometry_out>(
-            geometry1, geometry2,
-            range::back_inserter(output_collection));
+        apply(geometry1, geometry2, output_collection, strategy_type());
     }
 };
 
diff --git a/include/boost/geometry/algorithms/sym_difference.hpp b/include/boost/geometry/algorithms/sym_difference.hpp
index a281bf1db5..4d8550bcb1 100644
--- a/include/boost/geometry/algorithms/sym_difference.hpp
+++ b/include/boost/geometry/algorithms/sym_difference.hpp
@@ -2,8 +2,8 @@
 
 // Copyright (c) 2007-2015 Barend Gehrels, Amsterdam, the Netherlands.
 
-// This file was modified by Oracle on 2015, 2017, 2019.
-// Modifications copyright (c) 2015-2019 Oracle and/or its affiliates.
+// This file was modified by Oracle on 2015, 2017, 2019, 2020.
+// Modifications copyright (c) 2015-2020 Oracle and/or its affiliates.
 
 // Contributed and/or modified by Menelaos Karavelas, on behalf of Oracle
 // Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle
@@ -29,6 +29,7 @@
 #include <boost/geometry/geometries/multi_polygon.hpp>
 #include <boost/geometry/policies/robustness/get_rescale_policy.hpp>
 #include <boost/geometry/strategies/default_strategy.hpp>
+#include <boost/geometry/strategies/relate.hpp>
 #include <boost/geometry/util/range.hpp>
 
 
@@ -150,6 +151,124 @@ struct sym_difference_areal_areal
 };
 
 
+template
+<
+    typename GeometryOut,
+    typename SingleTag,
+    template <typename, typename, typename> class Algorithm
+>
+struct sym_difference_same_inputs_tupled_output
+{
+    template
+    <
+        typename Geometry1,
+        typename Geometry2,
+        typename RobustPolicy,
+        typename OutputIterator,
+        typename Strategy
+    >
+    static inline OutputIterator apply(Geometry1 const& geometry1,
+                                       Geometry2 const& geometry2,
+                                       RobustPolicy const& robust_policy,
+                                       OutputIterator out,
+                                       Strategy const& strategy)
+    {
+        typedef typename geometry::detail::output_geometry_access
+            <
+                GeometryOut, SingleTag, SingleTag
+            > access;
+
+        access::get(out) = Algorithm
+            <
+                typename access::type, Geometry1, Geometry2
+            >::apply(geometry1, geometry2, robust_policy, access::get(out), strategy);
+
+        return out;
+    }
+};
+
+
+template
+<
+    typename GeometryOut,
+    typename SingleTag1,
+    typename SingleTag2,
+    bool Reverse = (geometry::core_dispatch::top_dim<SingleTag1>::value
+                    > geometry::core_dispatch::top_dim<SingleTag2>::value)
+>
+struct sym_difference_different_inputs_tupled_output
+{
+    template
+    <
+        typename Geometry1,
+        typename Geometry2,
+        typename RobustPolicy,
+        typename OutputIterator,
+        typename Strategy
+    >
+    static inline OutputIterator apply(Geometry1 const& geometry1,
+                                       Geometry2 const& geometry2,
+                                       RobustPolicy const& robust_policy,
+                                       OutputIterator out,
+                                       Strategy const& strategy)
+    {
+        return sym_difference_different_inputs_tupled_output
+            <
+                GeometryOut, SingleTag2, SingleTag1
+            >::apply(geometry2, geometry1, robust_policy, out, strategy);
+    }
+};
+
+template
+<
+    typename GeometryOut,
+    typename SingleTag1,
+    typename SingleTag2
+>
+struct sym_difference_different_inputs_tupled_output
+    <
+        GeometryOut, SingleTag1, SingleTag2, false
+    >
+{
+    template
+    <
+        typename Geometry1,
+        typename Geometry2,
+        typename RobustPolicy,
+        typename OutputIterator,
+        typename Strategy
+    >
+    static inline OutputIterator apply(Geometry1 const& geometry1,
+                                       Geometry2 const& geometry2,
+                                       RobustPolicy const& robust_policy,
+                                       OutputIterator out,
+                                       Strategy const& strategy)
+    {
+        typedef typename geometry::detail::output_geometry_access
+            <
+                GeometryOut, SingleTag1, SingleTag1
+            > access1;
+        typedef typename geometry::detail::output_geometry_access
+            <
+                GeometryOut, SingleTag2, SingleTag2
+            > access2;
+
+        access1::get(out) = compute_difference
+            <
+                typename access1::type
+            >::apply(geometry1, geometry2, robust_policy, access1::get(out), strategy);
+
+        access2::get(out) = geometry::detail::convert_to_output
+            <
+                Geometry2,
+                typename access2::type
+            >::apply(geometry2, access2::get(out));
+
+        return out;
+    }
+};
+
+
 }} // namespace detail::sym_difference
 #endif // DOXYGEN_NO_DETAIL
 
@@ -167,13 +286,13 @@ template
     typename GeometryOut,
     typename TagIn1 = typename geometry::tag_cast
         <
-            typename tag<Geometry1>::type, areal_tag
+            typename tag<Geometry1>::type, pointlike_tag, linear_tag, areal_tag
         >::type,
     typename TagIn2 = typename geometry::tag_cast
         <
-            typename tag<Geometry2>::type, areal_tag
+            typename tag<Geometry2>::type, pointlike_tag, linear_tag, areal_tag
         >::type,
-    typename TagOut = typename geometry::tag<GeometryOut>::type
+    typename TagOut = typename detail::setop_insert_output_tag<GeometryOut>::type
 >
 struct sym_difference_insert
     : detail::sym_difference::sym_difference_generic
@@ -201,6 +320,95 @@ struct sym_difference_insert
 {};
 
 
+
+template
+<
+    typename PointLike1,
+    typename PointLike2,
+    typename GeometryOut
+>
+struct sym_difference_insert
+    <
+        PointLike1, PointLike2, GeometryOut,
+        pointlike_tag, pointlike_tag, detail::tupled_output_tag
+    >
+    : detail::expect_output<PointLike1, PointLike2, GeometryOut, point_tag>
+    , detail::sym_difference::sym_difference_same_inputs_tupled_output
+        <
+            GeometryOut,
+            point_tag,
+            detail::sym_difference::sym_difference_generic
+        >
+{};
+
+template
+<
+    typename Linear1,
+    typename Linear2,
+    typename GeometryOut
+>
+struct sym_difference_insert
+    <
+        Linear1, Linear2, GeometryOut,
+        linear_tag, linear_tag, detail::tupled_output_tag
+    >
+    : detail::expect_output<Linear1, Linear2, GeometryOut, linestring_tag>
+    , detail::sym_difference::sym_difference_same_inputs_tupled_output
+        <
+            GeometryOut,
+            linestring_tag,
+            detail::sym_difference::sym_difference_generic
+        >
+{};
+
+template
+<
+    typename Areal1,
+    typename Areal2,
+    typename GeometryOut
+>
+struct sym_difference_insert
+    <
+        Areal1, Areal2, GeometryOut,
+        areal_tag, areal_tag, detail::tupled_output_tag
+    >
+    : detail::expect_output<Areal1, Areal2, GeometryOut, polygon_tag>
+    , detail::sym_difference::sym_difference_same_inputs_tupled_output
+        <
+            GeometryOut,
+            polygon_tag,
+            detail::sym_difference::sym_difference_areal_areal
+        >
+{};
+
+template
+<
+    typename Geometry1,
+    typename Geometry2,
+    typename GeometryOut,
+    typename TagIn1,
+    typename TagIn2
+>
+struct sym_difference_insert
+    <
+        Geometry1, Geometry2, GeometryOut,
+        TagIn1, TagIn2, detail::tupled_output_tag
+    >
+    : detail::expect_output
+        <
+            Geometry1, Geometry2, GeometryOut,
+            typename detail::single_tag_from_base_tag<TagIn1>::type,
+            typename detail::single_tag_from_base_tag<TagIn2>::type
+        >
+    , detail::sym_difference::sym_difference_different_inputs_tupled_output
+        <
+            GeometryOut,
+            typename detail::single_tag_from_base_tag<TagIn1>::type,
+            typename detail::single_tag_from_base_tag<TagIn2>::type
+        >
+{};
+
+
 } // namespace dispatch
 #endif // DOXYGEN_NO_DISPATCH
 
@@ -244,7 +452,8 @@ inline OutputIterator sym_difference_insert(Geometry1 const& geometry1,
 {
     concepts::check<Geometry1 const>();
     concepts::check<Geometry2 const>();
-    concepts::check<GeometryOut>();
+    //concepts::check<GeometryOut>();
+    geometry::detail::output_geometry_concept_check<GeometryOut>::apply();
 
     typedef typename geometry::rescale_overlay_policy_type
         <
@@ -288,10 +497,6 @@ template
 inline OutputIterator sym_difference_insert(Geometry1 const& geometry1,
             Geometry2 const& geometry2, OutputIterator out)
 {
-    concepts::check<Geometry1 const>();
-    concepts::check<Geometry2 const>();
-    concepts::check<GeometryOut>();
-
     typedef typename strategy::intersection::services::default_strategy
         <
             typename cs_tag<GeometryOut>::type
@@ -320,11 +525,14 @@ struct sym_difference
                              Collection & output_collection,
                              Strategy const& strategy)
     {
-        typedef typename boost::range_value<Collection>::type geometry_out;
+        typedef typename geometry::detail::output_geometry_value
+            <
+                Collection
+            >::type single_out;
 
-        detail::sym_difference::sym_difference_insert<geometry_out>(
+        detail::sym_difference::sym_difference_insert<single_out>(
             geometry1, geometry2,
-            range::back_inserter(output_collection),
+            geometry::detail::output_geometry_back_inserter(output_collection),
             strategy);
     }
 
@@ -339,11 +547,12 @@ struct sym_difference
                              Collection & output_collection,
                              default_strategy)
     {
-        typedef typename boost::range_value<Collection>::type geometry_out;
-        
-        detail::sym_difference::sym_difference_insert<geometry_out>(
-            geometry1, geometry2,
-            range::back_inserter(output_collection));
+        typedef typename strategy::relate::services::default_strategy
+            <
+                Geometry1, Geometry2
+            >::type strategy_type;
+
+        apply(geometry1, geometry2, output_collection, strategy_type());
     }
 };
 
diff --git a/include/boost/geometry/algorithms/union.hpp b/include/boost/geometry/algorithms/union.hpp
index de47795629..cfdd1c4aed 100644
--- a/include/boost/geometry/algorithms/union.hpp
+++ b/include/boost/geometry/algorithms/union.hpp
@@ -2,8 +2,8 @@
 
 // Copyright (c) 2007-2014 Barend Gehrels, Amsterdam, the Netherlands.
 
-// This file was modified by Oracle on 2014, 2017, 2018, 2019.
-// Modifications copyright (c) 2014-2019 Oracle and/or its affiliates.
+// This file was modified by Oracle on 2014, 2017, 2018, 2019, 2020.
+// Modifications copyright (c) 2014-2020 Oracle and/or its affiliates.
 
 // Contributed and/or modified by Menelaos Karavelas, on behalf of Oracle
 // Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle
@@ -18,7 +18,6 @@
 
 #include <boost/range/metafunctions.hpp>
 
-#include <boost/geometry/core/is_areal.hpp>
 #include <boost/geometry/core/point_order.hpp>
 #include <boost/geometry/core/reverse_dispatch.hpp>
 #include <boost/geometry/geometries/concepts/check.hpp>
@@ -28,6 +27,8 @@
 #include <boost/geometry/strategies/default_strategy.hpp>
 #include <boost/geometry/util/range.hpp>
 
+#include <boost/geometry/algorithms/detail/intersection/multi.hpp>
+#include <boost/geometry/algorithms/detail/overlay/intersection_insert.hpp>
 #include <boost/geometry/algorithms/detail/overlay/linear_linear.hpp>
 #include <boost/geometry/algorithms/detail/overlay/pointlike_pointlike.hpp>
 
@@ -44,13 +45,10 @@ template
     typename Geometry1, typename Geometry2, typename GeometryOut,
     typename TagIn1 = typename tag<Geometry1>::type,
     typename TagIn2 = typename tag<Geometry2>::type,
-    typename TagOut = typename tag<GeometryOut>::type,
-    bool Areal1 = geometry::is_areal<Geometry1>::value,
-    bool Areal2 = geometry::is_areal<Geometry2>::value,
-    bool ArealOut = geometry::is_areal<GeometryOut>::value,
-    bool Reverse1 = detail::overlay::do_reverse<geometry::point_order<Geometry1>::value>::value,
-    bool Reverse2 = detail::overlay::do_reverse<geometry::point_order<Geometry2>::value>::value,
-    bool ReverseOut = detail::overlay::do_reverse<geometry::point_order<GeometryOut>::value>::value,
+    typename TagOut = typename detail::setop_insert_output_tag<GeometryOut>::type,
+    typename CastedTagIn1 = typename geometry::tag_cast<TagIn1, areal_tag, linear_tag, pointlike_tag>::type,
+    typename CastedTagIn2 = typename geometry::tag_cast<TagIn2, areal_tag, linear_tag, pointlike_tag>::type,
+    typename CastedTagOut = typename geometry::tag_cast<TagOut, areal_tag, linear_tag, pointlike_tag>::type,
     bool Reverse = geometry::reverse_dispatch<Geometry1, Geometry2>::type::value
 >
 struct union_insert: not_implemented<TagIn1, TagIn2, TagOut>
@@ -63,24 +61,22 @@ template
 <
     typename Geometry1, typename Geometry2, typename GeometryOut,
     typename TagIn1, typename TagIn2, typename TagOut,
-    bool Areal1, bool Areal2, bool ArealOut,
-    bool Reverse1, bool Reverse2, bool ReverseOut
+    typename CastedTagIn1, typename CastedTagIn2, typename CastedTagOut
 >
 struct union_insert
     <
         Geometry1, Geometry2, GeometryOut,
         TagIn1, TagIn2, TagOut,
-        Areal1, Areal2, ArealOut,
-        Reverse1, Reverse2, ReverseOut,
+        CastedTagIn1, CastedTagIn2, CastedTagOut,
         true
-    >: union_insert<Geometry2, Geometry1, GeometryOut>
+    >
 {
     template <typename RobustPolicy, typename OutputIterator, typename Strategy>
     static inline OutputIterator apply(Geometry1 const& g1,
-            Geometry2 const& g2,
-            RobustPolicy const& robust_policy,
-            OutputIterator out,
-            Strategy const& strategy)
+                                       Geometry2 const& g2,
+                                       RobustPolicy const& robust_policy,
+                                       OutputIterator out,
+                                       Strategy const& strategy)
     {
         return union_insert
             <
@@ -93,44 +89,22 @@ struct union_insert
 template
 <
     typename Geometry1, typename Geometry2, typename GeometryOut,
-    typename TagIn1, typename TagIn2, typename TagOut,
-    bool Reverse1, bool Reverse2, bool ReverseOut
+    typename TagIn1, typename TagIn2, typename TagOut
 >
 struct union_insert
     <
         Geometry1, Geometry2, GeometryOut,
         TagIn1, TagIn2, TagOut,
-        true, true, true,
-        Reverse1, Reverse2, ReverseOut,
+        areal_tag, areal_tag, areal_tag,
         false
     > : detail::overlay::overlay
-        <Geometry1, Geometry2, Reverse1, Reverse2, ReverseOut, GeometryOut, overlay_union>
-{};
-
-
-// dispatch for union of non-areal geometries
-template
-<
-    typename Geometry1, typename Geometry2, typename GeometryOut,
-    typename TagIn1, typename TagIn2, typename TagOut,
-    bool Reverse1, bool Reverse2, bool ReverseOut
->
-struct union_insert
-    <
-        Geometry1, Geometry2, GeometryOut,
-        TagIn1, TagIn2, TagOut,
-        false, false, false,
-        Reverse1, Reverse2, ReverseOut,
-        false
-    > : union_insert
         <
-            Geometry1, Geometry2, GeometryOut,
-            typename tag_cast<TagIn1, pointlike_tag, linear_tag>::type,
-            typename tag_cast<TagIn2, pointlike_tag, linear_tag>::type,
-            TagOut,
-            false, false, false,
-            Reverse1, Reverse2, ReverseOut,
-            false
+            Geometry1, Geometry2,
+            detail::overlay::do_reverse<geometry::point_order<Geometry1>::value>::value,
+            detail::overlay::do_reverse<geometry::point_order<Geometry2>::value>::value,
+            detail::overlay::do_reverse<geometry::point_order<GeometryOut>::value>::value,
+            GeometryOut,
+            overlay_union
         >
 {};
 
@@ -139,14 +113,13 @@ struct union_insert
 template
 <
     typename Linear1, typename Linear2, typename LineStringOut,
-    bool Reverse1, bool Reverse2, bool ReverseOut
+    typename TagIn1, typename TagIn2
 >
 struct union_insert
     <
         Linear1, Linear2, LineStringOut,
-        linear_tag, linear_tag, linestring_tag,
-        false, false, false,
-        Reverse1, Reverse2, ReverseOut,
+        TagIn1, TagIn2, linestring_tag,
+        linear_tag, linear_tag, linear_tag,
         false
     > : detail::overlay::linear_linear_linestring
         <
@@ -159,14 +132,13 @@ struct union_insert
 template
 <
     typename PointLike1, typename PointLike2, typename PointOut,
-    bool Reverse1, bool Reverse2, bool ReverseOut
+    typename TagIn1, typename TagIn2
 >
 struct union_insert
     <
         PointLike1, PointLike2, PointOut,
-        pointlike_tag, pointlike_tag, point_tag,
-        false, false, false,
-        Reverse1, Reverse2, ReverseOut,
+        TagIn1, TagIn2, point_tag,
+        pointlike_tag, pointlike_tag, pointlike_tag,
         false
     > : detail::overlay::union_pointlike_pointlike_point
         <
@@ -175,6 +147,172 @@ struct union_insert
 {};
 
 
+template
+<
+    typename Geometry1, typename Geometry2, typename SingleTupledOut,
+    typename TagIn1, typename TagIn2,
+    typename CastedTagIn
+>
+struct union_insert
+    <
+        Geometry1, Geometry2, SingleTupledOut,
+        TagIn1, TagIn2, detail::tupled_output_tag,
+        CastedTagIn, CastedTagIn, detail::tupled_output_tag,
+        false
+    >
+{
+    typedef typename geometry::detail::single_tag_from_base_tag
+        <
+            CastedTagIn
+        >::type single_tag;
+
+    typedef detail::expect_output
+        <
+            Geometry1, Geometry2, SingleTupledOut, single_tag
+        > expect_check;
+
+    typedef typename geometry::detail::output_geometry_access
+        <
+            SingleTupledOut, single_tag, single_tag
+        > access;
+
+    template <typename RobustPolicy, typename OutputIterator, typename Strategy>
+    static inline OutputIterator apply(Geometry1 const& g1,
+                                       Geometry2 const& g2,
+                                       RobustPolicy const& robust_policy,
+                                       OutputIterator out,
+                                       Strategy const& strategy)
+    {
+        access::get(out) = union_insert
+            <
+                Geometry2, Geometry1, typename access::type
+            >::apply(g2, g1, robust_policy, access::get(out), strategy);
+
+        return out;
+    }
+};
+
+
+template
+<
+    typename Geometry1, typename Geometry2, typename SingleTupledOut,
+    typename SingleTag1, typename SingleTag2,
+    bool Geometry1LesserTopoDim = (topological_dimension<Geometry1>::value
+                                    < topological_dimension<Geometry2>::value)
+>
+struct union_insert_tupled_different
+{
+    typedef typename geometry::detail::output_geometry_access
+        <
+            SingleTupledOut, SingleTag1, SingleTag1
+        > access1;
+
+    typedef typename geometry::detail::output_geometry_access
+        <
+            SingleTupledOut, SingleTag2, SingleTag2
+        > access2;
+
+    template <typename RobustPolicy, typename OutputIterator, typename Strategy>
+    static inline OutputIterator apply(Geometry1 const& g1,
+                                       Geometry2 const& g2,
+                                       RobustPolicy const& robust_policy,
+                                       OutputIterator out,
+                                       Strategy const& strategy)
+    {
+        access1::get(out) = geometry::dispatch::intersection_insert
+            <
+                Geometry1, Geometry2,
+                typename access1::type,
+                overlay_difference,
+                geometry::detail::overlay::do_reverse<geometry::point_order<Geometry1>::value>::value,
+                geometry::detail::overlay::do_reverse<geometry::point_order<Geometry2>::value, true>::value
+            >::apply(g1, g2, robust_policy, access1::get(out), strategy);
+
+        access2::get(out) = geometry::detail::convert_to_output
+            <
+                Geometry2,
+                typename access2::type
+            >::apply(g2, access2::get(out));
+
+        return out;
+    }
+};
+
+
+template
+<
+    typename Geometry1, typename Geometry2, typename SingleTupledOut,
+    typename SingleTag1, typename SingleTag2
+>
+struct union_insert_tupled_different
+    <
+        Geometry1, Geometry2, SingleTupledOut, SingleTag1, SingleTag2, false
+    >
+{
+    template <typename RobustPolicy, typename OutputIterator, typename Strategy>
+    static inline OutputIterator apply(Geometry1 const& g1,
+                                       Geometry2 const& g2,
+                                       RobustPolicy const& robust_policy,
+                                       OutputIterator out,
+                                       Strategy const& strategy)
+    {
+        return union_insert_tupled_different
+            <
+                Geometry2, Geometry1, SingleTupledOut, SingleTag2, SingleTag1, true
+            >::apply(g2, g1, robust_policy, out, strategy);
+    }
+};
+
+
+template
+<
+    typename Geometry1, typename Geometry2, typename SingleTupledOut,
+    typename TagIn1, typename TagIn2,
+    typename CastedTagIn1, typename CastedTagIn2
+>
+struct union_insert
+    <
+        Geometry1, Geometry2, SingleTupledOut,
+        TagIn1, TagIn2, detail::tupled_output_tag,
+        CastedTagIn1, CastedTagIn2, detail::tupled_output_tag,
+        false
+    >
+{
+    typedef typename geometry::detail::single_tag_from_base_tag
+        <
+            CastedTagIn1
+        >::type single_tag1;
+
+    typedef detail::expect_output
+        <
+            Geometry1, Geometry2, SingleTupledOut, single_tag1
+        > expect_check1;
+
+    typedef typename geometry::detail::single_tag_from_base_tag
+        <
+            CastedTagIn2
+        >::type single_tag2;
+
+    typedef detail::expect_output
+        <
+            Geometry1, Geometry2, SingleTupledOut, single_tag2
+        > expect_check2;
+
+    template <typename RobustPolicy, typename OutputIterator, typename Strategy>
+    static inline OutputIterator apply(Geometry1 const& g1,
+                                       Geometry2 const& g2,
+                                       RobustPolicy const& robust_policy,
+                                       OutputIterator out,
+                                       Strategy const& strategy)
+    {
+        return union_insert_tupled_different
+            <
+                Geometry1, Geometry2, SingleTupledOut, single_tag1, single_tag2
+            >::apply(g1, g2, robust_policy, out, strategy);
+    }
+};
+
+
 } // namespace dispatch
 #endif // DOXYGEN_NO_DISPATCH
 
@@ -209,7 +347,7 @@ inline OutputIterator union_insert(Geometry1 const& geometry1,
 {
     concepts::check<Geometry1 const>();
     concepts::check<Geometry2 const>();
-    concepts::check<GeometryOut>();
+    geometry::detail::output_geometry_concept_check<GeometryOut>::apply();
 
     typename strategy::intersection::services::default_strategy
         <
@@ -253,7 +391,10 @@ struct union_
                              Collection & output_collection,
                              Strategy const& strategy)
     {
-        typedef typename boost::range_value<Collection>::type geometry_out;
+        typedef typename geometry::detail::output_geometry_value
+            <
+                Collection
+            >::type single_out;
 
         typedef typename geometry::rescale_overlay_policy_type
             <
@@ -268,9 +409,9 @@ struct union_
 
         dispatch::union_insert
            <
-               Geometry1, Geometry2, geometry_out
+               Geometry1, Geometry2, single_out
            >::apply(geometry1, geometry2, robust_policy,
-                    range::back_inserter(output_collection),
+                    geometry::detail::output_geometry_back_inserter(output_collection),
                     strategy);
     }
 
@@ -312,7 +453,14 @@ struct union_
     {
         concepts::check<Geometry1 const>();
         concepts::check<Geometry2 const>();
-        concepts::check<typename boost::range_value<Collection>::type>();
+        //concepts::check<typename boost::range_value<Collection>::type>();
+        geometry::detail::output_geometry_concept_check
+            <
+                typename geometry::detail::output_geometry_value
+                    <
+                        Collection
+                    >::type
+            >::apply();
 
         resolve_strategy::union_::apply(geometry1, geometry2,
                                         output_collection,
diff --git a/include/boost/geometry/arithmetic/infinite_line_functions.hpp b/include/boost/geometry/arithmetic/infinite_line_functions.hpp
index 0e82814831..a4496b9b04 100644
--- a/include/boost/geometry/arithmetic/infinite_line_functions.hpp
+++ b/include/boost/geometry/arithmetic/infinite_line_functions.hpp
@@ -76,19 +76,6 @@ side_value(model::infinite_line<Type> const& line, Point const& p)
     return side_value(line, geometry::get<0>(p), geometry::get<1>(p));
 }
 
-// Returns true for two lines which are supposed to be (close to) collinear
-// (which is not checked) and have a similar direction
-// (in practice up to 45 degrees, TO BE VERIFIED)
-// true: -----------------> p -----------------> q
-// false: -----------------> p <----------------- q
-template <typename Type>
-inline
-bool similar_direction(const model::infinite_line<Type>& p,
-                       const model::infinite_line<Type>& q)
-{
-    return p.a * q.a >= 0 && p.b * q.b >= 0;
-}
-
 template <typename Type>
 inline bool is_degenerate(const model::infinite_line<Type>& line)
 {
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 <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
+{
+
+template <typename Filter, std::size_t N, std::size_t End>
+struct make_filter_impl
+{
+    template <typename ExtremaArray, typename ...Reals>
+    static Filter apply(const ExtremaArray& extrema, const Reals&... args)
+    {
+        return make_filter_impl<Filter, N + 1, End>
+            ::apply(extrema, args..., extrema[N]);
+    }
+};
+
+template <typename Filter, std::size_t End>
+struct make_filter_impl<Filter, End, End>
+{
+    template <typename ExtremaArray, typename ...Reals>
+    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<Expression>;
+    using ct = CalculationType;
+    using extrema_array = std::array<ct, 2 * expression_max_argn::value>;
+    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<ct>::infinity());
+        std::fill(m_extrema.begin() + m_extrema.size() / 2,
+                  m_extrema.end(),
+                  std::numeric_limits<ct>::infinity());
+    }
+    template <typename ...Reals>
+    int apply(const Reals&... args) const
+    {
+        return m_filter.apply(args...);
+    }
+
+    template <typename ...Reals>
+    inline void update_extrema(const Reals&... args)
+    {
+        std::array<ct, sizeof...(Reals)> input {{ static_cast<ct>(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 <typename ...Reals>
+    inline bool update_extrema_check(const Reals&... args)
+    {
+        bool changed = false;
+        std::array<ct, sizeof...(Reals)> input {{ static_cast<ct>(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 <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
+{
+
+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<Node, 0, Real, InputArr>
+{
+    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<Node, Node::argn, Real, InputArr>::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<All, Node>::value];
+    }
+};
+
+template
+<
+    typename All,
+    typename Node,
+    typename Real,
+    typename Arr,
+    typename InputArr
+>
+struct get_approx_impl<All, Node, Real, Arr, true, InputArr>
+{
+    static inline Real apply(Arr& interim_results, const InputArr& input)
+    {
+        return get_nth_real<Node, Node::argn, Real, InputArr>(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<Node>::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<Remaining>;
+        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<Remaining>::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<Remaining>;
+        interim_results[boost::mp11::mp_find<All, node>::value] =
+                  get_approx<All, typename node::left, Real>(interim_results,
+                                                             input)
+                * get_approx<All, typename node::right, Real>(interim_results,
+                                                              input);
+        approximate_remainder
+            <
+                All,
+                boost::mp11::mp_pop_front<Remaining>,
+                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<Remaining>;
+        interim_results[boost::mp11::mp_find<All, node>::value] = std::max(
+                  get_approx<All, typename node::left, Real>(interim_results,
+                                                             input),
+                  get_approx<All, typename node::right, Real>(interim_results,
+                                                              input));
+        approximate_remainder
+            <
+                All,
+                boost::mp11::mp_pop_front<Remaining>,
+                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<Remaining>;
+        interim_results[boost::mp11::mp_find<All, node>::value] = std::min(
+                  get_approx<All, typename node::left, Real>(interim_results,
+                                                             input),
+                  get_approx<All, typename node::right, Real>(interim_results,
+                                                              input));
+        approximate_remainder
+            <
+                All,
+                boost::mp11::mp_pop_front<Remaining>,
+                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<Remaining>;
+        interim_results[boost::mp11::mp_find<All, node>::value] =
+                  get_approx<All, typename node::left, Real>(interim_results,
+                                                             input)
+                + get_approx<All, typename node::right, Real>(interim_results,
+                                                              input);
+        approximate_remainder
+            <
+                All,
+                boost::mp11::mp_pop_front<Remaining>,
+                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<Remaining>;
+        interim_results[boost::mp11::mp_find<All, node>::value] =
+                  get_approx<All, typename node::left, Real>(interim_results,
+                                                             input)
+                - get_approx<All, typename node::right, Real>(interim_results,
+                                                              input);
+        approximate_remainder
+            <
+                All,
+                boost::mp11::mp_pop_front<Remaining>,
+                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<Remaining>;
+        interim_results[boost::mp11::mp_find<All, node>::value] =
+            std::abs(get_approx
+                <
+                    All,
+                    typename node::child,
+                    Real
+                >(interim_results, input));
+        approximate_remainder
+            <
+                All,
+                boost::mp11::mp_pop_front<Remaining>,
+                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<Remaining>,
+                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<typename Expression, typename Real, typename InputArr>
+inline Real approximate_value(const InputArr& input)
+{
+    using stack = typename boost::mp11::mp_unique<post_order<Expression>>;
+    using evals = typename boost::mp11::mp_remove_if<stack, is_leaf>;
+    std::array<Real, boost::mp11::mp_size<evals>::value> results;
+        approximate_interim<evals, evals, Real>(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 <cstddef>
+
+#include <boost/mp11/integral.hpp>
+#include <boost/mp11/list.hpp>
+
+#include <boost/geometry/extensions/generic_robust_predicates/strategies/cartesian/detail/static_util.hpp>
+
+namespace boost { namespace geometry
+{
+
+namespace detail { namespace generic_robust_predicates
+{
+
+//only meant to be used in assertions
+template <typename T> using is_mp_int =
+    boost::mp11::mp_same<T, boost::mp11::mp_int<T::value>>;
+
+template <typename L> using is_coeff_list = boost::mp11::mp_and
+    <
+        boost::mp11::mp_all_of<L, is_mp_int>,
+        boost::mp11::mp_similar<L, boost::mp11::mp_list<>>
+    >;
+
+template <typename L>
+struct coeff_truncate_impl
+{
+private:
+    static_assert(is_coeff_list<L>::value, "L must be a coefficient list");
+    using first_nz = boost::mp11::mp_find_if<L, is_not_zero>;
+    using after_second_nz =
+        boost::mp11::mp_min<
+            boost::mp11::mp_size_t<first_nz::value + 2>,
+            boost::mp11::mp_size<L>
+        >;
+    using tail = boost::mp11::mp_erase_c<L, 0, after_second_nz::value>;
+    using head = boost::mp11::mp_erase
+            <
+                L,
+                after_second_nz,
+                boost::mp11::mp_size<L>
+            >;
+    using round_up = boost::mp11::mp_less
+            <
+                boost::mp11::mp_find_if<tail, is_not_zero>,
+                boost::mp11::mp_size<tail>
+            >;
+public:
+    using type = boost::mp11::mp_if
+            <
+                round_up,
+                boost::mp11::mp_push_back<head, boost::mp11::mp_int<1>>,
+                head
+            >;
+    static_assert(is_coeff_list<L>::value, "type should be a coefficient list");
+};
+
+template <typename L> using coeff_truncate = typename coeff_truncate_impl<L>::type;
+
+template <typename L> using app_zero_b =
+    boost::mp11::mp_push_front<L, boost::mp11::mp_int<0>>;
+template <typename L> using app_zero_f =
+    boost::mp11::mp_push_back<L, boost::mp11::mp_int<0>>;
+template <typename L> using mult_by_1_p_eps = coeff_truncate
+        <
+            boost::mp11::mp_transform
+                <
+                    boost::mp11::mp_plus,
+                    app_zero_b<L>,
+                    app_zero_f<L>
+                >
+        >;
+
+template
+<
+    typename L,
+    typename N,
+    typename done = boost::mp11::mp_bool<N::value == 0>
+>
+struct mult_by_1_p_eps_pow_impl
+{
+private:
+    using next = mult_by_1_p_eps<L>;
+public:
+    using type =
+        typename mult_by_1_p_eps_pow_impl
+            <
+                next,
+                boost::mp11::mp_int<N::value - 1>
+            >::type;
+};
+
+template <typename L, typename N>
+struct mult_by_1_p_eps_pow_impl<L, N, boost::mp11::mp_true>
+{
+public:
+    using type = L;
+};
+
+template <typename L, typename N> using mult_by_1_p_eps_pow = coeff_truncate
+        <
+            typename mult_by_1_p_eps_pow_impl<L, N>::type
+        >;
+
+template <typename L> using div_by_1_m_eps_helper =
+    boost::mp11::mp_partial_sum
+        <
+            L,
+            boost::mp11::mp_int<0>,
+            boost::mp11::mp_plus
+        >;
+
+template <typename L> using div_by_1_m_eps = coeff_truncate
+    <
+        boost::mp11::mp_push_back
+            <
+                boost::mp11::mp_pop_back<div_by_1_m_eps_helper<L>>,
+                inc<boost::mp11::mp_back<div_by_1_m_eps_helper<L>>>
+            >
+    >;
+
+template
+<
+    typename L1,
+    typename L2,
+    typename L,
+    typename L1_Empty = typename empty_or_void<L1>::type,
+    typename L2_Empty = typename empty_or_void<L2>::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<L1>,
+            boost::mp11::mp_pop_front<L2>,
+            boost::mp11::mp_push_back
+                <
+                    L,
+                    boost::mp11::mp_plus
+                        <
+                            boost::mp11::mp_front<L1>,
+                            boost::mp11::mp_front<L2>
+                        >
+                >
+        >::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<L, L2>;
+};
+
+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<L, L1>;
+};
+
+template <typename L1, typename L2> 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<L1>::type,
+    typename L2_Empty = typename empty_or_void<L2>::type
+>
+struct coeff_max_impl {};
+
+template
+<
+    typename L1,
+    typename L2,
+    typename L
+>
+struct coeff_max_impl<L1, L2, L, boost::mp11::mp_false, boost::mp11::mp_false>
+{
+    using type = typename coeff_max_impl
+        <
+            boost::mp11::mp_if
+                <
+                    boost::mp11::mp_less
+                        <
+                            boost::mp11::mp_front<L1>,
+                            boost::mp11::mp_front<L2>
+                        >,
+                    boost::mp11::mp_list<>,
+                    boost::mp11::mp_pop_front<L1>
+                >,
+            boost::mp11::mp_if
+                <
+                    boost::mp11::mp_less
+                        <
+                            boost::mp11::mp_front<L2>,
+                            boost::mp11::mp_front<L1>
+                        >,
+                    boost::mp11::mp_list<>,
+                    boost::mp11::mp_pop_front<L2>
+                >,
+            boost::mp11::mp_push_back
+                <
+                    L,
+                    boost::mp11::mp_max
+                        <
+                            boost::mp11::mp_front<L1>,
+                            boost::mp11::mp_front<L2>
+                        >
+                >
+        >::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<L, L2>;
+};
+
+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<L, L1>;
+};
+
+template <typename L1, typename L2> 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<L>::value
+        - boost::mp11::mp_find_if<L, is_not_zero>::value
+>
+struct coeff_round_impl
+{
+private:
+    using first_nz = boost::mp11::mp_find_if<L, is_not_zero>;
+    using tail = boost::mp11::mp_erase_c<L, 0, first_nz::value + 2>;
+    using zero_tail = boost::mp11::mp_same
+        <
+            boost::mp11::mp_find_if<tail, is_not_zero>,
+            boost::mp11::mp_size<tail>
+        >;
+    using head = boost::mp11::mp_erase_c
+        <
+            L,
+            first_nz::value + 2,
+            boost::mp11::mp_size<L>::value
+        >;
+    using minor = boost::mp11::mp_if
+        <
+            zero_tail,
+            boost::mp11::mp_back<head>,
+            inc<boost::mp11::mp_back<head>>
+        >;
+    using major = boost::mp11::mp_at<head, first_nz>;
+    using major_rounded = boost::mp11::mp_int< 1 << log_2_ceil<major>::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<head>,
+            minor_rounded
+        >;
+};
+
+template <typename L>
+struct coeff_round_impl<L, 0> { using type = L; };
+
+template <typename L>
+struct coeff_round_impl<L, 1> { using type = L; };
+
+template <typename L> using coeff_round = typename coeff_round_impl<L>::type;
+
+template <typename IV1, typename IV2> using indexed_value_product =
+    boost::mp11::mp_list
+        <
+            boost::mp11::mp_plus
+                <
+                    boost::mp11::mp_first<IV1>,
+                    boost::mp11::mp_first<IV2>
+                >,
+            boost::mp11::mp_int
+                <
+                      boost::mp11::mp_second<IV1>::value
+                    * boost::mp11::mp_second<IV2>::value
+                >
+        >;
+
+template <typename V1> struct add_nested
+{
+    template <typename K, typename V2> using fn =
+        boost::mp11::mp_plus<V1, V2>;
+};
+
+template <typename M, typename IV> using list_product_fold =
+    boost::mp11::mp_map_update_q
+        <
+            M,
+            IV,
+            add_nested< boost::mp11::mp_second<IV> >
+        >;
+
+template<typename IV1, typename IV2> using index_less =
+    boost::mp11::mp_less
+        <
+            boost::mp11::mp_first<IV1>,
+            boost::mp11::mp_first<IV2>
+        >;
+
+template <typename L1, typename L2> using list_product =
+    strip_index
+        <
+            boost::mp11::mp_sort
+                <
+                    boost::mp11::mp_fold
+                        <
+                            boost::mp11::mp_product
+                                <
+                                    indexed_value_product,
+                                    indexed<L1>,
+                                    indexed<L2>
+                                >,
+                            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 <limits>
+
+#include <boost/mp11/integral.hpp>
+#include <boost/mp11/list.hpp>
+#include <boost/mp11/map.hpp>
+
+#include <boost/geometry/extensions/generic_robust_predicates/strategies/cartesian/detail/coefficient_list.hpp>
+
+namespace boost { namespace geometry
+{
+
+namespace detail { namespace generic_robust_predicates
+{
+
+template <typename KV> using second_is_coeff_list =
+    is_coeff_list< boost::mp11::mp_second<KV> >;
+
+template <typename EM> using is_error_map = boost::mp11::mp_and
+    <
+        boost::mp11::mp_is_map<EM>,
+        boost::mp11::mp_all_of<EM, second_is_coeff_list>,
+        boost::mp11::mp_similar<EM, boost::mp11::mp_list<>
+    >
+    //TODO: It would be desirable to also validate that the keys are good
+>;
+
+template <typename KV, typename M>
+struct error_map_insert_impl
+{
+private:
+    using key = boost::mp11::mp_front<KV>;
+    using value = boost::mp11::mp_second<KV>;
+    using other_value = typename mp_map_at_second_or_void<M, key>::type;
+    using merged_value = coeff_merge<value, other_value>;
+    using nkv = boost::mp11::mp_list<key, merged_value>;
+public:
+    using type = boost::mp11::mp_map_replace<M, nkv>;
+};
+
+template <typename KV, typename M> using error_map_insert =
+    typename error_map_insert_impl<KV, M>::type;
+
+template <typename M1, typename M2>
+struct add_fold_operator
+{
+public:
+    template <typename M, typename K> using fn = boost::mp11::mp_map_insert
+        <
+            M,
+            boost::mp11::mp_list
+                <
+                    K,
+                    coeff_merge
+                        <
+                            typename mp_map_at_second_or_void<M1, K>::type,
+                            typename mp_map_at_second_or_void<M2, K>::type
+                        >
+                >
+        >;
+};
+
+template <typename M1, typename M2> using add_children = boost::mp11::mp_fold
+    <
+        boost::mp11::mp_set_union
+            <
+                boost::mp11::mp_map_keys<M1>,
+                boost::mp11::mp_map_keys<M2>
+            >,
+        boost::mp11::mp_list<>,
+        add_fold_operator<M1, M2>::template fn
+    >;
+
+template
+<
+    typename Exp,
+    typename EM,
+    typename Out,
+    typename LErr = typename val_or_empty_list<EM, typename Exp::left>::type,
+    typename RErr = typename val_or_empty_list<EM, typename Exp::right>::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<boost::mp11::mp_int<1>>
+                >,
+            Out
+        >;
+};
+
+template
+<
+    typename Exp,
+    typename EM,
+    typename Out
+>
+using decompose_add = typename decompose_add_impl<Exp, EM, Out>::type;
+
+template
+<
+    typename Exp,
+    typename EM,
+    typename LErr = typename val_or_empty_list<EM, typename Exp::left>::type,
+    typename RErr = typename val_or_empty_list<EM, typename Exp::right>::type,
+    typename Children_Empty =
+        boost::mp11::mp_and
+            <
+                typename empty_or_void<LErr>::type,
+                typename empty_or_void<RErr>::type
+            >
+>
+struct sum_err_impl
+{
+private:
+    static_assert(is_error_map<LErr>::value, "LErr needs to be a valid error map.");
+    static_assert(is_error_map<RErr>::value, "RErr needs to be a valid error map.");
+    using children = add_children<LErr, RErr>;
+    static_assert(
+            is_error_map<children>::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>>>
+    >;*/
+    using type = decompose_add<Exp, EM, children>;
+    static_assert(is_error_map<type>::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<LErr>::value, "LErr needs to be a valid error map.");
+    static_assert(is_error_map<RErr>::value, "RErr needs to be a valid error map.");
+    using type = boost::mp11::mp_list
+        <
+            boost::mp11::mp_list
+                <
+                    Exp,
+                    boost::mp11::mp_list<boost::mp11::mp_int<1>>
+                >
+        >;
+    static_assert(is_error_map<type>::value, "type needs to be a valid error map.");
+};
+
+template <typename Exp, typename EM> using sum_err =
+    typename sum_err_impl<Exp, EM>::type;
+
+template <typename L> using pad_second = boost::mp11::mp_list
+    <
+        boost::mp11::mp_front<L>,
+        boost::mp11::mp_push_front
+            <
+                boost::mp11::mp_second<L>,
+                boost::mp11::mp_int<0>
+            >
+    >;
+
+template <typename L> using pop_front_second = boost::mp11::mp_list
+    <
+        boost::mp11::mp_front<L>,
+        boost::mp11::mp_pop_front<boost::mp11::mp_second<L>>
+    >;
+
+template <typename K, typename V> using increment_first_of_second =
+    boost::mp11::mp_transform_front<V, inc>;
+
+template <typename KV1, typename KV2> using prod_entry_merge =
+    boost::mp11::mp_list
+        <
+            boost::mp11::mp_flatten
+                <
+                    boost::mp11::mp_list
+                        <
+                            boost::mp11::mp_first<KV1>,
+                            boost::mp11::mp_first<KV2>
+                        >
+                >,
+            list_product
+                <
+                    boost::mp11::mp_second<KV1>,
+                    boost::mp11::mp_second<KV2>
+                >
+        >;
+
+template
+<
+    typename Exp,
+    typename EM,
+    typename LErr = typename val_or_empty_list<EM, typename Exp::left>::type,
+    typename RErr = typename val_or_empty_list<EM, typename Exp::right>::type
+>
+struct prod_children_impl
+{
+private:
+    static_assert(is_error_map<LErr>::value, "LErr needs to be a valid error map.");
+    static_assert(is_error_map<RErr>::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<pad_second, LErr>;
+    using added_left_lerr = decompose_add<left, EM, padded_lerr>;
+    using padded_rerr = boost::mp11::mp_transform<pad_second, RErr>;
+    using added_right_rerr = decompose_add<right, EM, padded_rerr>;
+    using prod = boost::mp11::mp_product
+        <
+            prod_entry_merge,
+            added_left_lerr,
+            added_right_rerr
+        >;
+    using stripped_prod = boost::mp11::mp_transform<pop_front_second, prod>;
+public:
+    using type = stripped_prod;
+    static_assert(is_error_map<type>::value, "type needs to be a valid error map.");
+};
+
+template <typename Exp, typename EM> using prod_children =
+    typename prod_children_impl<Exp, EM>::type;
+
+template
+<
+    typename Exp,
+    typename EM,
+    typename LErr = typename val_or_empty_list<EM, typename Exp::left>::type,
+    typename RErr = typename val_or_empty_list<EM, typename Exp::right>::type
+>
+struct product_err_impl
+{
+private:
+    static_assert(is_error_map<LErr>::value, "LErr needs to be a valid error map.");
+    static_assert(is_error_map<RErr>::value, "RErr needs to be a valid error map.");
+    using children = prod_children<Exp, EM>;
+    static_assert(is_error_map<children>::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<type>::value, "type needs to be a valid error map.");
+};
+
+template <typename Exp, typename EM> using product_err =
+    typename product_err_impl<Exp, EM>::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<Exp, Errors>;
+};
+
+template
+<
+    typename Errors,
+    typename Exp
+>
+struct sum_or_product_err_impl<Errors, Exp, boost::mp11::mp_false>
+{
+    using type = product_err<Exp, Errors>;
+};
+
+template <typename Errors, typename Exp> using sum_or_product_err =
+    typename sum_or_product_err_impl<Errors, Exp>::type;
+
+template
+<
+    typename Errors,
+    typename Exp
+>
+struct error_fold_impl
+{
+private:
+    using lerr = typename val_or_empty_list<Errors, typename Exp::left>::type;
+    using rerr = typename val_or_empty_list<Errors, typename Exp::right>::type;
+    using err = sum_or_product_err<Errors, Exp>;
+public:
+    using type = boost::mp11::mp_map_insert
+        <
+            Errors,
+            boost::mp11::mp_list<Exp, err>
+        >;
+};
+
+
+template <typename Errors, typename Exp> using error_fold =
+    typename error_fold_impl<Errors, Exp>::type;
+
+template <typename Evals> using evals_error =
+    boost::mp11::mp_fold
+        <
+            Evals,
+            boost::mp11::mp_list<>,
+            error_fold
+        >;
+
+template <typename T> using is_mp_list = boost::mp11::mp_similar
+    <
+        boost::mp11::mp_list<>,
+        T
+    >;
+
+template <typename KV>
+struct list_to_product_impl
+{
+private:
+    using key = boost::mp11::mp_front<KV>;
+    using value = boost::mp11::mp_second<KV>;
+    using multiplications = boost::mp11::mp_int<boost::mp11::mp_size<key>::value - 1>;
+    using nvalue = mult_by_1_p_eps_pow<value, multiplications>;
+    using nkey = boost::mp11::mp_reverse_fold<
+        boost::mp11::mp_pop_back<key>,
+        boost::mp11::mp_back<key>,
+        product
+    >;
+public:
+    using type = boost::mp11::mp_list<nkey, nvalue>;
+};
+
+template <typename KV> using list_to_product = typename list_to_product_impl<KV>::type;
+
+template
+<
+    typename M,
+    typename KV,
+    typename KeyMPList = is_mp_list< boost::mp11::mp_front<KV> >
+>
+struct error_map_list_to_product_fold_impl
+{
+    using type = error_map_insert<list_to_product<KV>, M>;
+};
+
+template <typename M, typename KV>
+struct error_map_list_to_product_fold_impl
+    <
+        M,
+        KV,
+        boost::mp11::mp_false
+    >
+{
+    using type = error_map_insert<KV, M>;
+};
+
+template <typename M, typename KV> using error_map_list_to_product_fold =
+    typename error_map_list_to_product_fold_impl<M, KV>::type;
+
+template <typename M>
+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 <typename M> using error_map_list_to_product =
+    typename error_map_list_to_product_impl<M>::type;
+
+template <typename Expression>
+struct abs_unless_non_negative
+{
+    using type = boost::mp11::mp_if_c
+        <
+            Expression::non_negative,
+            Expression,
+            abs<Expression>
+        >;
+};
+
+template <typename EM, typename KV> using abs_fold =
+    boost::mp11::mp_push_back
+        <
+            EM,
+            boost::mp11::mp_list
+                <
+                    typename abs_unless_non_negative<boost::mp11::mp_front<KV>>
+                        ::type,
+                    boost::mp11::mp_second<KV>
+                >
+        >;
+
+template<typename EM> 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<KV1>;
+    using key2 = boost::mp11::mp_front<KV2>;
+    using nkey = sum<key1, key2>;
+    using val1 = boost::mp11::mp_second<KV1>;
+    using val2 = boost::mp11::mp_second<KV2>;
+    using mval = coeff_max<val1, val2>;
+    static_assert(is_coeff_list<mval>::value, "merged coefficient list fails coefficient list check");
+    using nval = mult_by_1_p_eps<mval>;
+public:
+    using type = boost::mp11::mp_list<nkey, nval>;
+};
+
+template<typename KV1, typename KV2> using abs_sum_error_term = typename abs_sum_error_term_impl<KV1, KV2>::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<typename M> using error_map_sum_up =
+    boost::mp11::mp_fold
+        <
+            boost::mp11::mp_pop_front<M>,
+            boost::mp11::mp_first<M>,
+            abs_sum_error_term
+        >;
+
+template<typename M, std::size_t MS = boost::mp11::mp_size<M>::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<M, 0, mid>
+        >::type;
+    using right = typename error_map_balanced_sum_up_impl
+        <
+            boost::mp11::mp_erase_c<M, mid, MS>
+        >::type;
+public:
+    using type = abs_sum_error_term<left, right>;
+};
+
+template<typename M>
+struct error_map_balanced_sum_up_impl<M, 1>
+{
+    using type = boost::mp11::mp_front<M>;
+};
+
+template<typename M> using error_map_balanced_sum_up =
+    typename error_map_balanced_sum_up_impl<M>::type;
+
+template
+<
+    typename Real,
+    typename Exp
+>
+struct eps_pow
+{
+    static constexpr Real value =
+          std::numeric_limits<Real>::epsilon()/2.0
+        * eps_pow<Real, boost::mp11::mp_size_t< Exp::value - 1 >>::value;
+};
+
+template <typename Real> struct eps_pow<Real, boost::mp11::mp_size_t<0>>
+{
+    static constexpr Real value = 1.0;
+};
+
+template
+<
+    typename Real,
+    typename L,
+    typename S = boost::mp11::mp_size<L>
+>
+struct eval_eps_polynomial
+{
+private:
+    using last = boost::mp11::mp_back<L>;
+    using s2last = boost::mp11::mp_back< boost::mp11::mp_pop_back<L> >;
+public:
+    static constexpr Real value =
+          s2last::value
+        * eps_pow<Real, boost::mp11::mp_size_t<S::value - 1>>::value
+        + last::value * eps_pow<Real, S>::value;
+};
+
+template
+<
+    typename Real,
+    typename L
+>
+struct eval_eps_polynomial<Real, L, boost::mp11::mp_size_t<1>>
+{
+    static constexpr Real value =
+          boost::mp11::mp_front<L>::value
+        * std::numeric_limits<Real>::epsilon() / 2.0;
+};
+
+template
+<
+    typename Real,
+    typename L
+>
+struct eval_eps_polynomial<Real, L, boost::mp11::mp_size_t<0>>
+{
+    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 <cstddef>
+#include <type_traits>
+
+#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
+{
+
+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;
+    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>;
+
+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 Node, typename IsLeaf = is_leaf<Node>>
+struct max_argn_impl;
+
+template <typename Node> using max_argn = typename max_argn_impl<Node>::type;
+
+template <typename Node>
+struct max_argn_impl<Node, boost::mp11::mp_false>
+{
+private:
+    using children_list = boost::mp11::mp_rename<Node, boost::mp11::mp_list>;
+    using children_max_argn =
+        boost::mp11::mp_transform<max_argn, children_list>;
+public:
+    using type = boost::mp11::mp_max_element
+        <
+            children_max_argn,
+            boost::mp11::mp_less
+        >;
+};
+
+template <typename Node>
+struct max_argn_impl<Node, boost::mp11::mp_true>
+{
+    using type = boost::mp11::mp_size_t<Node::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>;
+
+}} // 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 <boost/geometry/extensions/generic_robust_predicates/strategies/cartesian/detail/expression_tree.hpp>
+
+namespace boost { namespace geometry
+{
+
+namespace detail { namespace generic_robust_predicates
+{
+
+template
+<
+    typename A11, typename A12,
+    typename A21, typename A22
+>
+using det2x2 = difference
+    <
+        product<A11, A22>,
+        product<A12, A21>
+    >;
+
+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<A11, det2x2<A22, A23, A32, A33>>;
+    using minor2 = product<A21, det2x2<A12, A13, A32, A33>>;
+    using minor3 = product<A31, det2x2<A12, A13, A22, A23>>;
+public:
+    using type = sum<minor1, sum<minor2, minor3>>;
+};
+
+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<adx, adx>, product<ady, ady>>;
+    using blift = sum<product<bdx, bdx>, product<bdy, bdy>>;
+    using clift = sum<product<cdx, cdx>, product<cdy, cdy>>;
+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 <boost/geometry/extensions/generic_robust_predicates/strategies/cartesian/detail/expression_tree.hpp>
+
+namespace boost { namespace geometry
+{
+
+namespace detail { namespace generic_robust_predicates
+{
+
+template <typename Expression, std::size_t max_argn>
+struct interval_min_impl
+{
+    using type = Expression;
+};
+
+template <typename Expression, std::size_t max_argn>
+struct interval_max_impl
+{
+    using type = Expression;
+};
+
+template <typename Left, typename Right, std::size_t max_argn>
+struct interval_min_impl<difference<Left, Right>, max_argn>
+{
+private:
+    using min_left = typename interval_min_impl<Left, max_argn>::type;
+    using max_right = typename interval_max_impl<Right, max_argn>::type;
+public:
+    using type = difference<min_left, max_right>;
+};
+
+template <typename Left, typename Right, std::size_t max_argn>
+struct interval_min_impl<sum<Left, Right>, max_argn>
+{
+private:
+    using min_left = typename interval_min_impl<Left, max_argn>::type;
+    using min_right = typename interval_min_impl<Right, max_argn>::type;
+public:
+    using type = sum<min_left, min_right>;
+};
+
+template <typename Left, typename Right, std::size_t max_argn>
+struct interval_max_impl<difference<Left, Right>, max_argn>
+{
+private:
+    using max_left = typename interval_max_impl<Left, max_argn>::type;
+    using min_right = typename interval_min_impl<Right, max_argn>::type;
+public:
+    using type = difference<max_left, min_right>;
+};
+
+template <typename Left, typename Right, std::size_t max_argn>
+struct interval_max_impl<sum<Left, Right>, max_argn>
+{
+private:
+    using max_left = typename interval_max_impl<Left, max_argn>::type;
+    using max_right = typename interval_max_impl<Right, max_argn>::type;
+public:
+    using type = sum<max_left, max_right>;
+};
+
+template <typename Left, typename Right, std::size_t max_argn>
+struct interval_min_impl<product<Left, Right>, max_argn>
+{
+private:
+    using min_left = typename interval_min_impl<Left, max_argn>::type;
+    using max_left = typename interval_max_impl<Left, max_argn>::type;
+    using min_right = typename interval_min_impl<Right, max_argn>::type;
+    using max_right = typename interval_max_impl<Right, max_argn>::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 max_argn>
+struct interval_min_impl<product<Child, Child>, max_argn>
+{
+private:
+    using min_child = typename interval_min_impl<Child, max_argn>::type;
+    using max_child = typename interval_max_impl<Child, max_argn>::type;
+public:
+    using type = min
+        <
+            product<min_child, min_child>,
+            product<max_child, max_child>
+        >;
+};
+
+template <typename Left, typename Right, std::size_t max_argn>
+struct interval_max_impl<product<Left, Right>, max_argn>
+{
+private:
+    using min_left = typename interval_min_impl<Left, max_argn>::type;
+    using max_left = typename interval_max_impl<Left, max_argn>::type;
+    using min_right = typename interval_min_impl<Right, max_argn>::type;
+    using max_right = typename interval_max_impl<Right, max_argn>::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 max_argn>
+struct interval_max_impl<product<Child, Child>, max_argn>
+{
+private:
+    using min_child = typename interval_min_impl<Child, max_argn>::type;
+    using max_child = typename interval_max_impl<Child, max_argn>::type;
+public:
+    using type = max
+        <
+            product<min_child, min_child>,
+            product<max_child, max_child>
+        >;
+};
+
+template <typename Child, std::size_t max_argn>
+struct interval_min_impl<abs<Child>, max_argn>
+{
+private:
+    using min_child = typename interval_min_impl<Child, max_argn>::type;
+    using max_child = typename interval_max_impl<Child, max_argn>::type;
+public:
+    using type = min<abs<min_child>, abs<max_child>>;
+};
+
+template <typename Child, std::size_t max_argn>
+struct interval_max_impl<abs<Child>, max_argn>
+{
+private:
+    using min_child = typename interval_min_impl<Child, max_argn>::type;
+    using max_child = typename interval_max_impl<Child, max_argn>::type;
+public:
+    using type = max<abs<min_child>, abs<max_child>>;
+};
+
+template <std::size_t argn, std::size_t max_argn>
+struct interval_min_impl<argument<argn>, max_argn>
+{
+    using type = argument<argn + max_argn>;
+};
+
+template <std::size_t argn, std::size_t max_argn>
+struct interval_max_impl<argument<argn>, max_argn>
+{
+    using type = argument<argn>;
+};
+
+template <typename Expression, std::size_t max_argn>
+struct interval_impl
+{
+    using type = Expression;
+};
+
+template <typename Child, std::size_t max_argn>
+struct interval_impl<abs<Child>, max_argn>
+{
+private:
+    using min_child = typename interval_min_impl<Child, max_argn>::type;
+    using max_child = typename interval_max_impl<Child, max_argn>::type;
+public:
+    using type = max<abs<min_child>, abs<max_child>>;
+};
+
+template <typename Left, typename Right, std::size_t max_argn>
+struct interval_impl<difference<Left, Right>, max_argn>
+{
+private:
+    using left = typename interval_impl<Left, max_argn>::type;
+    using right = typename interval_impl<Right, max_argn>::type;
+public:
+    using type = difference<left, right>;
+};
+
+template <typename Left, typename Right, std::size_t max_argn>
+struct interval_impl<sum<Left, Right>, max_argn>
+{
+private:
+    using left = typename interval_impl<Left, max_argn>::type;
+    using right = typename interval_impl<Right, max_argn>::type;
+public:
+    using type = sum<left, right>;
+};
+
+template <typename Left, typename Right, std::size_t max_argn>
+struct interval_impl<product<Left, Right>, max_argn>
+{
+private:
+    using left = typename interval_impl<Left, max_argn>::type;
+    using right = typename interval_impl<Right, max_argn>::type;
+public:
+    using type = product<left, right>;
+};
+
+template <typename Expression, std::size_t max_argn>
+using interval_min = typename interval_min_impl<Expression, max_argn>::type;
+
+template <typename Expression, std::size_t max_argn>
+using interval_max = typename interval_max_impl<Expression, max_argn>::type;
+
+template <typename Expression>
+using interval =
+    typename interval_impl<Expression, max_argn<Expression>::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 <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/approximate.hpp>
+
+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<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:
+    template <typename ...Reals>
+    static inline int apply(const Reals&... args)
+    {
+        std::array<CalculationType, sizeof...(Reals)> input
+            {{ static_cast<ct>(args)... }};
+        std::array<ct, boost::mp11::mp_size<all_evals>::value> results;
+        approximate_interim<all_evals, all_evals, ct>(results, input);
+        const ct error_bound =
+            get_approx<all_evals, ErrorExpression, ct>(results, input);
+        const ct det =
+            get_approx<all_evals, Expression, ct>(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 <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/coefficient_list.hpp>
+#include <boost/geometry/extensions/generic_robust_predicates/strategies/cartesian/detail/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
+{
+
+template <typename Expression, typename CalculationType>
+struct stage_a_error_bound_impl
+{
+private:
+    using ct = CalculationType;
+    using root = Expression;
+    using stack = typename boost::mp11::mp_unique<post_order<Expression>>;
+    using evals = typename boost::mp11::mp_remove_if<stack, is_leaf>;
+    using interim_evals = typename boost::mp11::mp_remove
+        <
+            boost::mp11::mp_remove_if<stack, is_leaf>,
+            root
+        >;
+    using interim_errors = evals_error<interim_evals>;
+    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<final_children>;
+    using final_children_ltp_abs = abs_all<final_children_ltp>;
+    using final_children_sum_fold = error_map_balanced_sum_up<final_children_ltp_abs>;
+    using final_coeff = coeff_round
+        <
+            div_by_1_m_eps
+                <
+                    mult_by_1_p_eps
+                        <
+                            boost::mp11::mp_second<final_children_sum_fold>
+                        >
+                >
+        >;
+    static constexpr ct final_coeff_value = eval_eps_polynomial<ct, final_coeff>::value;
+
+    struct final_coeff_constant : public static_constant_interface<ct>
+    {
+        static constexpr ct value = final_coeff_value;
+        static constexpr bool non_negative = true;
+    };
+
+    using error_expression_variable = boost::mp11::mp_front<final_children_sum_fold>;
+public:
+    using type = product<final_coeff_constant, error_expression_variable>;
+};
+
+template
+<
+    typename Expression,
+    typename CalculationType
+>
+using stage_a_error_bound =
+    typename stage_a_error_bound_impl<Expression, CalculationType>::type;
+
+template
+<
+    typename Expression,
+    typename CalculationType
+>
+using stage_a_semi_static = semi_static_filter
+        <
+            Expression,
+            CalculationType,
+            stage_a_error_bound<Expression, CalculationType>
+        >;
+
+template
+<
+    typename Expression,
+    typename CalculationType
+>
+using stage_a_static = static_filter
+        <
+            Expression,
+            CalculationType,
+            interval<stage_a_error_bound<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
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 <array>
+
+#include <boost/mp11/list.hpp>
+#include <boost/mp11/algorithm.hpp>
+#include <boost/mp11/map.hpp>
+#include <boost/mp11/set.hpp>
+
+#include <boost/geometry/extensions/generic_robust_predicates/strategies/cartesian/detail/stage_a.hpp>
+
+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<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>
+        >;
+    ct m_error_bound;
+public:
+    inline static_filter() {}
+
+    ct error_bound() const { return m_error_bound; }
+
+    template <typename ...Reals>
+    inline static_filter(const Reals&... args)
+        : m_error_bound(approximate_value<ErrorExpression, ct>(
+                std::array<ct, sizeof...(Reals)>
+                    {static_cast<ct>(args)...}))
+    {
+        static_assert(sizeof...(Reals) == max_argn<ErrorExpression>::value,
+                      "Number of constructor arguments is incompatible with error expression.");
+    }
+
+    template <typename ...Reals>
+    inline int apply(const Reals&... args) const
+    {
+        std::array<ct, sizeof...(Reals)> input {static_cast<ct>(args)...};
+        std::array<ct, boost::mp11::mp_size<evals>::value> results;
+        approximate_interim<evals, evals, ct>(results, input);
+        const ct det = get_approx<evals, Expression, ct>(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 <typename ...Reals>
+    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 <boost/mp11/integral.hpp>
+#include <boost/mp11/list.hpp>
+#include <boost/mp11/map.hpp>
+
+
+namespace boost { namespace geometry
+{
+
+namespace detail { namespace generic_robust_predicates
+{
+
+template
+<
+    typename LOV,
+    typename Void = boost::mp11::mp_same<LOV, void>
+>
+struct empty_or_void
+{
+    using type = boost::mp11::mp_empty<LOV>;
+};
+
+template <typename LOV>
+struct empty_or_void<LOV, boost::mp11::mp_true>
+{
+    using type = boost::mp11::mp_true;
+};
+
+template <typename N>
+struct log_2_floor_impl
+{
+    using type = typename boost::mp11::mp_plus
+        <
+            typename log_2_floor_impl<boost::mp11::mp_int<N::value / 2>>::type,
+            boost::mp11::mp_int<1>
+        >;
+};
+
+template <>
+struct log_2_floor_impl<boost::mp11::mp_int<1>>
+{
+    using type = boost::mp11::mp_int<0>;
+};
+
+template <typename N> using log_2_floor = typename log_2_floor_impl<N>::type;
+
+template <typename N>
+struct log_2_ceil_impl
+{
+private:
+    using floor = log_2_floor<N>;
+public:
+    using type = boost::mp11::mp_int
+        <
+            (1 << floor::value) == N::value ? floor::value : floor::value + 1
+        >;
+};
+
+template <typename N> using log_2_ceil = typename log_2_ceil_impl<N>::type;
+
+template <typename S> using is_not_zero = boost::mp11::mp_bool<S::value != 0>;
+
+template
+<
+    typename M,
+    typename K,
+    typename Contained = boost::mp11::mp_map_contains<M, K>
+>
+struct mp_map_at_second_or_void
+{
+    using type = void;
+};
+
+template<typename M, typename K>
+struct mp_map_at_second_or_void<M, K, boost::mp11::mp_true>
+{
+    using type = boost::mp11::mp_second<boost::mp11::mp_map_find<M, K>>;
+};
+
+template <typename T> using inc = boost::mp11::mp_int<T::value + 1>;
+
+template <typename L> using indexed = boost::mp11::mp_transform
+    <
+        boost::mp11::mp_list,
+        boost::mp11::mp_iota<boost::mp11::mp_size<L>>,
+        L
+    >;
+
+template <typename L> using strip_index =
+    boost::mp11::mp_transform<boost::mp11::mp_second, L>;
+
+template
+<
+    typename Map,
+    typename Key,
+    typename Contains = boost::mp11::mp_map_contains<Map, Key>
+>
+struct val_or_empty_list
+{
+    using type = boost::mp11::mp_second<boost::mp11::mp_map_find<Map, Key>>;
+};
+
+template
+<
+    typename Map,
+    typename Key
+>
+struct val_or_empty_list<Map, Key, boost::mp11::mp_false>
+{
+    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
diff --git a/include/boost/geometry/extensions/triangulation/strategies/cartesian/detail/precise_math.hpp b/include/boost/geometry/extensions/triangulation/strategies/cartesian/detail/precise_math.hpp
index 7c4a416272..c707620970 100644
--- a/include/boost/geometry/extensions/triangulation/strategies/cartesian/detail/precise_math.hpp
+++ b/include/boost/geometry/extensions/triangulation/strategies/cartesian/detail/precise_math.hpp
@@ -88,42 +88,6 @@ inline std::array<RealNumber, 2> two_diff(RealNumber const a,
     return {{ x, y }};
 }
 
-// constexpr power-method, helper for splitter
-template
-<
-    typename RealNumber
->
-constexpr RealNumber int_pow(RealNumber const base,
-                             int exp,
-                             RealNumber out = 1.0)
-{
-    return exp < 1 ? out :
-        int_pow<RealNumber>(base*base, exp/2, (exp % 2) ? out*base : out);
-}
-
-// consexpr method to compute 2^s + 1 as in Theorem 17, page 18
-template
-<
-    typename RealNumber
->
-constexpr RealNumber splitter()
-{
-    return int_pow<RealNumber>(2.0,
-            (std::numeric_limits<RealNumber>::digits + 1) / 2) + 1;
-}
-
-// see theorem 17, page 18
-template
-<
-    typename RealNumber
->
-inline std::array<RealNumber, 2> split(RealNumber const a) {
-    RealNumber c = splitter<RealNumber>() * a;
-    RealNumber a_big = c - a;
-    RealNumber a_hi = c - a_big;
-    return {{ a_hi, a - a_hi }};
-}
-
 // see theorem 18, page 19
 template
 <
@@ -133,12 +97,7 @@ inline RealNumber two_product_tail(RealNumber const a,
                                    RealNumber const b,
                                    RealNumber const x)
 {
-    std::array<RealNumber, 2> a_expansion = split(a);
-    std::array<RealNumber, 2> b_expansion = split(b);
-    RealNumber err1 = x - (a_expansion[0] * b_expansion[0]);
-    RealNumber err2 = err1 - (a_expansion[1] * b_expansion[0]);
-    RealNumber err3 = err2 - (a_expansion[0] * b_expansion[1]);
-    return (a_expansion[1] * b_expansion[1]) - err3;
+    return std::fma(a, b, -x);
 }
 
 // see theorem 18, page 19
@@ -199,45 +158,63 @@ inline int fast_expansion_sum_zeroelim(
 {
     std::array<RealNumber, 2> Qh;
     int i_e = 0, i_f = 0, i_h = 0;
-    if (std::abs(f[0]) > std::abs(e[0])) {
+    if (std::abs(f[0]) > std::abs(e[0]))
+    {
         Qh[0] = e[i_e++];
-    } else {
+    }
+    else
+    {
         Qh[0] = f[i_f++];
     }
     i_h = 0;
-    if ((i_e < m) && (i_f < n)) {
-        if (std::abs(f[i_f]) > std::abs(e[i_e])) {
+    if ((i_e < m) && (i_f < n))
+    {
+        if (std::abs(f[i_f]) > std::abs(e[i_e]))
+        {
             Qh = fast_two_sum(e[i_e++], Qh[0]);
-        } else {
+        }
+        else
+        {
             Qh = fast_two_sum(f[i_f++], Qh[0]);
         }
-        if (Qh[1] != 0.0) {
+        if (Qh[1] != 0.0)
+        {
             h[i_h++] = Qh[1];
         }
-        while ((i_e < m) && (i_f < n)) {
-            if (std::abs(f[i_f]) > std::abs(e[i_e])) {
+        while ((i_e < m) && (i_f < n))
+        {
+            if (std::abs(f[i_f]) > std::abs(e[i_e]))
+            {
                 Qh = two_sum(Qh[0], e[i_e++]);
-            } else {
+            }
+            else
+            {
                 Qh = two_sum(Qh[0], f[i_f++]);
             }
-            if (Qh[1] != 0.0) {
+            if (Qh[1] != 0.0)
+            {
                 h[i_h++] = Qh[1];
             }
         }
     }
-    while (i_e < m) {
+    while (i_e < m)
+    {
         Qh = two_sum(Qh[0], e[i_e++]);
-        if (Qh[1] != 0.0) {
+        if (Qh[1] != 0.0)
+        {
             h[i_h++] = Qh[1];
         }
     }
-    while (i_f < n) {
+    while (i_f < n)
+    {
         Qh = two_sum(Qh[0], f[i_f++]);
-        if (Qh[1] != 0.0) {
+        if (Qh[1] != 0.0)
+        {
             h[i_h++] = Qh[1];
         }
     }
-    if ((Qh[0] != 0.0) || (i_h == 0)) {
+    if ((Qh[0] != 0.0) || (i_h == 0))
+    {
         h[i_h++] = Qh[0];
     }
     return i_h;
@@ -259,83 +236,76 @@ inline int scale_expansion_zeroelim(
 {
     std::array<RealNumber, 2> Qh = two_product(e[0], b);
     int i_h = 0;
-    if (Qh[1] != 0) {
+    if (Qh[1] != 0)
+    {
         h[i_h++] = Qh[1];
     }
-    for (int i_e = 1; i_e < e_non_zeros; i_e++) {
+    for (int i_e = 1; i_e < e_non_zeros; i_e++)
+    {
         std::array<RealNumber, 2> Tt = two_product(e[i_e], b);
         Qh = two_sum(Qh[0], Tt[1]);
-        if (Qh[1] != 0) {
+        if (Qh[1] != 0)
+        {
             h[i_h++] = Qh[1];
         }
         Qh = fast_two_sum(Tt[0], Qh[0]);
-        if (Qh[1] != 0) {
+        if (Qh[1] != 0)
+        {
             h[i_h++] = Qh[1];
         }
     }
-    if ((Qh[0] != 0.0) || (i_h == 0)) {
+    if ((Qh[0] != 0.0) || (i_h == 0))
+    {
         h[i_h++] = Qh[0];
     }
     return i_h;
 }
 
-// see page 38, Figure 21 for the calculations, notation follows the notation in the figure.
+template<typename RealNumber>
+struct vec2d
+{
+    RealNumber x;
+    RealNumber y;
+};
+
 template
 <
     typename RealNumber,
-    std::size_t Robustness = 3
+    std::size_t Robustness
 >
-inline RealNumber orient2d(std::array<RealNumber, 2> const& p1,
-                           std::array<RealNumber, 2> const& p2,
-                           std::array<RealNumber, 2> const& p3)
+inline RealNumber orient2dtail(vec2d<RealNumber> const& p1,
+                               vec2d<RealNumber> const& p2,
+                               vec2d<RealNumber> const& p3,
+                               std::array<RealNumber, 2>& t1,
+                               std::array<RealNumber, 2>& t2,
+                               std::array<RealNumber, 2>& t3,
+                               std::array<RealNumber, 2>& t4,
+                               std::array<RealNumber, 2>& t5_01,
+                               std::array<RealNumber, 2>& t6_01,
+                               RealNumber const& magnitude
+                               )
 {
-    if(Robustness == 0) {
-        return (p1[0]-p3[0])*(p2[1]-p3[1]) - (p1[1]-p3[1])*(p2[0] - p3[0]);
-    }
-    std::array<RealNumber, 2> t1, t2, t3, t4;
-    t1[0] = p1[0] - p3[0];
-    t2[0] = p2[1] - p3[1];
-    t3[0] = p1[1] - p3[1];
-    t4[0] = p2[0] - p3[0];
-    std::array<RealNumber, 2> t5_01, t6_01;
-    t5_01[0] = t1[0] * t2[0];
-    t6_01[0] = t3[0] * t4[0];
-    RealNumber det = t5_01[0] - t6_01[0];
-    if ( (t5_01[0] > 0 && t6_01[0] <= 0) || (t5_01[0] < 0 && t6_01[0] >= 0) ) {
-        //if diagonal and antidiagonal have different sign, the sign of det is
-        //obvious
-        return det;
-    }
-    RealNumber const magnitude = std::abs(t5_01[0]) + std::abs(t6_01[0]);
-
-    // see p.39, mind the different definition of epsilon for error bound
-    RealNumber const A_relative_bound =
-          (1.5 + 4 * std::numeric_limits<RealNumber>::epsilon())
-        * std::numeric_limits<RealNumber>::epsilon();
-    RealNumber absolute_bound = A_relative_bound * magnitude;
-    if ( std::abs(det) >= absolute_bound ) {
-        return det; //A estimate
-    }
-
     t5_01[1] = two_product_tail(t1[0], t2[0], t5_01[0]);
     t6_01[1] = two_product_tail(t3[0], t4[0], t6_01[0]);
     std::array<RealNumber, 4> tA_03 = two_two_expansion_diff(t5_01, t6_01);
-    det = std::accumulate(tA_03.begin(), tA_03.end(), static_cast<RealNumber>(0));
+    RealNumber det = std::accumulate(tA_03.begin(), tA_03.end(), static_cast<RealNumber>(0));
     if(Robustness == 1) return det;
     // see p.39, mind the different definition of epsilon for error bound
     RealNumber B_relative_bound =
           (1 + 3 * std::numeric_limits<RealNumber>::epsilon())
         * std::numeric_limits<RealNumber>::epsilon();
-    absolute_bound = B_relative_bound * magnitude;
-    if (std::abs(det) >= absolute_bound) {
+    RealNumber absolute_bound = B_relative_bound * magnitude;
+    if (std::abs(det) >= absolute_bound)
+    {
         return det; //B estimate
     }
-    t1[1] = two_diff_tail(p1[0], p3[0], t1[0]);
-    t2[1] = two_diff_tail(p2[1], p3[1], t2[0]);
-    t3[1] = two_diff_tail(p1[1], p3[1], t3[0]);
-    t4[1] = two_diff_tail(p2[0], p3[0], t4[0]);
+    t1[1] = two_diff_tail(p1.x, p3.x, t1[0]);
+    t2[1] = two_diff_tail(p2.y, p3.y, t2[0]);
+    t3[1] = two_diff_tail(p1.y, p3.y, t3[0]);
+    t4[1] = two_diff_tail(p2.x, p3.x, t4[0]);
 
-    if ((t1[1] == 0) && (t3[1] == 0) && (t2[1] == 0) && (t4[1] == 0)) {
+    if ((t1[1] == 0) && (t3[1] == 0) && (t2[1] == 0) && (t4[1] == 0))
+    {
         return det; //If all tails are zero, there is noething else to compute
     }
     RealNumber sub_bound =
@@ -348,7 +318,8 @@ inline RealNumber orient2d(std::array<RealNumber, 2> const& p1,
         * std::numeric_limits<RealNumber>::epsilon();
     absolute_bound = C_relative_bound * magnitude + sub_bound * std::abs(det);
     det += (t1[0] * t2[1] + t2[0] * t1[1]) - (t3[0] * t4[1] + t4[0] * t3[1]);
-    if (Robustness == 2 || std::abs(det) >= absolute_bound) {
+    if (Robustness == 2 || std::abs(det) >= absolute_bound)
+    {
         return det; //C estimate
     }
     std::array<RealNumber, 8> D_left;
@@ -376,6 +347,50 @@ inline RealNumber orient2d(std::array<RealNumber, 2> const& p1,
     return(D[D_nz - 1]);
 }
 
+// see page 38, Figure 21 for the calculations, notation follows the notation in the figure.
+template
+<
+    typename RealNumber,
+    std::size_t Robustness = 3
+>
+inline RealNumber orient2d(vec2d<RealNumber> const& p1,
+                           vec2d<RealNumber> const& p2,
+                           vec2d<RealNumber> const& p3)
+{
+    if(Robustness == 0)
+    {
+        return (p1.x - p3.x) * (p2.y - p3.y) - (p1.y - p3.y) * (p2.x - p3.x);
+    }
+    std::array<RealNumber, 2> t1, t2, t3, t4;
+    t1[0] = p1.x - p3.x;
+    t2[0] = p2.y - p3.y;
+    t3[0] = p1.y - p3.y;
+    t4[0] = p2.x - p3.x;
+    std::array<RealNumber, 2> t5_01, t6_01;
+    t5_01[0] = t1[0] * t2[0];
+    t6_01[0] = t3[0] * t4[0];
+    RealNumber det = t5_01[0] - t6_01[0];
+    RealNumber const magnitude = std::abs(t5_01[0]) + std::abs(t6_01[0]);
+
+    // see p.39, mind the different definition of epsilon for error bound
+    RealNumber const A_relative_bound =
+          (1.5 + 4 * std::numeric_limits<RealNumber>::epsilon())
+        * std::numeric_limits<RealNumber>::epsilon();
+    RealNumber absolute_bound = A_relative_bound * magnitude;
+    if ( std::abs(det) >= absolute_bound )
+    {
+        return det; //A estimate
+    }
+
+    if ( (t5_01[0] > 0 && t6_01[0] <= 0) || (t5_01[0] < 0 && t6_01[0] >= 0) )
+    {
+        //if diagonal and antidiagonal have different sign, the sign of det is
+        //obvious
+        return det;
+    }
+    return orient2dtail<RealNumber, Robustness>(p1, p2, p3, t1, t2, t3, t4, t5_01, t6_01, magnitude);
+}
+
 // This method adaptively computes increasingly precise approximations of the following
 // determinant using Laplace expansion along the last column.
 // det A =
@@ -433,7 +448,8 @@ RealNumber incircle(std::array<RealNumber, 2> const& p1,
           (5 + 24 * std::numeric_limits<RealNumber>::epsilon())
         * std::numeric_limits<RealNumber>::epsilon();
     RealNumber absolute_bound = A_relative_bound * magnitude;
-    if (std::abs(det) > absolute_bound) {
+    if (std::abs(det) > absolute_bound)
+    {
         return det;
     }
     // (p2_x - p4_x) * (p3_y - p4_y)
@@ -492,16 +508,16 @@ RealNumber incircle(std::array<RealNumber, 2> const& p1,
     std::array<RealNumber, 16> C_23_x_A_22_sq;
     // ( (cx - dx) * (ay - dy) - (ax - dx) * (cy - dy) ) * ( by - dy ) * ( by - dy )
     int C_23_x_A_22_sq_nz = scale_expansion_zeroelim(C_23_x_A_22, A_22,
-													 C_23_x_A_22_sq,
-													 C_23_x_A_22_nz);
+                                                     C_23_x_A_22_sq,
+                                                     C_23_x_A_22_nz);
     std::array<RealNumber, 32> A_23_x_C_23;
-	//   ( (cx - dx) * (ay - dy) - (ax - dx) * (cy - dy) )
-	// * ( ( bx - dx ) * ( bx - dx ) + ( by - dy ) * ( by - dy ) )
+    //   ( (cx - dx) * (ay - dy) - (ax - dx) * (cy - dy) )
+    // * ( ( bx - dx ) * ( bx - dx ) + ( by - dy ) * ( by - dy ) )
     int A_23_x_C_23_nz = fast_expansion_sum_zeroelim(C_23_x_A_21_sq,
-													 C_23_x_A_22_sq,
-													 A_23_x_C_23,
-													 C_23_x_A_21_sq_nz,
-													 C_23_x_A_22_sq_nz);
+                                                     C_23_x_A_22_sq,
+                                                     A_23_x_C_23,
+                                                     C_23_x_A_21_sq_nz,
+                                                     C_23_x_A_22_sq_nz);
 
     // (ax - dx) * (by - dy)
     A_11_x_A_22[1] = two_product_tail(A_11, A_22, A_11_x_A_22[0]);
@@ -555,7 +571,8 @@ RealNumber incircle(std::array<RealNumber, 2> const& p1,
           (2 + 12 * std::numeric_limits<RealNumber>::epsilon())
         * std::numeric_limits<RealNumber>::epsilon();
     absolute_bound = B_relative_bound * magnitude;
-    if (std::abs(det) >= absolute_bound) {
+    if (std::abs(det) >= absolute_bound)
+    {
         return det;
     }
     RealNumber A_11tail = two_diff_tail(p1[0], p4[0], A_11);
@@ -565,7 +582,8 @@ RealNumber incircle(std::array<RealNumber, 2> const& p1,
     RealNumber A_31tail = two_diff_tail(p3[0], p4[0], A_31);
     RealNumber A_32tail = two_diff_tail(p3[1], p4[1], A_32);
     if ((A_11tail == 0) && (A_21tail == 0) && (A_31tail == 0)
-        && (A_12tail == 0) && (A_22tail == 0) && (A_32tail == 0)) {
+        && (A_12tail == 0) && (A_22tail == 0) && (A_32tail == 0))
+    {
         return det;
     }
     //  RealNumber sub_bound =  (1.5 + 2.0 * std::numeric_limits<RealNumber>::epsilon())
@@ -583,7 +601,8 @@ RealNumber incircle(std::array<RealNumber, 2> const& p1,
     + ((A_31 * A_31 + A_32 * A_32) * ((A_11 * A_22tail + A_22 * A_11tail)
         - (A_12 * A_21tail + A_21 * A_12tail))
     + 2 * (A_31 * A_31tail + A_32 * A_32tail) * (A_11 * A_22 - A_12 * A_21));
-    //if (std::abs(det) >= absolute_bound) {
+    //if (std::abs(det) >= absolute_bound)
+    //{
     return det;
     //}
 }
diff --git a/include/boost/geometry/extensions/triangulation/strategies/cartesian/side_robust.hpp b/include/boost/geometry/extensions/triangulation/strategies/cartesian/side_robust.hpp
index 3c7d4f8a69..59228c0603 100644
--- a/include/boost/geometry/extensions/triangulation/strategies/cartesian/side_robust.hpp
+++ b/include/boost/geometry/extensions/triangulation/strategies/cartesian/side_robust.hpp
@@ -49,9 +49,10 @@ struct side_robust
     static inline PromotedType side_value(P1 const& p1, P2 const& p2,
         P const& p)
     {
-        std::array<PromotedType, 2> pa {{ get<0>(p1), get<1>(p1) }};
-        std::array<PromotedType, 2> pb {{ get<0>(p2), get<1>(p2) }};
-        std::array<PromotedType, 2> pc {{ get<0>(p), get<1>(p) }};
+        typedef ::boost::geometry::detail::precise_math::vec2d<PromotedType> vec2d;
+        vec2d pa { get<0>(p1), get<1>(p1) };
+        vec2d pb { get<0>(p2), get<1>(p2) };
+        vec2d pc { get<0>(p), get<1>(p) };
         return ::boost::geometry::detail::precise_math::orient2d
             <PromotedType, Robustness>(pa, pb, pc);
     }
diff --git a/include/boost/geometry/index/detail/varray_detail.hpp b/include/boost/geometry/index/detail/varray_detail.hpp
index 31b77c40fe..00aab324a4 100644
--- a/include/boost/geometry/index/detail/varray_detail.hpp
+++ b/include/boost/geometry/index/detail/varray_detail.hpp
@@ -2,7 +2,7 @@
 //
 // varray details
 //
-// Copyright (c) 2012-2015 Adam Wulkiewicz, Lodz, Poland.
+// Copyright (c) 2012-2020 Adam Wulkiewicz, Lodz, Poland.
 // Copyright (c) 2011-2013 Andrew Hundt.
 //
 // Use, modification and distribution is subject to the Boost Software License,
@@ -17,6 +17,8 @@
 #include <memory>
 #include <limits>
 
+#include <boost/config.hpp>
+
 #include <boost/mpl/if.hpp>
 #include <boost/mpl/and.hpp>
 #include <boost/mpl/or.hpp>
@@ -36,8 +38,11 @@
 //#include <boost/type_traits/has_nothrow_assign.hpp>
 //#include <boost/type_traits/has_nothrow_destructor.hpp>
 
+#if ! defined(BOOST_NO_CXX11_HDR_TYPE_TRAITS)
+#include <type_traits>
+#endif
+
 #include <boost/detail/no_exceptions_support.hpp>
-#include <boost/config.hpp>
 #include <boost/move/move.hpp>
 #include <boost/core/addressof.hpp>
 #include <boost/iterator/iterator_traits.hpp>
@@ -53,7 +58,28 @@
 #include <boost/container/vector.hpp>
 #endif // BOOST_GEOMETRY_INDEX_DETAIL_VARRAY_ENABLE_VECTOR_OPTIMIZATION && !BOOST_NO_EXCEPTIONS
 
-namespace boost { namespace geometry { namespace index { namespace detail { namespace varray_detail {
+namespace boost { namespace geometry { namespace index { namespace detail { namespace varray_detail
+{
+
+#if ! defined(BOOST_NO_CXX11_HDR_TYPE_TRAITS)
+
+using std::is_trivially_copyable;
+
+#else
+
+template <typename T>
+struct is_trivially_copyable
+    : boost::integral_constant
+        <
+            bool,
+            boost::has_trivial_copy<T>::value
+         && boost::has_trivial_assign<T>::value
+         && boost::has_trivial_destructor<T>::value
+        >
+{};
+
+#endif
+
 
 template <typename I>
 struct are_elements_contiguous : boost::is_pointer<I>
@@ -164,7 +190,7 @@ template <typename I>
 void destroy(I first, I last)
 {
     typedef typename boost::iterator_value<I>::type value_type;
-    destroy_dispatch(first, last, has_trivial_destructor<value_type>());
+    destroy_dispatch(first, last, boost::has_trivial_destructor<value_type>());
 }
 
 // destroy(I)
@@ -186,7 +212,7 @@ template <typename I>
 void destroy(I pos)
 {
     typedef typename boost::iterator_value<I>::type value_type;
-    destroy_dispatch(pos, has_trivial_destructor<value_type>());
+    destroy_dispatch(pos, boost::has_trivial_destructor<value_type>());
 }
 
 // copy(I, I, O)
@@ -215,7 +241,7 @@ inline O copy(I first, I last, O dst)
     typedef typename
     ::boost::mpl::and_<
         are_corresponding<I, O>,
-        ::boost::has_trivial_assign<
+        is_trivially_copyable<
             typename ::boost::iterator_value<O>::type
         >
     >::type
@@ -253,7 +279,7 @@ F uninitialized_copy(I first, I last, F dst)
     typedef typename
     ::boost::mpl::and_<
         are_corresponding<I, F>,
-        ::boost::has_trivial_copy<
+        is_trivially_copyable<
             typename ::boost::iterator_value<F>::type
         >
     >::type
@@ -308,7 +334,7 @@ O uninitialized_move(I first, I last, O dst)
     typedef typename
     ::boost::mpl::and_<
         are_corresponding<I, O>,
-        ::boost::has_trivial_copy<
+        is_trivially_copyable<
             typename ::boost::iterator_value<O>::type
         >
     >::type
@@ -348,7 +374,7 @@ O move(I first, I last, O dst)
     typedef typename
     ::boost::mpl::and_<
         are_corresponding<I, O>,
-        ::boost::has_trivial_assign<
+        is_trivially_copyable<
             typename ::boost::iterator_value<O>::type
         >
     >::type
@@ -387,7 +413,7 @@ BDO move_backward(BDI first, BDI last, BDO dst)
     typedef typename
     ::boost::mpl::and_<
         are_corresponding<BDI, BDO>,
-        ::boost::has_trivial_assign<
+        is_trivially_copyable<
             typename ::boost::iterator_value<BDO>::type
         >
     >::type
@@ -561,7 +587,7 @@ void construct(DisableTrivialInit const&,
     typedef typename
     ::boost::mpl::and_<
         is_corresponding_value<I, P>,
-        ::boost::has_trivial_copy<P>
+        is_trivially_copyable<P>
     >::type
     use_memcpy;
 
@@ -668,7 +694,7 @@ void assign(I pos, V const& v)
     typedef typename
     ::boost::mpl::and_<
         is_corresponding_value<I, V>,
-        ::boost::has_trivial_assign<V>
+        is_trivially_copyable<V>
     >::type
     use_memcpy;
 
diff --git a/test/algorithms/buffer/buffer_multi_polygon.cpp b/test/algorithms/buffer/buffer_multi_polygon.cpp
index cd9ee11d64..db0780610e 100644
--- a/test/algorithms/buffer/buffer_multi_polygon.cpp
+++ b/test/algorithms/buffer/buffer_multi_polygon.cpp
@@ -281,6 +281,19 @@ static std::string const rt_u12
 static std::string const rt_u13
     = "MULTIPOLYGON(((6 4,6 5,7 5,6 4)),((3 2,3 3,4 3,3 2)),((7 8,7 9,8 9,8 8,7 8)),((4 9,4 10,5 10,4 9)),((7 7,7 8,8 7,7 7)),((2 6,2 7,3 7,2 6)),((0 1,1 2,1 1,0 1)),((3 1,4 2,4 1,3 1)),((2 5,2 6,3 6,2 5)),((3 5,4 4,3 4,2 4,3 5)),((4 1,5 2,5 1,4 1)),((2 0,2 1,3 1,2 0)),((5 7,5 8,6 7,5 7)),((0 2,0 3,1 3,0 2)),((9 8,9 9,10 9,10 8,9 8)),((7 5,7 6,8 5,7 5)),((5 6,5 7,6 6,5 6)),((0 6,0 7,1 7,1 6,0 6)),((5 0,5 1,6 1,5 0)),((8 7,8 8,9 8,8 7)),((4.5 4.5,5 4,4 4,4 5,5 5,4.5 4.5)),((6 2,5 2,5 3,6 3,7 3,8 2,7 2,6 2)),((8 6,8 7,9 7,9 6,9 5,8 5,8 6)),((8 1,9 0,8 0,7 0,8 1)))";
 
+static std::string const rt_v1
+    = "MULTIPOLYGON(((2 8,2 9,3 9,2 8)),((3 6,3 7,4 7,3 6)),((9 9,9 10,10 10,10 9,9 9)),((0 7,1 8,1 7,0 7)))";
+
+static std::string const rt_v2
+    = "MULTIPOLYGON(((8 4,8 5,9 5,9 4,8 4)),((2 5,2 6,3 5,2 5)),((9 7,9 8,10 8,10 7,9 7)),((2 2,2 3,3 2,2 2)),((6 6,7 5,6 5,6 6)),((6 6,6 7,7 7,7 6,6 6)),((8 9,9 9,8 8,8 9)),((8 9,7 9,8 10,8 9)))";
+
+static std::string const rt_v3
+    = "MULTIPOLYGON(((7 0,7 1,8 1,8 0,7 0)),((6 2,6 3,7 2,6 2)),((9 3,8 3,8 4,9 3)),((9 3,9 4,10 3,9 3)))";
+
+static std::string const rt_v4
+    = "MULTIPOLYGON(((5 4,5 5,6 5,6 4,5 4)),((7 1,6 1,7 2,7 1)),((7 1,8 1,8 0,7 0,7 1)),((6 1,5 1,5 2,6 1)))";
+
+
 static std::string const neighbouring
     = "MULTIPOLYGON(((0 0,0 10,10 10,10 0,0 0)),((10 10,10 20,20 20,20 10,10 10)))";
 
@@ -502,6 +515,11 @@ void test_all()
     test_one<multi_polygon_type, polygon_type>("rt_u12", rt_u12, join_miter, end_flat, 142.1348, 1.0);
     test_one<multi_polygon_type, polygon_type>("rt_u13", rt_u13, join_miter, end_flat, 115.4853, 1.0);
 
+    test_one<multi_polygon_type, polygon_type>("rt_v1", rt_v1, join_round32, end_flat, 26.9994, 1.0);
+    test_one<multi_polygon_type, polygon_type>("rt_v2", rt_v2, join_round32, end_flat, 47.3510, 1.0);
+    test_one<multi_polygon_type, polygon_type>("rt_v3", rt_v3, join_round32, end_flat, 22.9158, 1.0);
+    test_one<multi_polygon_type, polygon_type>("rt_v4", rt_v4, join_round32, end_flat, 23.4146, 1.0);
+
     test_one<multi_polygon_type, polygon_type>("neighbouring_small",
         neighbouring,
         join_round32, end_round32, 128.0, -1.0);
diff --git a/test/algorithms/buffer/test_buffer_svg.hpp b/test/algorithms/buffer/test_buffer_svg.hpp
index b4d20df566..87739e68e5 100644
--- a/test/algorithms/buffer/test_buffer_svg.hpp
+++ b/test/algorithms/buffer/test_buffer_svg.hpp
@@ -157,9 +157,8 @@ private :
                     << "/" << it->operations[1].enriched.get_next_turn_index()
                     //<< " frac " << it->operations[0].fraction
 
-    //                If you want to see (robust)point coordinates (e.g. to find duplicates)
+    //                If you want to see point coordinates (e.g. to find duplicates)
                     << std::endl << std::setprecision(16) << bg::wkt(it->point)
-                    << std::endl  << bg::wkt(it->robust_point)
 
                     << std::endl;
                 out << " " << bg::method_char(it->method)
diff --git a/test/algorithms/set_operations/difference/Jamfile b/test/algorithms/set_operations/difference/Jamfile
index f107240217..c719158c8c 100644
--- a/test/algorithms/set_operations/difference/Jamfile
+++ b/test/algorithms/set_operations/difference/Jamfile
@@ -18,15 +18,16 @@ test-suite boost-geometry-algorithms-difference
     :
     [ run difference.cpp                      : : : <define>BOOST_GEOMETRY_TEST_ONLY_ONE_TYPE
                                                   : algorithms_difference ]
+    [ run difference_areal_linear.cpp       : : : : algorithms_difference_areal_linear ]
     [ run difference_l_a_sph.cpp            : : : : algorithms_difference_l_a_sph ]
     [ run difference_linear_linear.cpp      : : : : algorithms_difference_linear_linear ]
-    [ run difference_areal_linear.cpp       : : : : algorithms_difference_areal_linear ]
-    [ run difference_pl_a.cpp               : : : : algorithms_difference_pl_a ]
-    [ run difference_pl_l.cpp               : : : : algorithms_difference_pl_l ]
-    [ run difference_pl_pl.cpp              : : : : algorithms_difference_pl_pl ]
     [ run difference_multi.cpp                : : : <define>BOOST_GEOMETRY_TEST_ONLY_ONE_TYPE
                                                   : algorithms_difference_multi ]
     [ run difference_multi_areal_linear.cpp : : : : algorithms_difference_multi_areal_linear ]
     [ run difference_multi_spike.cpp          : : : <define>BOOST_GEOMETRY_TEST_ONLY_ONE_TYPE
                                                   : algorithms_difference_multi_spike ]
+    [ run difference_pl_a.cpp               : : : : algorithms_difference_pl_a ]
+    [ run difference_pl_l.cpp               : : : : algorithms_difference_pl_l ]
+    [ run difference_pl_pl.cpp              : : : : algorithms_difference_pl_pl ]
+    [ run difference_tupled.cpp             : : : : algorithms_difference_tupled ]
     ;
diff --git a/test/algorithms/set_operations/difference/difference_tupled.cpp b/test/algorithms/set_operations/difference/difference_tupled.cpp
new file mode 100644
index 0000000000..cf7afb1ccf
--- /dev/null
+++ b/test/algorithms/set_operations/difference/difference_tupled.cpp
@@ -0,0 +1,407 @@
+// Boost.Geometry (aka GGL, Generic Geometry Library)
+
+// Copyright (c) 2020, Oracle and/or its affiliates.
+
+// Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle
+
+// Licensed under the Boost Software License version 1.0.
+// http://www.boost.org/users/license.html
+
+#include <geometry_test_common.hpp>
+
+#include <boost/geometry/algorithms/correct.hpp>
+#include <boost/geometry/algorithms/equals.hpp>
+#include <boost/geometry/algorithms/difference.hpp>
+#include <boost/geometry/geometries/geometries.hpp>
+#include <boost/geometry/io/wkt/wkt.hpp>
+#include <boost/geometry/strategies/cartesian/area.hpp>
+#include <boost/geometry/strategies/cartesian/intersection.hpp>
+#include <boost/geometry/strategies/cartesian/point_in_poly_winding.hpp>
+#include <boost/geometry/strategies/cartesian/point_in_point.hpp>
+
+#include <boost/tuple/tuple.hpp>
+
+typedef bg::model::point<double, 2, bg::cs::cartesian> Pt;
+typedef bg::model::linestring<Pt> Ls;
+typedef bg::model::polygon<Pt> Po;
+typedef bg::model::ring<Pt> R;
+typedef bg::model::multi_point<Pt> MPt;
+typedef bg::model::multi_linestring<Ls> MLs;
+typedef bg::model::multi_polygon<Po> MPo;
+
+#ifdef BOOST_GEOMETRY_CXX11_TUPLE
+
+#include <tuple>
+
+#endif
+
+template <typename G>
+inline void check(std::string const& wkt1,
+                  std::string const& wkt2,
+                  G const& g,
+                  std::string const& expected)
+{
+    G expect;
+    bg::read_wkt(expected, expect);
+    bg::correct(expect);
+    if (! boost::empty(g) || ! boost::empty(expect))
+    {
+        BOOST_CHECK_MESSAGE(
+            // Commented out becasue the output in reversed case may be slightly different
+            //   e.g. different number of duplicated points in MultiPoint
+            //boost::size(g) == boost::size(expect) &&
+            bg::equals(g, expect),
+            wkt1 << " \\ " << wkt2 << " -> " << bg::wkt(g)
+                 << " different than expected: " << expected
+        );
+    }
+}
+
+template <int I>
+inline void check(std::string const& wkt1,
+                  std::string const& wkt2,
+                  boost::tuple<MPt, MLs, MPo> const& tup,
+                  std::string const& out_str)
+{
+    check(wkt1, wkt2, boost::get<I>(tup), out_str);
+}
+
+template <int I>
+inline void check(std::string const& wkt1,
+                  std::string const& wkt2,
+                  std::pair<MPt, MLs> const& pair,
+                  std::string const& out_str)
+{
+    if (BOOST_GEOMETRY_CONDITION(I == 0))
+        check(wkt1, wkt2, pair.first, out_str);
+    else
+        check(wkt1, wkt2, pair.second, out_str);
+}
+
+#ifdef BOOST_GEOMETRY_CXX11_TUPLE
+
+template <int I>
+inline void check(std::string const& wkt1,
+                  std::string const& wkt2,
+                  std::tuple<MPt, MLs, MPo> const& tup,
+                  std::string const& out_str)
+{
+    check(wkt1, wkt2, std::get<I>(tup), out_str);
+}
+
+#endif
+
+template <typename Geometry>
+struct out_id
+    : boost::mpl::if_c
+        <
+            boost::is_base_of<bg::pointlike_tag, typename bg::tag<Geometry>::type>::value,
+            boost::mpl::int_<0>,
+            typename boost::mpl::if_c
+                <
+                    boost::is_base_of<bg::linear_tag, typename bg::tag<Geometry>::type>::value,
+                    boost::mpl::int_<1>,
+                    boost::mpl::int_<2>
+                >::type
+        >::type
+{};
+
+template <typename In1, typename In2, typename Tup>
+inline void test_one(std::string const& in1_str,
+                     std::string const& in2_str,
+                     std::string const& out1_str,
+                     std::string const& out2_str)
+{
+    In1 in1;
+    bg::read_wkt(in1_str, in1);
+    bg::correct(in1);
+
+    In2 in2;
+    bg::read_wkt(in2_str, in2);
+    bg::correct(in2);
+
+    {
+        Tup result;
+        bg::difference(in1, in2, result);
+        check<out_id<In1>::value>(in1_str, in2_str, result, out1_str);
+    }
+    {
+        Tup result;
+        bg::difference(in2, in1, result);
+        check<out_id<In2>::value>(in2_str, in1_str, result, out2_str);
+    }
+}
+
+template <typename Tup>
+inline void test_pp()
+{
+    test_one<Pt, Pt, Tup>(
+        "POINT(0 0)",
+        "POINT(0 0)",
+        "MULTIPOINT()",
+        "MULTIPOINT()");
+
+    test_one<Pt, Pt, Tup>(
+        "POINT(0 0)",
+        "POINT(1 1)",
+        "MULTIPOINT(0 0)",
+        "MULTIPOINT(1 1)");
+
+    test_one<Pt, MPt, Tup>(
+        "POINT(0 0)",
+        "MULTIPOINT(0 0, 1 1)",
+        "MULTIPOINT()",
+        "MULTIPOINT(1 1)");
+
+    test_one<Pt, MPt, Tup>(
+        "POINT(2 2)",
+        "MULTIPOINT(0 0, 1 1)",
+        "MULTIPOINT(2 2)",
+        "MULTIPOINT(0 0, 1 1)");
+
+    test_one<MPt, MPt, Tup>(
+        "MULTIPOINT(0 0, 1 1, 2 2)",
+        "MULTIPOINT(1 1, 3 3, 4 4)",
+        "MULTIPOINT(0 0, 2 2)",
+        "MULTIPOINT(3 3, 4 4)");
+}
+
+template <typename Tup>
+inline void test_pl()
+{
+    test_one<Pt, Ls, Tup>(
+        "POINT(0 0)",
+        "LINESTRING(0 0, 1 1)",
+        "MULTIPOINT()",
+        "MULTILINESTRING((0 0, 1 1))");
+
+    test_one<Pt, MLs, Tup>(
+        "POINT(0 1)",
+        "MULTILINESTRING((0 0, 1 1),(1 1, 2 2),(4 4, 5 5))",
+        "MULTIPOINT(0 1)",
+        "MULTILINESTRING((0 0, 1 1),(1 1, 2 2),(4 4, 5 5))");
+
+    test_one<MPt, Ls, Tup>(
+        "MULTIPOINT(0 0, 1 1, 2 2, 3 3)",
+        "LINESTRING(0 0, 1 1)",
+        "MULTIPOINT(2 2, 3 3)",
+        "MULTILINESTRING((0 0, 1 1))");
+
+    test_one<MPt, MLs, Tup>(
+        "MULTIPOINT(0 0, 1 1, 2 2, 3 3)",
+        "MULTILINESTRING((0 0, 1 1),(1 1, 2 2),(4 4, 5 5))",
+        "MULTIPOINT(3 3)",
+        "MULTILINESTRING((0 0, 1 1),(1 1, 2 2),(4 4, 5 5))");
+}
+
+template <typename Tup>
+inline void test_pa()
+{
+    test_one<Pt, R, Tup>(
+        "POINT(0 0)",
+        "POLYGON((0 0, 0 5, 5 5, 5 0, 0 0))",
+        "MULTIPOINT()",
+        "MULTIPOLYGON(((0 0, 0 5, 5 5, 5 0, 0 0)))");
+
+    test_one<Pt, Po, Tup>(
+        "POINT(0 0)",
+        "POLYGON((0 0, 0 5, 5 5, 5 0, 0 0),(0 0, 4 1, 4 4, 1 4, 0 0))",
+        "MULTIPOINT()",
+        "MULTIPOLYGON(((0 0, 0 5, 5 5, 5 0, 0 0),(0 0, 4 1, 4 4, 1 4, 0 0)))");
+
+    test_one<Pt, Po, Tup>(
+        "POINT(3 3)",
+        "POLYGON((0 0, 0 5, 5 5, 5 0, 0 0),(0 0, 4 1, 4 4, 1 4, 0 0))",
+        "MULTIPOINT(3 3)",
+        "MULTIPOLYGON(((0 0, 0 5, 5 5, 5 0, 0 0),(0 0, 4 1, 4 4, 1 4, 0 0)))");
+
+    test_one<Pt, MPo, Tup>(
+        "POINT(2 2)",
+        "MULTIPOLYGON(((0 0, 0 5, 5 5, 5 0, 0 0),(0 0, 4 1, 4 4, 1 4, 0 0)),((0 0, 1 2, 2 2, 2 1, 0 0)))",
+        "MULTIPOINT()",
+        "MULTIPOLYGON(((0 0, 0 5, 5 5, 5 0, 0 0),(0 0, 4 1, 4 4, 1 4, 0 0)),((0 0, 1 2, 2 2, 2 1, 0 0)))");
+
+    test_one<Pt, MPo, Tup>(
+        "POINT(3 3)",
+        "MULTIPOLYGON(((0 0, 0 5, 5 5, 5 0, 0 0),(0 0, 4 1, 4 4, 1 4, 0 0)),((0 0, 1 2, 2 2, 2 1, 0 0)))",
+        "MULTIPOINT(3 3)",
+        "MULTIPOLYGON(((0 0, 0 5, 5 5, 5 0, 0 0),(0 0, 4 1, 4 4, 1 4, 0 0)),((0 0, 1 2, 2 2, 2 1, 0 0)))");
+
+    test_one<Pt, MPo, Tup>(
+        "POINT(6 6)",
+        "MULTIPOLYGON(((0 0, 0 5, 5 5, 5 0, 0 0),(0 0, 4 1, 4 4, 1 4, 0 0)),((0 0, 1 2, 2 2, 2 1, 0 0)))",
+        "MULTIPOINT(6 6)",
+        "MULTIPOLYGON(((0 0, 0 5, 5 5, 5 0, 0 0),(0 0, 4 1, 4 4, 1 4, 0 0)),((0 0, 1 2, 2 2, 2 1, 0 0)))");
+
+    test_one<MPt, R, Tup>(
+        "MULTIPOINT(0 0, 1 1, 2 2, 3 3, 4 4, 5 5, 6 6)",
+        "POLYGON((0 0, 0 5, 5 5, 5 0, 0 0))",
+        "MULTIPOINT(6 6)",
+        "MULTIPOLYGON(((0 0, 0 5, 5 5, 5 0, 0 0)))");
+
+    test_one<MPt, Po, Tup>(
+        "MULTIPOINT(0 0, 1 1, 2 2, 3 3, 4 4, 5 5, 6 6)",
+        "POLYGON((0 0, 0 5, 5 5, 5 0, 0 0),(0 0, 4 1, 4 4, 1 4, 0 0))",
+        "MULTIPOINT(1 1, 2 2, 3 3, 6 6)",
+        "MULTIPOLYGON(((0 0, 0 5, 5 5, 5 0, 0 0),(0 0, 4 1, 4 4, 1 4, 0 0)))");
+
+    test_one<MPt, MPo, Tup>(
+        "MULTIPOINT(0 0, 1 1, 2 2, 3 3, 4 4, 5 5, 6 6)",
+        "MULTIPOLYGON(((0 0, 0 5, 5 5, 5 0, 0 0),(0 0, 4 1, 4 4, 1 4, 0 0)),((0 0, 1 2, 2 2, 2 1, 0 0)))",
+        "MULTIPOINT(3 3, 6 6)",
+        "MULTIPOLYGON(((0 0, 0 5, 5 5, 5 0, 0 0),(0 0, 4 1, 4 4, 1 4, 0 0)),((0 0, 1 2, 2 2, 2 1, 0 0)))");
+}
+
+template <typename Tup>
+inline void test_ll()
+{
+    test_one<Ls, Ls, Tup>(
+        "LINESTRING(0 0, 1 0, 2 1, 3 0)",
+        "LINESTRING(0 0, 1 0, 3 0, 4 0)",
+        "MULTILINESTRING((1 0, 2 1, 3 0))",
+        "MULTILINESTRING((1 0, 3 0, 4 0))");
+
+    test_one<Ls, MLs, Tup>(
+        "LINESTRING(0 0, 1 0, 2 1, 3 0)",
+        "MULTILINESTRING((0 0, 1 0, 3 0),(2 1, 2 2))",
+        "MULTILINESTRING((1 0, 2 1, 3 0))",
+        "MULTILINESTRING((1 0, 3 0),(2 1, 2 2))");
+
+    test_one<MLs, MLs, Tup>(
+        "MULTILINESTRING((0 0, 1 0, 2 1),(2 1, 3 0))",
+        "MULTILINESTRING((0 0, 1 0, 3 0),(2 1, 2 2))",
+        "MULTILINESTRING((1 0, 2 1),(2 1, 3 0))",
+        "MULTILINESTRING((1 0, 3 0),(2 1, 2 2))");
+
+    test_one<Ls, Ls, Tup>(
+        "LINESTRING(0 0, 0 5, 5 5, 5 0, 0 0)",
+        "LINESTRING(0 0, 0 1, 6 1, 5 2, 5 5, 5 6, 4 5, 4 7, 7 7, 7 0, 0 0)",
+        "MULTILINESTRING((0 1, 0 5, 5 5),(5 2, 5 0))",
+        "MULTILINESTRING((0 1, 6 1, 5 2),(5 5, 5 6, 4 5, 4 7, 7 7, 7 0, 5 0))");
+
+    test_one<MLs, MLs, Tup>(
+        "MULTILINESTRING((0 0, 0 5, 5 5, 5 0, 0 0),(0 0, 4 1, 4 4, 1 4, 0 0))",
+        "MULTILINESTRING((0 0, 1 4, 5 4, 5 1, 4 1, 0 0),(0 0, 2 1, 2 2, 1 2, 0 0))",
+        "MULTILINESTRING((0 0, 0 5, 5 5, 5 4),(5 1, 5 0, 0 0),(4 4, 4 1))",
+        "MULTILINESTRING((4 4, 5 4),(5 1, 4 1),(0 0, 2 1, 2 2, 1 2, 0 0))");
+}
+
+template <typename Tup>
+inline void test_la()
+{
+    test_one<Ls, R, Tup>(
+        "LINESTRING(0 2, -4 1, 0 0, 5 0, 9 1, 5 2, 9 3, 5 5, 4 9, 4 5, 3 3, 2 5, 2 9, 0 5)",
+        "POLYGON((0 0, 0 5, 5 5, 5 0, 0 0))",
+        "MULTILINESTRING((0 2, -4 1, 0 0),(5 0, 9 1, 5 2, 9 3, 5 5, 4 9, 4 5),(2 5, 2 9, 0 5))",
+        "MULTIPOLYGON(((0 0, 0 5, 5 5, 5 0, 0 0)))");
+
+    test_one<Ls, Po, Tup>(
+        "LINESTRING(1 4, -4 1, 0 0, 5 0, 9 1, 5 2, 9 3, 5 5, 4 9, 4 5, 3 3, 2 5, 2 9, 0 5)",
+        "POLYGON((0 0, 0 5, 5 5, 5 0, 0 0),(0 0, 4 1, 4 4, 1 4, 0 0))",
+        "MULTILINESTRING((0 3.4, -4 1, 0 0),(5 0, 9 1, 5 2, 9 3, 5 5, 4 9, 4 5),(3.5 4, 3 3, 2.5 4),(2 5, 2 9, 0 5))",
+        "MULTIPOLYGON(((0 0, 0 5, 5 5, 5 0, 0 0),(0 0, 4 1, 4 4, 1 4, 0 0)))");
+
+    test_one<MLs, R, Tup>(
+        "MULTILINESTRING((0 2, -4 1, 0 0, 5 0, 9 1, 5 2, 9 3, 5 5, 4 9), (4 9, 4 5, 3 3, 2 5, 2 9, 0 5))",
+        "POLYGON((0 0, 0 5, 5 5, 5 0, 0 0))",
+        "MULTILINESTRING((0 2, -4 1, 0 0),(5 0, 9 1, 5 2, 9 3, 5 5, 4 9),(4 9, 4 5),(2 5, 2 9, 0 5))",
+        "MULTIPOLYGON(((0 0, 0 5, 5 5, 5 0, 0 0)))");
+
+    test_one<MLs, Po, Tup>(
+        "MULTILINESTRING((1 4, -4 1, 0 0, 5 0, 9 1, 5 2, 9 3, 5 5, 4 9), (4 9, 4 5, 3 3, 2 5, 2 9, 0 5))",
+        "POLYGON((0 0, 0 5, 5 5, 5 0, 0 0),(0 0, 4 1, 4 4, 1 4, 0 0))",
+        "MULTILINESTRING((0 3.4, -4 1, 0 0),(5 0, 9 1, 5 2, 9 3, 5 5, 4 9),(4 9, 4 5),(3.5 4, 3 3, 2.5 4),(2 5, 2 9, 0 5))",
+        "MULTIPOLYGON(((0 0, 0 5, 5 5, 5 0, 0 0),(0 0, 4 1, 4 4, 1 4, 0 0)))");
+
+    test_one<MLs, MPo, Tup>(
+        "MULTILINESTRING((1 4, -4 1, 0 0, 5 0, 9 1, 5 2, 9 3, 5 5, 4 9), (4 9, 4 5, 4 4, 2 2, 2 5, 1 9, 0 5))",
+        "MULTIPOLYGON(((0 0, 0 5, 5 5, 5 0, 0 0),(0 0, 4 1, 4 4, 1 4, 0 0)),((0 0, 1 2, 2 2, 2 1, 0 0)))",
+        "MULTILINESTRING((0 3.4, -4 1, 0 0),(5 0, 9 1, 5 2, 9 3, 5 5, 4 9),(4 9, 4 5),(4 4, 2 2, 2 4),(2 5, 1 9, 0 5))",
+        "MULTIPOLYGON(((0 0, 0 5, 5 5, 5 0, 0 0),(0 0, 4 1, 4 4, 1 4, 0 0)),((0 0, 1 2, 2 2, 2 1, 0 0)))");
+
+    test_one<Ls, R, Tup>(
+        "LINESTRING(0 0, 0 5, 5 5, 5 0, 0 0)",
+        "POLYGON((0 0, 0 1, 6 1, 5 2, 5 5, 5 6, 4 5, 4 7, 7 7, 7 0, 0 0))",
+        "MULTILINESTRING((0 1, 0 5, 5 5),(5 2, 5 1))",
+        "MULTIPOLYGON(((0 0, 0 1, 6 1, 5 2, 5 5, 5 6, 4 5, 4 7, 7 7, 7 0, 0 0)))");
+
+    test_one<MLs, Po, Tup>(
+        "MULTILINESTRING((0 0, 0 5, 5 5, 5 0, 0 0),(0 0, 4 1, 4 4, 1 4, 0 0))",
+        "POLYGON((0 0, 1 4, 5 4, 5 1, 4 1, 0 0),(0 0, 2 1, 2 2, 1 2, 0 0))",
+        "MULTILINESTRING((0 0, 0 5, 5 5, 5 4),(5 1, 5 0, 0 0))",
+        "MULTIPOLYGON(((0 0, 1 4, 5 4, 5 1, 4 1, 0 0),(0 0, 2 1, 2 2, 1 2, 0 0)))");
+}
+
+template <typename Tup>
+inline void test_aa()
+{
+    test_one<R, R, Tup>(
+        "POLYGON((0 0, 0 5, 5 5, 5 0, 0 0))",
+        "POLYGON((0 0, 0 1, 6 1, 5 2, 5 5, 5 6, 4 5, 4 7, 7 7, 7 0, 0 0))",
+        "MULTIPOLYGON(((0 1,0 5,5 5,5 1,0 1)))",
+        "MULTIPOLYGON(((5 1,6 1,5 2,5 5,5 6,4 5,4 7,7 7,7 0,5 0,5 1)))");
+
+    test_one<R, MPo, Tup>(
+        "POLYGON((0 0, 0 5, 5 5, 5 0, 0 0))",
+        "MULTIPOLYGON(((0 0, 0 1, 6 1, 6 0, 0 0)),"
+                     "((6 1, 5 2, 5 5, 5 6, 4 5, 4 7, 7 7, 7 1, 6 1)))",
+        "MULTIPOLYGON(((0 1,0 5,5 5,5 1,0 1)))",
+        "MULTIPOLYGON(((5 1,6 1,6 0,5 0,5 1)),((5 2,5 5,5 6,4 5,4 7,7 7,7 1,6 1,5 2)))");
+
+    test_one<Po, Po, Tup>(
+        "POLYGON((0 0, 0 5, 5 5, 5 0, 0 0),(0 0, 4 1, 4 4, 1 4, 0 0))",
+        "POLYGON((0 0, 1 4, 5 4, 5 1, 4 1, 0 0),(0 0, 2 1, 2 2, 1 2, 0 0))",
+        "MULTIPOLYGON(((5 1,5 0,0 0,4 1,5 1)),((5 4,1 4,0 0,0 5,5 5,5 4)))",
+        "MULTIPOLYGON(((1 4,4 4,4 1,0 0,1 4),(0 0,2 1,2 2,1 2,0 0)))");
+
+    test_one<Po, MPo, Tup>(
+        "POLYGON((0 0, 0 5, 5 5, 5 0, 0 0),(0 0, 4 1, 4 4, 1 4, 0 0))",
+        "MULTIPOLYGON(((0 0, 1 4, 5 4, 5 1, 4 1, 0 0),(0 0, 2 1, 2 2, 1 2, 0 0)),"
+                     "((5 0, 5 1, 6 1, 6 4, 5 4, 3 6, 2 5, 2 7, 7 7, 7 0 5 0)))",
+        "MULTIPOLYGON(((4 5,5 4,1 4,0 0,0 5,4 5)),((5 1,5 0,0 0,4 1,5 1)))",
+        "MULTIPOLYGON(((5 0,5 1,6 1,6 4,5 4,5 5,4 5,3 6,2 5,2 7,7 7,7 0,5 0)),"
+                     "((1 4,4 4,4 1,0 0,1 4),(0 0,2 1,2 2,1 2,0 0)))");
+
+    test_one<MPo, MPo, Tup>(
+        "MULTIPOLYGON(((0 0, 0 5, 5 5, 5 0, 0 0),(0 0, 4 1, 4 4, 1 4, 0 0)),"
+                     "((2 6, 2 8, 8 8, 8 5, 7 5, 7 6, 2 6)))",
+        "MULTIPOLYGON(((0 0, 1 4, 5 4, 5 1, 4 1, 0 0),(0 0, 2 1, 2 2, 1 2, 0 0)),"
+                     "((5 0, 5 1, 6 1, 6 4, 5 4, 3 6, 2 5, 2 7, 7 7, 7 0 5 0)))",
+        "MULTIPOLYGON(((4 5,5 4,1 4,0 0,0 5,4 5)),"
+                     "((5 1,5 0,0 0,4 1,5 1)),"
+                     "((2 7,2 8,8 8,8 5,7 5,7 6,7 7,2 7)))",
+        "MULTIPOLYGON(((1 4,4 4,4 1,0 0,1 4),(0 0,2 1,2 2,1 2,0 0)),"
+                     "((5 1,6 1,6 4,5 4,5 5,4 5,3 6,7 6,7 5,7 0,5 0,5 1)),"
+                     "((3 6,2 5,2 6,3 6)))");
+}
+
+template <typename Tup>
+inline void test_pair()
+{
+    test_pp<Tup>();
+    test_pl<Tup>();
+    test_ll<Tup>();
+}
+
+template <typename Tup>
+inline void test_tuple()
+{
+    test_pp<Tup>();
+    test_pl<Tup>();
+    test_pa<Tup>();
+    test_ll<Tup>();
+    test_la<Tup>();
+    test_aa<Tup>();
+}
+
+int test_main(int, char* [])
+{
+    test_pair<std::pair<MPt, MLs> >();
+    test_tuple<boost::tuple<MPt, MLs, MPo> >();
+
+#ifdef BOOST_GEOMETRY_CXX11_TUPLE
+    test_tuple<std::tuple<MPt, MLs, MPo> >();
+#endif
+
+    return 0;
+}
diff --git a/test/algorithms/set_operations/sym_difference/Jamfile b/test/algorithms/set_operations/sym_difference/Jamfile
index 72bbeeab2b..7f3707ba1c 100644
--- a/test/algorithms/set_operations/sym_difference/Jamfile
+++ b/test/algorithms/set_operations/sym_difference/Jamfile
@@ -4,8 +4,8 @@
 # Copyright (c) 2008-2015 Bruno Lalande, Paris, France.
 # Copyright (c) 2009-2015 Mateusz Loskot, London, UK.
 #
-# This file was modified by Oracle on 2014, 2015.
-# Modifications copyright (c) 2014-2015, Oracle and/or its affiliates.
+# This file was modified by Oracle on 2014, 2015, 2020.
+# Modifications copyright (c) 2014-2020, Oracle and/or its affiliates.
 #
 # Contributed and/or modified by Menelaos Karavelas, on behalf of Oracle
 # Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle
@@ -16,6 +16,7 @@
 
 test-suite boost-geometry-algorithms-sym_difference
     :
-    [ run sym_difference_areal_areal.cpp : : : : algorithms_sym_difference_areal_areal ]
+    [ run sym_difference_areal_areal.cpp   : : : : algorithms_sym_difference_areal_areal ]
     [ run sym_difference_linear_linear.cpp : : : : algorithms_sym_difference_linear_linear ]
+	[ run sym_difference_tupled.cpp        : : : : algorithms_sym_difference_tupled ]
     ;
diff --git a/test/algorithms/set_operations/sym_difference/sym_difference_tupled.cpp b/test/algorithms/set_operations/sym_difference/sym_difference_tupled.cpp
new file mode 100644
index 0000000000..14864e1d6b
--- /dev/null
+++ b/test/algorithms/set_operations/sym_difference/sym_difference_tupled.cpp
@@ -0,0 +1,383 @@
+// Boost.Geometry (aka GGL, Generic Geometry Library)
+
+// Copyright (c) 2020, Oracle and/or its affiliates.
+
+// Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle
+
+// Licensed under the Boost Software License version 1.0.
+// http://www.boost.org/users/license.html
+
+#include <geometry_test_common.hpp>
+
+#include <boost/geometry/algorithms/correct.hpp>
+#include <boost/geometry/algorithms/equals.hpp>
+#include <boost/geometry/algorithms/sym_difference.hpp>
+#include <boost/geometry/geometries/geometries.hpp>
+#include <boost/geometry/io/wkt/wkt.hpp>
+#include <boost/geometry/strategies/cartesian/area.hpp>
+#include <boost/geometry/strategies/cartesian/intersection.hpp>
+#include <boost/geometry/strategies/cartesian/point_in_poly_winding.hpp>
+#include <boost/geometry/strategies/cartesian/point_in_point.hpp>
+
+#include <boost/tuple/tuple.hpp>
+
+typedef bg::model::point<double, 2, bg::cs::cartesian> Pt;
+typedef bg::model::linestring<Pt> Ls;
+typedef bg::model::polygon<Pt> Po;
+typedef bg::model::ring<Pt> R;
+typedef bg::model::multi_point<Pt> MPt;
+typedef bg::model::multi_linestring<Ls> MLs;
+typedef bg::model::multi_polygon<Po> MPo;
+
+#ifdef BOOST_GEOMETRY_CXX11_TUPLE
+
+#include <tuple>
+
+#endif
+
+template <typename G>
+inline void check(std::string const& wkt1,
+                  std::string const& wkt2,
+                  G const& g,
+                  std::string const& expected)
+{
+    G expect;
+    bg::read_wkt(expected, expect);
+    bg::correct(expect);
+    if (! boost::empty(g) || ! boost::empty(expect))
+    {
+        BOOST_CHECK_MESSAGE(
+            // Commented out becasue the output in reversed case may be slightly different
+            //   e.g. different number of duplicated points in MultiPoint
+            //boost::size(g) == boost::size(expect) &&
+            bg::equals(g, expect),
+            wkt1 << " ^ " << wkt2 << " -> " << bg::wkt(g)
+                 << " different than expected: " << expected
+        );
+    }
+}
+
+inline void check(std::string const& wkt1,
+                  std::string const& wkt2,
+                  boost::tuple<MPt, MLs, MPo> const& tup,
+                  std::string const& out_p_str,
+                  std::string const& out_l_str,
+                  std::string const& out_a_str)
+{
+    check(wkt1, wkt2, boost::get<0>(tup), out_p_str);
+    check(wkt1, wkt2, boost::get<1>(tup), out_l_str);
+    check(wkt1, wkt2, boost::get<2>(tup), out_a_str);
+}
+
+inline void check(std::string const& wkt1,
+                  std::string const& wkt2,
+                  std::pair<MPt, MLs> const& pair,
+                  std::string const& out_p_str,
+                  std::string const& out_l_str,
+                  std::string const& )
+{
+    check(wkt1, wkt2, pair.first, out_p_str);
+    check(wkt1, wkt2, pair.second, out_l_str);
+}
+
+#ifdef BOOST_GEOMETRY_CXX11_TUPLE
+
+inline void check(std::string const& wkt1,
+                  std::string const& wkt2,
+                  std::tuple<MPt, MLs, MPo> const& tup,
+                  std::string const& out_p_str,
+                  std::string const& out_l_str,
+                  std::string const& out_a_str)
+{
+    check(wkt1, wkt2, std::get<0>(tup), out_p_str);
+    check(wkt1, wkt2, std::get<1>(tup), out_l_str);
+    check(wkt1, wkt2, std::get<2>(tup), out_a_str);
+}
+
+#endif
+
+template <typename In1, typename In2, typename Tup>
+inline void test_one(std::string const& in1_str,
+                     std::string const& in2_str,
+                     std::string const& out_p_str = "MULTIPOINT()",
+                     std::string const& out_l_str = "MULTILINESTRING()",
+                     std::string const& out_a_str = "MULTIPOLYGON()")
+{
+    In1 in1;
+    bg::read_wkt(in1_str, in1);
+    bg::correct(in1);
+
+    In2 in2;
+    bg::read_wkt(in2_str, in2);
+    bg::correct(in2);
+
+    {
+        Tup result;
+        bg::sym_difference(in1, in2, result);
+        check(in1_str, in2_str, result, out_p_str, out_l_str, out_a_str);
+    }
+    {
+        Tup result;
+        bg::sym_difference(in2, in1, result);
+        check(in1_str, in2_str, result, out_p_str, out_l_str, out_a_str);
+    }
+}
+
+template <typename Tup>
+inline void test_pp()
+{
+    test_one<Pt, Pt, Tup>(
+        "POINT(0 0)",
+        "POINT(0 0)",
+        "MULTIPOINT()");
+
+    test_one<Pt, MPt, Tup>(
+        "POINT(0 0)",
+        "MULTIPOINT(0 0, 1 1)",
+        "MULTIPOINT(1 1)");
+
+    test_one<MPt, MPt, Tup>(
+        "MULTIPOINT(0 0, 1 1, 2 2)",
+        "MULTIPOINT(1 1, 2 2, 3 3)",
+        "MULTIPOINT(0 0, 3 3)");
+}
+
+template <typename Tup>
+inline void test_pl()
+{
+    test_one<Pt, Ls, Tup>(
+        "POINT(0 0)",
+        "LINESTRING(0 0, 1 1)",
+        "MULTIPOINT()",
+        "MULTILINESTRING((0 0, 1 1))");
+
+    test_one<Pt, MLs, Tup>(
+        "POINT(1 1)",
+        "MULTILINESTRING((0 0, 1 1),(1 1, 2 2),(4 4, 5 5))",
+        "MULTIPOINT()",
+        "MULTILINESTRING((0 0, 1 1),(1 1, 2 2),(4 4, 5 5))");
+
+    test_one<MPt, Ls, Tup>(
+        "MULTIPOINT(0 0, 1 1, 2 2, 3 3)",
+        "LINESTRING(0 0, 1 1)",
+        "MULTIPOINT(2 2, 3 3)",
+        "MULTILINESTRING((0 0, 1 1))");
+
+    test_one<MPt, MLs, Tup>(
+        "MULTIPOINT(0 0, 1 1, 2 2, 3 3)",
+        "MULTILINESTRING((0 0, 1 1),(1 1, 2 2),(4 4, 5 5))",
+        "MULTIPOINT(3 3)",
+        "MULTILINESTRING((0 0, 1 1),(1 1, 2 2),(4 4, 5 5))");
+}
+
+template <typename Tup>
+inline void test_pa()
+{
+    test_one<Pt, R, Tup>(
+        "POINT(0 0)",
+        "POLYGON((0 0, 0 5, 5 5, 5 0, 0 0))",
+        "MULTIPOINT()",
+        "MULTILINESTRING()",
+        "MULTIPOLYGON(((0 0, 0 5, 5 5, 5 0, 0 0)))");
+
+    test_one<Pt, Po, Tup>(
+        "POINT(0 0)",
+        "POLYGON((0 0, 0 5, 5 5, 5 0, 0 0),(0 0, 4 1, 4 4, 1 4, 0 0))",
+        "MULTIPOINT()",
+        "MULTILINESTRING()",
+        "MULTIPOLYGON(((0 0, 0 5, 5 5, 5 0, 0 0),(0 0, 4 1, 4 4, 1 4, 0 0)))");
+
+    test_one<Pt, MPo, Tup>(
+        "POINT(0 0)",
+        "MULTIPOLYGON(((0 0, 0 5, 5 5, 5 0, 0 0),(0 0, 4 1, 4 4, 1 4, 0 0)),((0 0, 1 2, 2 2, 2 1, 0 0)))",
+        "MULTIPOINT()",
+        "MULTILINESTRING()",
+        "MULTIPOLYGON(((0 0, 0 5, 5 5, 5 0, 0 0),(0 0, 4 1, 4 4, 1 4, 0 0)),((0 0, 1 2, 2 2, 2 1, 0 0)))");
+
+    test_one<MPt, R, Tup>(
+        "MULTIPOINT(0 0, 1 1, 2 2, 3 3, 4 4, 5 5, 6 6)",
+        "POLYGON((0 0, 0 5, 5 5, 5 0, 0 0))",
+        "MULTIPOINT(6 6)",
+        "MULTILINESTRING()",
+        "MULTIPOLYGON(((0 0, 0 5, 5 5, 5 0, 0 0)))");
+
+    test_one<MPt, Po, Tup>(
+        "MULTIPOINT(0 0, 1 1, 2 2, 3 3, 4 4, 5 5, 6 6)",
+        "POLYGON((0 0, 0 5, 5 5, 5 0, 0 0),(0 0, 4 1, 4 4, 1 4, 0 0))",
+        "MULTIPOINT(1 1, 2 2, 3 3, 6 6)",
+        "MULTILINESTRING()",
+        "MULTIPOLYGON(((0 0, 0 5, 5 5, 5 0, 0 0),(0 0, 4 1, 4 4, 1 4, 0 0)))");
+
+    test_one<MPt, MPo, Tup>(
+        "MULTIPOINT(0 0, 0 0, 1 1, 2 2, 3 3, 4 4, 5 5, 6 6)",
+        "MULTIPOLYGON(((0 0, 0 5, 5 5, 5 0, 0 0),(0 0, 4 1, 4 4, 1 4, 0 0)),((0 0, 1 2, 2 2, 2 1, 0 0)))",
+        "MULTIPOINT(3 3, 6 6)",
+        "MULTILINESTRING()",
+        "MULTIPOLYGON(((0 0, 0 5, 5 5, 5 0, 0 0),(0 0, 4 1, 4 4, 1 4, 0 0)),((0 0, 1 2, 2 2, 2 1, 0 0)))");
+}
+
+template <typename Tup>
+inline void test_ll()
+{
+    test_one<Ls, Ls, Tup>(
+        "LINESTRING(0 0, 1 0, 2 1, 3 0)",
+        "LINESTRING(0 0, 1 0, 3 0, 4 0)",
+        "MULTIPOINT()",
+        "MULTILINESTRING((1 0, 2 1, 3 0),(1 0, 3 0, 4 0))");
+
+    test_one<Ls, MLs, Tup>(
+        "LINESTRING(0 0, 1 0, 2 1, 3 0)",
+        "MULTILINESTRING((0 0, 1 0, 3 0),(2 1, 2 2))",
+        "MULTIPOINT()",
+        "MULTILINESTRING((1 0, 2 1, 3 0),(1 0, 3 0),(2 1, 2 2))");
+
+    test_one<MLs, MLs, Tup>(
+        "MULTILINESTRING((0 0, 1 0, 2 1),(2 1, 3 0))",
+        "MULTILINESTRING((0 0, 1 0, 3 0),(2 1, 2 2))",
+        "MULTIPOINT()",
+        "MULTILINESTRING((1 0, 2 1),(2 1, 3 0),(1 0, 3 0),(2 1, 2 2))");
+
+    test_one<Ls, Ls, Tup>(
+        "LINESTRING(0 0, 0 5, 5 5, 5 0, 0 0)",
+        "LINESTRING(0 0, 0 1, 6 1, 5 2, 5 5, 5 6, 4 5, 4 7, 7 7, 7 0, 0 0)",
+        "MULTIPOINT()",
+        "MULTILINESTRING((0 1, 0 5, 5 5),(5 2, 5 0),(0 1, 6 1, 5 2), (5 5, 5 6, 4 5, 4 7, 7 7, 7 0, 5 0))");
+
+    test_one<MLs, MLs, Tup>(
+        "MULTILINESTRING((0 0, 0 5, 5 5, 5 0, 0 0),(0 0, 4 1, 4 4, 1 4, 0 0))",
+        "MULTILINESTRING((0 0, 1 4, 5 4, 5 1, 4 1, 0 0),(0 0, 2 1, 2 2, 1 2, 0 0))",
+        "MULTIPOINT()",
+        "MULTILINESTRING((0 0, 0 5, 5 5, 5 4),(5 1, 5 0, 0 0),(4 1, 4 4),(4 4, 5 4),(5 1, 4 1),(0 0, 2 1, 2 2, 1 2, 0 0))");
+}
+
+template <typename Tup>
+inline void test_la()
+{
+    test_one<Ls, R, Tup>(
+        "LINESTRING(0 2, -4 1, 0 0, 5 0, 9 1, 5 2, 9 3, 5 5, 4 9, 4 5, 3 3, 2 5, 2 9, 0 5)",
+        "POLYGON((0 0, 0 5, 5 5, 5 0, 0 0))",
+        "MULTIPOINT()",
+        "MULTILINESTRING((0 2, -4 1, 0 0),(5 0, 9 1, 5 2, 9 3, 5 5, 4 9, 4 5),(2 5, 2 9, 0 5))",
+        "MULTIPOLYGON(((0 0, 0 5, 5 5, 5 0, 0 0)))");
+
+    test_one<Ls, Po, Tup>(
+        "LINESTRING(1 4, -4 1, 0 0, 5 0, 9 1, 5 2, 9 3, 5 5, 4 9, 4 5, 3 3, 2 5, 2 9, 0 5)",
+        "POLYGON((0 0, 0 5, 5 5, 5 0, 0 0),(0 0, 4 1, 4 4, 1 4, 0 0))",
+        "MULTIPOINT()",
+        "MULTILINESTRING((0 3.4, -4 1, 0 0),(5 0, 9 1, 5 2, 9 3, 5 5, 4 9, 4 5),(3.5 4, 3 3, 2.5 4),(2 5, 2 9, 0 5))",
+        "MULTIPOLYGON(((0 0, 0 5, 5 5, 5 0, 0 0),(0 0, 4 1, 4 4, 1 4, 0 0)))");
+
+    test_one<MLs, R, Tup>(
+        "MULTILINESTRING((0 2, -4 1, 0 0, 5 0, 9 1, 5 2, 9 3, 5 5, 4 9), (4 9, 4 5, 3 3, 2 5, 2 9, 0 5))",
+        "POLYGON((0 0, 0 5, 5 5, 5 0, 0 0))",
+        "MULTIPOINT()",
+        "MULTILINESTRING((0 2, -4 1, 0 0),(5 0, 9 1, 5 2, 9 3, 5 5, 4 9),(4 9, 4 5),(2 5, 2 9, 0 5))",
+        "MULTIPOLYGON(((0 0, 0 5, 5 5, 5 0, 0 0)))");
+
+    test_one<MLs, Po, Tup>(
+        "MULTILINESTRING((1 4, -4 1, 0 0, 5 0, 9 1, 5 2, 9 3, 5 5, 4 9), (4 9, 4 5, 3 3, 2 5, 2 9, 0 5))",
+        "POLYGON((0 0, 0 5, 5 5, 5 0, 0 0),(0 0, 4 1, 4 4, 1 4, 0 0))",
+        "MULTIPOINT()",
+        "MULTILINESTRING((0 3.4, -4 1, 0 0),(5 0, 9 1, 5 2, 9 3, 5 5, 4 9),(4 9, 4 5),(3.5 4, 3 3, 2.5 4),(2 5, 2 9, 0 5))",
+        "MULTIPOLYGON(((0 0, 0 5, 5 5, 5 0, 0 0),(0 0, 4 1, 4 4, 1 4, 0 0)))");
+
+    test_one<MLs, MPo, Tup>(
+        "MULTILINESTRING((1 4, -4 1, 0 0, 5 0, 9 1, 5 2, 9 3, 5 5, 4 9), (4 9, 4 5, 4 4, 2 2, 2 5, 1 9, 0 5))",
+        "MULTIPOLYGON(((0 0, 0 5, 5 5, 5 0, 0 0),(0 0, 4 1, 4 4, 1 4, 0 0)),((0 0, 1 2, 2 2, 2 1, 0 0)))",
+        "MULTIPOINT()",
+        "MULTILINESTRING((0 3.4, -4 1, 0 0),(5 0, 9 1, 5 2, 9 3, 5 5, 4 9),(4 9, 4 5),(4 4, 2 2, 2 4),(2 5, 1 9, 0 5))",
+        "MULTIPOLYGON(((0 0, 0 5, 5 5, 5 0, 0 0),(0 0, 4 1, 4 4, 1 4, 0 0)),((0 0, 1 2, 2 2, 2 1, 0 0)))");
+
+    test_one<Ls, R, Tup>(
+        "LINESTRING(0 0, 0 5, 5 5, 5 0, 0 0)",
+        "POLYGON((0 0, 0 1, 6 1, 5 2, 5 5, 5 6, 4 5, 4 7, 7 7, 7 0, 0 0))",
+        "MULTIPOINT()",
+        "MULTILINESTRING((0 1, 0 5, 5 5),(5 2, 5 1))",
+        "MULTIPOLYGON(((0 0, 0 1, 6 1, 5 2, 5 5, 5 6, 4 5, 4 7, 7 7, 7 0, 0 0)))");
+
+    test_one<MLs, Po, Tup>(
+        "MULTILINESTRING((0 0, 0 5, 5 5, 5 0, 0 0),(0 0, 4 1, 4 4, 1 4, 0 0))",
+        "POLYGON((0 0, 1 4, 5 4, 5 1, 4 1, 0 0),(0 0, 2 1, 2 2, 1 2, 0 0))",
+        "MULTIPOINT()",
+        "MULTILINESTRING((0 0, 0 5, 5 5, 5 4),(5 1, 5 0, 0 0))",
+        "MULTIPOLYGON(((0 0, 1 4, 5 4, 5 1, 4 1, 0 0),(0 0, 2 1, 2 2, 1 2, 0 0)))");
+}
+
+template <typename Tup>
+inline void test_aa()
+{
+    test_one<R, R, Tup>(
+        "POLYGON((0 0, 0 5, 5 5, 5 0, 0 0))",
+        "POLYGON((0 0, 0 1, 6 1, 5 2, 5 5, 5 6, 4 5, 4 7, 7 7, 7 0, 0 0))",
+        "MULTIPOINT()",
+        "MULTILINESTRING()",
+        "MULTIPOLYGON(((4 5,4 7,7 7,7 0,5 0,5 1,0 1,0 5,4 5),(5 5,5 6,4 5,5 5),(5 2,5 1,6 1,5 2)))");
+
+    test_one<R, MPo, Tup>(
+        "POLYGON((0 0, 0 5, 5 5, 5 0, 0 0))",
+        "MULTIPOLYGON(((0 0, 0 1, 6 1, 6 0, 0 0)),"
+                     "((6 1, 5 2, 5 5, 5 6, 4 5, 4 7, 7 7, 7 1, 6 1)))",
+        "MULTIPOINT()",
+        "MULTILINESTRING()",
+        "MULTIPOLYGON(((4 5,4 7,7 7,7 1,6 1,5 2,5 1,0 1,0 5,4 5),(5 5,5 6,4 5,5 5)),"
+                     "((5 1,6 1,6 0,5 0,5 1)))");
+
+    test_one<Po, Po, Tup>(
+        "POLYGON((0 0, 0 5, 5 5, 5 0, 0 0),(0 0, 4 1, 4 4, 1 4, 0 0))",
+        "POLYGON((0 0, 1 4, 5 4, 5 1, 4 1, 0 0),(0 0, 2 1, 2 2, 1 2, 0 0))",
+        "MULTIPOINT()",
+        "MULTILINESTRING()",
+        "MULTIPOLYGON(((4 1,5 1,5 0,0 0,0 5,5 5,5 4,4 4,4 1),(0 0,2 1,2 2,1 2,0 0)))");
+
+    test_one<Po, MPo, Tup>(
+        "POLYGON((0 0, 0 5, 5 5, 5 0, 0 0),(0 0, 4 1, 4 4, 1 4, 0 0))",
+        "MULTIPOLYGON(((0 0, 1 4, 5 4, 5 1, 4 1, 0 0),(0 0, 2 1, 2 2, 1 2, 0 0)),"
+                     "((5 0, 5 1, 6 1, 6 4, 5 4, 3 6, 2 5, 2 7, 7 7, 7 0 5 0)))",
+        "MULTIPOINT()",
+        "MULTILINESTRING()",
+        "MULTIPOLYGON(((0 0,0 5,2 5,2 7,7 7,7 0,0 0),"
+                        "(4 4,4 1,5 1,6 1,6 4,4 4),(5 4,5 5,4 5,5 4),"
+                        "(0 0,2 1,2 2,1 2,0 0),(4 5,3 6,2 5,4 5)))");
+
+    test_one<MPo, MPo, Tup>(
+        "MULTIPOLYGON(((0 0, 0 5, 5 5, 5 0, 0 0),(0 0, 4 1, 4 4, 1 4, 0 0)),"
+                     "((2 6, 2 8, 8 8, 8 5, 7 5, 7 6, 2 6)))",
+        "MULTIPOLYGON(((0 0, 1 4, 5 4, 5 1, 4 1, 0 0),(0 0, 2 1, 2 2, 1 2, 0 0)),"
+                     "((5 0, 5 1, 6 1, 6 4, 5 4, 3 6, 2 5, 2 7, 7 7, 7 0 5 0)))",
+        "MULTIPOINT()",
+        "MULTILINESTRING()",
+        "MULTIPOLYGON(((0 0,0 5,4 5,3 6,7 6,7 7,2 7,2 8,8 8,8 5,7 5,7 0,0 0),"
+                        "(4 4,4 1,5 1,6 1,6 4,4 4),(5 4,5 5,4 5,5 4),(0 0,2 1,2 2,1 2,0 0)),"
+                     "((2 5,2 6,3 6,2 5)))");
+}
+
+template <typename Tup>
+inline void test_pair()
+{
+    test_pp<Tup>();
+    test_pl<Tup>();
+    test_ll<Tup>();
+}
+
+template <typename Tup>
+inline void test_tuple()
+{
+    test_pp<Tup>();
+    test_pl<Tup>();
+    test_pa<Tup>();
+    test_ll<Tup>();
+    test_la<Tup>();
+    test_aa<Tup>();
+}
+
+int test_main(int, char* [])
+{
+    test_pair<std::pair<MPt, MLs> >();
+    test_tuple<boost::tuple<MPt, MLs, MPo> >();
+
+#ifdef BOOST_GEOMETRY_CXX11_TUPLE
+    test_tuple<std::tuple<MPt, MLs, MPo> >();
+#endif
+
+    return 0;
+}
diff --git a/test/algorithms/set_operations/union/Jamfile b/test/algorithms/set_operations/union/Jamfile
index e351a870df..8d184ffa98 100644
--- a/test/algorithms/set_operations/union/Jamfile
+++ b/test/algorithms/set_operations/union/Jamfile
@@ -4,8 +4,8 @@
 # Copyright (c) 2008-2015 Bruno Lalande, Paris, France.
 # Copyright (c) 2009-2015 Mateusz Loskot, London, UK.
 #
-# This file was modified by Oracle on 2014, 2015, 2017, 2018.
-# Modifications copyright (c) 2014-2018, Oracle and/or its affiliates.
+# This file was modified by Oracle on 2014, 2015, 2017, 2018, 2020.
+# Modifications copyright (c) 2014-2020, Oracle and/or its affiliates.
 #
 # Contributed and/or modified by Menelaos Karavelas, on behalf of Oracle
 # Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle
@@ -24,4 +24,5 @@ test-suite boost-geometry-algorithms-union
     [ run union_multi.cpp         : : : <define>BOOST_GEOMETRY_TEST_ONLY_ONE_TYPE
                                         : algorithms_union_multi ]
     [ run union_pl_pl.cpp         : : : : algorithms_union_pl_pl ]
+	[ run union_tupled.cpp        : : : : algorithms_union_tupled ]
     ;
diff --git a/test/algorithms/set_operations/union/union_tupled.cpp b/test/algorithms/set_operations/union/union_tupled.cpp
new file mode 100644
index 0000000000..5ca9741efc
--- /dev/null
+++ b/test/algorithms/set_operations/union/union_tupled.cpp
@@ -0,0 +1,388 @@
+// Boost.Geometry (aka GGL, Generic Geometry Library)
+
+// Copyright (c) 2020, Oracle and/or its affiliates.
+
+// Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle
+
+// Licensed under the Boost Software License version 1.0.
+// http://www.boost.org/users/license.html
+
+#include <geometry_test_common.hpp>
+
+#include <boost/geometry/algorithms/correct.hpp>
+#include <boost/geometry/algorithms/equals.hpp>
+#include <boost/geometry/algorithms/union.hpp>
+#include <boost/geometry/geometries/geometries.hpp>
+#include <boost/geometry/io/wkt/wkt.hpp>
+#include <boost/geometry/strategies/cartesian/area.hpp>
+#include <boost/geometry/strategies/cartesian/intersection.hpp>
+#include <boost/geometry/strategies/cartesian/point_in_poly_winding.hpp>
+#include <boost/geometry/strategies/cartesian/point_in_point.hpp>
+
+#include <boost/tuple/tuple.hpp>
+
+typedef bg::model::point<double, 2, bg::cs::cartesian> Pt;
+typedef bg::model::linestring<Pt> Ls;
+typedef bg::model::polygon<Pt> Po;
+typedef bg::model::ring<Pt> R;
+typedef bg::model::multi_point<Pt> MPt;
+typedef bg::model::multi_linestring<Ls> MLs;
+typedef bg::model::multi_polygon<Po> MPo;
+
+#ifdef BOOST_GEOMETRY_CXX11_TUPLE
+
+#include <tuple>
+
+#endif
+
+template <typename G>
+inline void check(std::string const& wkt1,
+                  std::string const& wkt2,
+                  G const& g,
+                  std::string const& expected)
+{
+    G expect;
+    bg::read_wkt(expected, expect);
+    bg::correct(expect);
+    if (! boost::empty(g) || ! boost::empty(expect))
+    {
+        BOOST_CHECK_MESSAGE(
+            // Commented out becasue the output in reversed case may be slightly different
+            //   e.g. different number of duplicated points in MultiPoint
+            //boost::size(g) == boost::size(expect) &&
+            bg::equals(g, expect),
+            wkt1 << " + " << wkt2 << " -> " << bg::wkt(g)
+                 << " different than expected: " << expected
+        );
+    }
+}
+
+inline void check(std::string const& wkt1,
+                  std::string const& wkt2,
+                  boost::tuple<MPt, MLs, MPo> const& tup,
+                  std::string const& out_p_str,
+                  std::string const& out_l_str,
+                  std::string const& out_a_str)
+{
+    check(wkt1, wkt2, boost::get<0>(tup), out_p_str);
+    check(wkt1, wkt2, boost::get<1>(tup), out_l_str);
+    check(wkt1, wkt2, boost::get<2>(tup), out_a_str);
+}
+
+inline void check(std::string const& wkt1,
+                  std::string const& wkt2,
+                  std::pair<MPt, MLs> const& pair,
+                  std::string const& out_p_str,
+                  std::string const& out_l_str,
+                  std::string const& )
+{
+    check(wkt1, wkt2, pair.first, out_p_str);
+    check(wkt1, wkt2, pair.second, out_l_str);
+}
+
+#ifdef BOOST_GEOMETRY_CXX11_TUPLE
+
+inline void check(std::string const& wkt1,
+                  std::string const& wkt2,
+                  std::tuple<MPt, MLs, MPo> const& tup,
+                  std::string const& out_p_str,
+                  std::string const& out_l_str,
+                  std::string const& out_a_str)
+{
+    check(wkt1, wkt2, std::get<0>(tup), out_p_str);
+    check(wkt1, wkt2, std::get<1>(tup), out_l_str);
+    check(wkt1, wkt2, std::get<2>(tup), out_a_str);
+}
+
+#endif
+
+template <typename In1, typename In2, typename Tup>
+inline void test_one(std::string const& in1_str,
+                     std::string const& in2_str,
+                     std::string const& out_p_str = "MULTIPOINT()",
+                     std::string const& out_l_str = "MULTILINESTRING()",
+                     std::string const& out_a_str = "MULTIPOLYGON()")
+{
+    In1 in1;
+    bg::read_wkt(in1_str, in1);
+    bg::correct(in1);
+
+    In2 in2;
+    bg::read_wkt(in2_str, in2);
+    bg::correct(in2);
+
+    {
+        Tup result;
+        bg::union_(in1, in2, result);
+        check(in1_str, in2_str, result, out_p_str, out_l_str, out_a_str);
+    }
+    {
+        Tup result;
+        bg::union_(in2, in1, result);
+        check(in1_str, in2_str, result, out_p_str, out_l_str, out_a_str);
+    }
+}
+
+template <typename Tup>
+inline void test_pp()
+{
+    test_one<Pt, Pt, Tup>(
+        "POINT(0 0)",
+        "POINT(0 0)",
+        "MULTIPOINT(0 0)");
+
+    test_one<Pt, MPt, Tup>(
+        "POINT(0 0)",
+        "MULTIPOINT(0 0, 1 1)",
+        "MULTIPOINT(0 0, 1 1)");
+
+    test_one<MPt, MPt, Tup>(
+        "MULTIPOINT(0 0, 1 1, 2 2)",
+        "MULTIPOINT(1 1, 2 2, 3 3)",
+        "MULTIPOINT(0 0, 1 1, 2 2, 3 3)");
+}
+
+template <typename Tup>
+inline void test_pl()
+{
+    test_one<Pt, Ls, Tup>(
+        "POINT(0 0)",
+        "LINESTRING(0 0, 1 1)",
+        "MULTIPOINT()",
+        "MULTILINESTRING((0 0, 1 1))");
+
+    test_one<Pt, Ls, Tup>(
+        "POINT(2 2)",
+        "LINESTRING(0 0, 1 1)",
+        "MULTIPOINT(2 2)",
+        "MULTILINESTRING((0 0, 1 1))");
+
+    test_one<Pt, MLs, Tup>(
+        "POINT(1 1)",
+        "MULTILINESTRING((0 0, 1 1),(1 1, 2 2),(4 4, 5 5))",
+        "MULTIPOINT()",
+        "MULTILINESTRING((0 0, 1 1),(1 1, 2 2),(4 4, 5 5))");
+
+    test_one<MPt, Ls, Tup>(
+        "MULTIPOINT(0 0, 1 1, 2 2, 3 3)",
+        "LINESTRING(0 0, 1 1)",
+        "MULTIPOINT(2 2, 3 3)",
+        "MULTILINESTRING((0 0, 1 1))");
+
+    test_one<MPt, MLs, Tup>(
+        "MULTIPOINT(0 0, 1 1, 2 2, 3 3)",
+        "MULTILINESTRING((0 0, 1 1),(1 1, 2 2),(4 4, 5 5))",
+        "MULTIPOINT(3 3)",
+        "MULTILINESTRING((0 0, 1 1),(1 1, 2 2),(4 4, 5 5))");
+}
+
+template <typename Tup>
+inline void test_pa()
+{
+    test_one<Pt, R, Tup>(
+        "POINT(0 0)",
+        "POLYGON((0 0, 0 5, 5 5, 5 0, 0 0))",
+        "MULTIPOINT()",
+        "MULTILINESTRING()",
+        "MULTIPOLYGON(((0 0, 0 5, 5 5, 5 0, 0 0)))");
+
+    test_one<Pt, Po, Tup>(
+        "POINT(2 2)",
+        "POLYGON((0 0, 0 5, 5 5, 5 0, 0 0),(0 0, 4 1, 4 4, 1 4, 0 0))",
+        "MULTIPOINT(2 2)",
+        "MULTILINESTRING()",
+        "MULTIPOLYGON(((0 0, 0 5, 5 5, 5 0, 0 0),(0 0, 4 1, 4 4, 1 4, 0 0)))");
+
+    test_one<Pt, MPo, Tup>(
+        "POINT(2 2)",
+        "MULTIPOLYGON(((0 0, 0 5, 5 5, 5 0, 0 0),(0 0, 4 1, 4 4, 1 4, 0 0)),((0 0, 1 2, 2 2, 2 1, 0 0)))",
+        "MULTIPOINT()",
+        "MULTILINESTRING()",
+        "MULTIPOLYGON(((0 0, 0 5, 5 5, 5 0, 0 0),(0 0, 4 1, 4 4, 1 4, 0 0)),((0 0, 1 2, 2 2, 2 1, 0 0)))");
+
+    test_one<MPt, R, Tup>(
+        "MULTIPOINT(0 0, 1 1, 2 2, 3 3, 4 4, 5 5, 6 6)",
+        "POLYGON((0 0, 0 5, 5 5, 5 0, 0 0))",
+        "MULTIPOINT(6 6)",
+        "MULTILINESTRING()",
+        "MULTIPOLYGON(((0 0, 0 5, 5 5, 5 0, 0 0)))");
+
+    test_one<MPt, Po, Tup>(
+        "MULTIPOINT(0 0, 1 1, 2 2, 3 3, 4 4, 5 5, 6 6)",
+        "POLYGON((0 0, 0 5, 5 5, 5 0, 0 0),(0 0, 4 1, 4 4, 1 4, 0 0))",
+        "MULTIPOINT(1 1, 2 2, 3 3, 6 6)",
+        "MULTILINESTRING()",
+        "MULTIPOLYGON(((0 0, 0 5, 5 5, 5 0, 0 0),(0 0, 4 1, 4 4, 1 4, 0 0)))");
+
+    test_one<MPt, MPo, Tup>(
+        "MULTIPOINT(0 0, 0 0, 1 1, 2 2, 3 3, 4 4, 5 5, 6 6)",
+        "MULTIPOLYGON(((0 0, 0 5, 5 5, 5 0, 0 0),(0 0, 4 1, 4 4, 1 4, 0 0)),((0 0, 1 2, 2 2, 2 1, 0 0)))",
+        "MULTIPOINT(3 3, 6 6)",
+        "MULTILINESTRING()",
+        "MULTIPOLYGON(((0 0, 0 5, 5 5, 5 0, 0 0),(0 0, 4 1, 4 4, 1 4, 0 0)),((0 0, 1 2, 2 2, 2 1, 0 0)))");
+}
+
+template <typename Tup>
+inline void test_ll()
+{
+    test_one<Ls, Ls, Tup>(
+        "LINESTRING(0 0, 1 0, 2 1, 3 0)",
+        "LINESTRING(0 0, 1 0, 3 0, 4 0)",
+        "MULTIPOINT()",
+        "MULTILINESTRING((0 0,1 0,2 1,3 0),(1 0,3 0,4 0))");
+
+    test_one<Ls, MLs, Tup>(
+        "LINESTRING(0 0, 1 0, 2 1, 3 0)",
+        "MULTILINESTRING((0 0, 1 0, 3 0),(2 1, 2 2))",
+        "MULTIPOINT()",
+        "MULTILINESTRING((0 0,1 0,2 1,3 0),(1 0,3 0),(2 1,2 2))");
+
+    test_one<MLs, MLs, Tup>(
+        "MULTILINESTRING((0 0, 1 0, 2 1),(2 1, 3 0))",
+        "MULTILINESTRING((0 0, 1 0, 3 0),(2 1, 2 2))",
+        "MULTIPOINT()",
+        "MULTILINESTRING((0 0,1 0,3 0),(2 1,2 2),(1 0,2 1),(2 1,3 0))");
+
+    test_one<Ls, Ls, Tup>(
+        "LINESTRING(0 0, 0 5, 5 5, 5 0, 0 0)",
+        "LINESTRING(0 0, 0 1, 6 1, 5 2, 5 5, 5 6, 4 5, 4 7, 7 7, 7 0, 0 0)",
+        "MULTIPOINT()",
+        "MULTILINESTRING((0 0,0 5,5 5,5 0,0 0),(0 1,6 1,5 2),"
+                        "(5 5,5 6,4 5,4 7,7 7,7 0,5 0))");
+
+    test_one<MLs, MLs, Tup>(
+        "MULTILINESTRING((0 0, 0 5, 5 5, 5 0, 0 0),(0 0, 4 1, 4 4, 1 4, 0 0))",
+        "MULTILINESTRING((0 0, 1 4, 5 4, 5 1, 4 1, 0 0),(0 0, 2 1, 2 2, 1 2, 0 0))",
+        "MULTIPOINT()",
+        "MULTILINESTRING((0 0,1 4,5 4,5 1,4 1,0 0),(0 0,2 1,2 2,1 2,0 0),"
+                        "(0 0,0 5,5 5,5 4),(5 1,5 0,0 0),(4 1,4 4))");
+}
+
+template <typename Tup>
+inline void test_la()
+{
+    test_one<Ls, R, Tup>(
+        "LINESTRING(0 2, -4 1, 0 0, 5 0, 9 1, 5 2, 9 3, 5 5, 4 9, 4 5, 3 3, 2 5, 2 9, 0 5)",
+        "POLYGON((0 0, 0 5, 5 5, 5 0, 0 0))",
+        "MULTIPOINT()",
+        "MULTILINESTRING((0 2, -4 1, 0 0),(5 0, 9 1, 5 2, 9 3, 5 5, 4 9, 4 5),(2 5, 2 9, 0 5))",
+        "MULTIPOLYGON(((0 0, 0 5, 5 5, 5 0, 0 0)))");
+
+    test_one<Ls, Po, Tup>(
+        "LINESTRING(1 4, -4 1, 0 0, 5 0, 9 1, 5 2, 9 3, 5 5, 4 9, 4 5, 3 3, 2 5, 2 9, 0 5)",
+        "POLYGON((0 0, 0 5, 5 5, 5 0, 0 0),(0 0, 4 1, 4 4, 1 4, 0 0))",
+        "MULTIPOINT()",
+        "MULTILINESTRING((0 3.4, -4 1, 0 0),(5 0, 9 1, 5 2, 9 3, 5 5, 4 9, 4 5),(3.5 4, 3 3, 2.5 4),(2 5, 2 9, 0 5))",
+        "MULTIPOLYGON(((0 0, 0 5, 5 5, 5 0, 0 0),(0 0, 4 1, 4 4, 1 4, 0 0)))");
+
+    test_one<MLs, R, Tup>(
+        "MULTILINESTRING((0 2, -4 1, 0 0, 5 0, 9 1, 5 2, 9 3, 5 5, 4 9), (4 9, 4 5, 3 3, 2 5, 2 9, 0 5))",
+        "POLYGON((0 0, 0 5, 5 5, 5 0, 0 0))",
+        "MULTIPOINT()",
+        "MULTILINESTRING((0 2, -4 1, 0 0),(5 0, 9 1, 5 2, 9 3, 5 5, 4 9),(4 9, 4 5),(2 5, 2 9, 0 5))",
+        "MULTIPOLYGON(((0 0, 0 5, 5 5, 5 0, 0 0)))");
+
+    test_one<MLs, Po, Tup>(
+        "MULTILINESTRING((1 4, -4 1, 0 0, 5 0, 9 1, 5 2, 9 3, 5 5, 4 9), (4 9, 4 5, 3 3, 2 5, 2 9, 0 5))",
+        "POLYGON((0 0, 0 5, 5 5, 5 0, 0 0),(0 0, 4 1, 4 4, 1 4, 0 0))",
+        "MULTIPOINT()",
+        "MULTILINESTRING((0 3.4, -4 1, 0 0),(5 0, 9 1, 5 2, 9 3, 5 5, 4 9),(4 9, 4 5),(3.5 4, 3 3, 2.5 4),(2 5, 2 9, 0 5))",
+        "MULTIPOLYGON(((0 0, 0 5, 5 5, 5 0, 0 0),(0 0, 4 1, 4 4, 1 4, 0 0)))");
+
+    test_one<MLs, MPo, Tup>(
+        "MULTILINESTRING((1 4, -4 1, 0 0, 5 0, 9 1, 5 2, 9 3, 5 5, 4 9), (4 9, 4 5, 4 4, 2 2, 2 5, 1 9, 0 5))",
+        "MULTIPOLYGON(((0 0, 0 5, 5 5, 5 0, 0 0),(0 0, 4 1, 4 4, 1 4, 0 0)),((0 0, 1 2, 2 2, 2 1, 0 0)))",
+        "MULTIPOINT()",
+        "MULTILINESTRING((0 3.4, -4 1, 0 0),(5 0, 9 1, 5 2, 9 3, 5 5, 4 9),(4 9, 4 5),(4 4, 2 2, 2 4),(2 5, 1 9, 0 5))",
+        "MULTIPOLYGON(((0 0, 0 5, 5 5, 5 0, 0 0),(0 0, 4 1, 4 4, 1 4, 0 0)),((0 0, 1 2, 2 2, 2 1, 0 0)))");
+
+    test_one<Ls, R, Tup>(
+        "LINESTRING(0 0, 0 5, 5 5, 5 0, 0 0)",
+        "POLYGON((0 0, 0 1, 6 1, 5 2, 5 5, 5 6, 4 5, 4 7, 7 7, 7 0, 0 0))",
+        "MULTIPOINT()",
+        "MULTILINESTRING((0 1, 0 5, 5 5),(5 2, 5 1))",
+        "MULTIPOLYGON(((0 0, 0 1, 6 1, 5 2, 5 5, 5 6, 4 5, 4 7, 7 7, 7 0, 0 0)))");
+
+    test_one<MLs, Po, Tup>(
+        "MULTILINESTRING((0 0, 0 5, 5 5, 5 0, 0 0),(0 0, 4 1, 4 4, 1 4, 0 0))",
+        "POLYGON((0 0, 1 4, 5 4, 5 1, 4 1, 0 0),(0 0, 2 1, 2 2, 1 2, 0 0))",
+        "MULTIPOINT()",
+        "MULTILINESTRING((0 0, 0 5, 5 5, 5 4),(5 1, 5 0, 0 0))",
+        "MULTIPOLYGON(((0 0, 1 4, 5 4, 5 1, 4 1, 0 0),(0 0, 2 1, 2 2, 1 2, 0 0)))");
+}
+
+template <typename Tup>
+inline void test_aa()
+{
+    test_one<R, R, Tup>(
+        "POLYGON((0 0, 0 5, 5 5, 5 0, 0 0))",
+        "POLYGON((0 0, 0 1, 6 1, 5 2, 5 5, 5 6, 4 5, 4 7, 7 7, 7 0, 0 0))",
+        "MULTIPOINT()",
+        "MULTILINESTRING()",
+        "MULTIPOLYGON(((0 5,4 5,4 7,7 7,7 0,0 0,0 5),(5 1,6 1,5 2,5 1),(5 5,5 6,4 5,5 5)))");
+
+    test_one<R, MPo, Tup>(
+        "POLYGON((0 0, 0 5, 5 5, 5 0, 0 0))",
+        "MULTIPOLYGON(((0 0, 0 1, 6 1, 6 0, 0 0)),"
+                     "((6 1, 5 2, 5 5, 5 6, 4 5, 4 7, 7 7, 7 1, 6 1)))",
+        "MULTIPOINT()",
+        "MULTILINESTRING()",
+        "MULTIPOLYGON(((0 5,4 5,4 7,7 7,7 1,6 1,6 0,0 0,0 5),(5 1,6 1,5 2,5 1),(5 5,5 6,4 5,5 5)))");
+
+    test_one<Po, Po, Tup>(
+        "POLYGON((0 0, 0 5, 5 5, 5 0, 0 0),(0 0, 4 1, 4 4, 1 4, 0 0))",
+        "POLYGON((0 0, 1 4, 5 4, 5 1, 4 1, 0 0),(0 0, 2 1, 2 2, 1 2, 0 0))",
+        "MULTIPOINT()",
+        "MULTILINESTRING()",
+        "MULTIPOLYGON(((5 0,0 0,0 5,5 5,5 0),(0 0,2 1,2 2,1 2,0 0)))");
+
+    test_one<Po, MPo, Tup>(
+        "POLYGON((0 0, 0 5, 5 5, 5 0, 0 0),(0 0, 4 1, 4 4, 1 4, 0 0))",
+        "MULTIPOLYGON(((0 0, 1 4, 5 4, 5 1, 4 1, 0 0),(0 0, 2 1, 2 2, 1 2, 0 0)),"
+                     "((5 0, 5 1, 6 1, 6 4, 5 4, 3 6, 2 5, 2 7, 7 7, 7 0 5 0)))",
+        "MULTIPOINT()",
+        "MULTILINESTRING()",
+        "MULTIPOLYGON(((2 5,2 7,7 7,7 0,0 0,0 5,2 5),(4 5,3 6,2 5,4 5),"
+                      "(5 4,5 1,6 1,6 4,5 4),(0 0,2 1,2 2,1 2,0 0)))");
+
+    test_one<MPo, MPo, Tup>(
+        "MULTIPOLYGON(((0 0, 0 5, 5 5, 5 0, 0 0),(0 0, 4 1, 4 4, 1 4, 0 0)),"
+                     "((2 6, 2 8, 8 8, 8 5, 7 5, 7 6, 2 6)))",
+        "MULTIPOLYGON(((0 0, 1 4, 5 4, 5 1, 4 1, 0 0),(0 0, 2 1, 2 2, 1 2, 0 0)),"
+                     "((5 0, 5 1, 6 1, 6 4, 5 4, 3 6, 2 5, 2 7, 7 7, 7 0 5 0)))",
+        "MULTIPOINT()",
+        "MULTILINESTRING()",
+        "MULTIPOLYGON(((0 0,0 5,2 5,2 7,2 8,8 8,8 5,7 5,7 0,0 0),"
+                      "(5 4,5 1,6 1,6 4,5 4),(0 0,2 1,2 2,1 2,0 0),(4 5,3 6,2 5,4 5)))");
+}
+
+template <typename Tup>
+inline void test_pair()
+{
+    test_pp<Tup>();
+    test_pl<Tup>();
+    test_ll<Tup>();
+}
+
+template <typename Tup>
+inline void test_tuple()
+{
+    test_pp<Tup>();
+    test_pl<Tup>();
+    test_pa<Tup>();
+    test_ll<Tup>();
+    test_la<Tup>();
+    test_aa<Tup>();
+}
+
+int test_main(int, char* [])
+{
+    test_pair<std::pair<MPt, MLs> >();
+    test_tuple<boost::tuple<MPt, MLs, MPo> >();
+
+#ifdef BOOST_GEOMETRY_CXX11_TUPLE
+    test_tuple<std::tuple<MPt, MLs, MPo> >();
+#endif
+
+    return 0;
+}
diff --git a/test/arithmetic/infinite_line_functions.cpp b/test/arithmetic/infinite_line_functions.cpp
index 9c879178cc..845513415b 100755
--- a/test/arithmetic/infinite_line_functions.cpp
+++ b/test/arithmetic/infinite_line_functions.cpp
@@ -96,77 +96,6 @@ void test_get_intersection()
     verify_point_on_line(q, bg::get<0>(ip), bg::get<1>(ip));
 }
 
-template <typename T>
-void test_same_direction()
-{
-    bg::model::infinite_line<T> p, q;
-
-    // Exactly opposite, diagonal
-    p = bg::detail::make::make_infinite_line<T>(2, 1, 12, 11);
-    q = bg::detail::make::make_infinite_line<T>(12, 11, 2, 1);
-    BOOST_CHECK(! bg::arithmetic::similar_direction(p, q));
-
-    // Exactly opposite, horizontal
-    p = bg::detail::make::make_infinite_line<T>(0, 0, 10, 0);
-    q = bg::detail::make::make_infinite_line<T>(10, 0, 0, 0);
-    BOOST_CHECK(! bg::arithmetic::similar_direction(p, q));
-
-    // Exactly opposite, vertical
-    p = bg::detail::make::make_infinite_line<T>(0, 0, 0, 10);
-    q = bg::detail::make::make_infinite_line<T>(0, 10, 0, 0);
-    BOOST_CHECK(! bg::arithmetic::similar_direction(p, q));
-
-    // Exactly equal, diagonal
-    p = bg::detail::make::make_infinite_line<T>(0, 0, 10, 10);
-    q = bg::detail::make::make_infinite_line<T>(0, 0, 10, 10);
-    BOOST_CHECK(bg::arithmetic::similar_direction(p, q));
-
-    // Exactly equal, horizontal
-    p = bg::detail::make::make_infinite_line<T>(0, 0, 10, 0);
-    q = bg::detail::make::make_infinite_line<T>(0, 0, 10, 0);
-    BOOST_CHECK(bg::arithmetic::similar_direction(p, q));
-
-    // Exactly equal, vertical
-    p = bg::detail::make::make_infinite_line<T>(0, 0, 0, 10);
-    q = bg::detail::make::make_infinite_line<T>(0, 0, 0, 10);
-    BOOST_CHECK(bg::arithmetic::similar_direction(p, q));
-
-    // Coming together, diagonal
-    p = bg::detail::make::make_infinite_line<T>(0, 0, 10, 10);
-    q = bg::detail::make::make_infinite_line<T>(20, 20, 10, 10);
-    BOOST_CHECK(! bg::arithmetic::similar_direction(p, q));
-
-    // Leaving from common point, diagonal
-    p = bg::detail::make::make_infinite_line<T>(10, 10, 0, 0);
-    q = bg::detail::make::make_infinite_line<T>(0, 0, 10, 10);
-    BOOST_CHECK(! bg::arithmetic::similar_direction(p, q));
-
-    // Continuing each other, diagonal
-    p = bg::detail::make::make_infinite_line<T>(0, 0, 10, 10);
-    q = bg::detail::make::make_infinite_line<T>(10, 10, 20, 20);
-    BOOST_CHECK(bg::arithmetic::similar_direction(p, q));
-
-    // (Nearly) perpendicular
-    p = bg::detail::make::make_infinite_line<T>(0, 0, 10, 10);
-    q = bg::detail::make::make_infinite_line<T>(0, 0, -10, 10);
-    BOOST_CHECK(! bg::arithmetic::similar_direction(p, q));
-
-    // 45 deg
-    p = bg::detail::make::make_infinite_line<T>(0, 0, 10, 10);
-    q = bg::detail::make::make_infinite_line<T>(0, 0, 0, 10);
-    BOOST_CHECK(bg::arithmetic::similar_direction(p, q));
-
-    // a bit more than 45 deg
-    p = bg::detail::make::make_infinite_line<T>(0, 0, 10, 10);
-    q = bg::detail::make::make_infinite_line<T>(0, 0, -1, 10);
-    BOOST_CHECK(! bg::arithmetic::similar_direction(p, q));
-
-    // 135 deg
-    p = bg::detail::make::make_infinite_line<T>(0, 0, 10, 10);
-    q = bg::detail::make::make_infinite_line<T>(0, 0, -10, 0);
-    BOOST_CHECK(! bg::arithmetic::similar_direction(p, q));
-}
-
 template <typename T>
 void test_degenerate()
 {
diff --git a/test/strategies/Jamfile b/test/strategies/Jamfile
index d47835650e..e8b8031cd5 100644
--- a/test/strategies/Jamfile
+++ b/test/strategies/Jamfile
@@ -35,7 +35,6 @@ test-suite boost-geometry-strategies
     [ run segment_intersection_collinear.cpp : : : : strategies_segment_intersection_collinear ]
     [ run segment_intersection_geo.cpp       : : : : strategies_segment_intersection_geo ]
     [ run segment_intersection_sph.cpp       : : : : strategies_segment_intersection_sph ]
-#    [ run side_of_intersection.cpp           : : : : strategies_side_of_intersection ]
     [ run spherical_side.cpp                 : : : : strategies_spherical_side ]
     [ run thomas.cpp                         : : : : strategies_thomas ]
     [ run transform_cs.cpp                   : : : : strategies_transform_cs ]
diff --git a/test/strategies/side_of_intersection.cpp b/test/strategies/side_of_intersection.cpp
deleted file mode 100644
index b42270e3d1..0000000000
--- a/test/strategies/side_of_intersection.cpp
+++ /dev/null
@@ -1,108 +0,0 @@
-// Boost.Geometry (aka GGL, Generic Geometry Library)
-// Unit Test
-
-// Copyright (c) 2011-2015 Barend Gehrels, Amsterdam, the Netherlands.
-
-// 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/strategies/cartesian/side_of_intersection.hpp>
-#include <boost/geometry/geometries/point_xy.hpp>
-#include <boost/geometry/geometries/segment.hpp>
-
-
-namespace bg = boost::geometry;
-
-int test_main(int, char* [])
-{
-    typedef bg::model::d2::point_xy<boost::long_long_type> point;
-    typedef bg::model::segment<point> segment;
-
-    typedef bg::strategy::side::side_of_intersection side;
-
-    point no_fb(-99, -99);
-
-    segment a(point(20, 10), point(10, 20));
-
-    segment b1(point(11, 16), point(20, 14));  // IP with a: (14.857, 15.143)
-    segment b2(point(10, 16), point(20, 14));  // IP with a: (15, 15)
-
-    segment c1(point(15, 16), point(13, 8));
-    segment c2(point(15, 16), point(14, 8));
-    segment c3(point(15, 16), point(15, 8));
-
-
-    BOOST_CHECK_EQUAL( 1, side::apply(a, b1, c1, no_fb));
-    BOOST_CHECK_EQUAL(-1, side::apply(a, b1, c2, no_fb));
-    BOOST_CHECK_EQUAL(-1, side::apply(a, b1, c3, no_fb));
-
-    BOOST_CHECK_EQUAL( 1, side::apply(a, b2, c1, no_fb));
-    BOOST_CHECK_EQUAL( 1, side::apply(a, b2, c2, no_fb));
-    BOOST_CHECK_EQUAL( 0, side::apply(a, b2, c3, no_fb));
-
-    // Collinear cases
-    // 1: segments intersecting are collinear (with a fallback point):
-    {
-        point fb(5, 5);
-        segment col1(point(0, 5), point(5, 5));
-        segment col2(point(5, 5), point(10, 5)); // One IP with col1 at (5,5)
-        segment col3(point(5, 0), point(5, 5));
-        BOOST_CHECK_EQUAL( 0, side::apply(col1, col2, col3, fb));
-    }
-    // 2: segment of side calculation collinear with one of the segments
-    {
-        point fb(5, 5);
-        segment col1(point(0, 5), point(10, 5));
-        segment col2(point(5, 5), point(5, 12)); // IP with col1 at (5,5)
-        segment col3(point(5, 1), point(5, 5)); // collinear with col2
-        BOOST_CHECK_EQUAL( 0, side::apply(col1, col2, col3, fb));
-    }
-    {
-        point fb(5, 5);
-        segment col1(point(10, 5), point(0, 5));
-        segment col2(point(5, 5), point(5, 12)); // IP with col1 at (5,5)
-        segment col3(point(5, 1), point(5, 5)); // collinear with col2
-        BOOST_CHECK_EQUAL( 0, side::apply(col1, col2, col3, fb));
-    }
-    {
-        point fb(5, 5);
-        segment col1(point(0, 5), point(10, 5));
-        segment col2(point(5, 12), point(5, 5)); // IP with col1 at (5,5)
-        segment col3(point(5, 1), point(5, 5)); // collinear with col2
-        BOOST_CHECK_EQUAL( 0, side::apply(col1, col2, col3, fb));
-    }
-
-    {
-        point fb(517248, -517236);
-        segment col1(point(-172408, -517236), point(862076, -517236));
-        segment col2(point(517248, -862064), point(517248, -172408));
-        segment col3(point(517248, -172408), point(517248, -517236));
-        BOOST_CHECK_EQUAL( 0, side::apply(col1, col2, col3, fb));
-    }
-    {
-        point fb(-221647, -830336);
-        segment col1(point(-153817, -837972), point(-222457, -830244));
-        segment col2(point(-221139, -833615), point(-290654, -384388));
-        segment col3(point(-255266, -814663), point(-264389, -811197));
-        BOOST_CHECK_EQUAL(1, side::apply(col1, col2, col3, fb)); // Left of segment...
-    }
-
-
-    {
-        point fb(27671131, 30809240);
-        segment col1(point(27671116, 30809247), point(27675474, 30807351));
-        segment col2(point(27666779, 30811130), point(27671139, 30809237));
-        segment col3(point(27671122, 30809244), point(27675480, 30807348));
-        BOOST_CHECK_EQUAL(1, side::apply(col1, col2, col3, fb)); // Left of segment...
-    }
-
-    // TODO: we might add a check calculating the IP, determining the side
-    // with the normal side strategy, and verify the results are equal
-
-    return 0;
-}
-