fix(rtree): clamp content differences to prevent FMA-induced numerical instability

Might fix #1452

Author: barendgehrels

include/boost/geometry/index/detail/algorithms/content.hpp

@@ -112,6 +112,19 @@ typename default_content_result<Indexable>::type content(Indexable const& b)
             >::apply(b);
 }
 
+// Returns the content increase when expanding 'original' to 'expanded'.
+// Precondition: 'original' is covered by 'expanded'. Under this precondition
+// the result is mathematically non-negative, so clamping to 0 is valid and
+// guards against floating-point contraction artifacts (e.g. GCC FMA fusion
+// across function boundaries) that can make the difference negative or
+// spuriously non-zero for equal boxes. See GitHub issue #1452.
+template <typename Box>
+typename default_content_result<Box>::type content_diff(Box const& expanded, Box const& original)
+{
+    using content_type = typename default_content_result<Box>::type;
+    return (std::max)(content_type(0), content(expanded) - content(original));
+}
+
 #if defined(BOOST_GCC)
 #pragma GCC pop_options
 #endif

include/boost/geometry/index/detail/rtree/linear/redistribute_elements.hpp

@@ -16,6 +16,7 @@
 #ifndef BOOST_GEOMETRY_INDEX_DETAIL_RTREE_LINEAR_REDISTRIBUTE_ELEMENTS_HPP
 #define BOOST_GEOMETRY_INDEX_DETAIL_RTREE_LINEAR_REDISTRIBUTE_ELEMENTS_HPP
 
+#include <algorithm>
 #include <type_traits>
 
 #include <boost/core/ignore_unused.hpp>
@@ -432,8 +433,10 @@ struct redistribute_elements<MembersHolder, linear_tag>
                         content_type enlarged_content1 = index::detail::content(enlarged_box1);
                         content_type enlarged_content2 = index::detail::content(enlarged_box2);
 
-                        content_type content_increase1 = enlarged_content1 - content1;
-                        content_type content_increase2 = enlarged_content2 - content2;
+                        content_type const content_increase1
+                            = (std::max)(content_type(0), enlarged_content1 - content1);
+                        content_type const content_increase2
+                            = (std::max)(content_type(0), enlarged_content2 - content2);
 
                         // choose group which box content have to be enlarged least or has smaller content or has fewer elements
                         if ( content_increase1 < content_increase2 ||

include/boost/geometry/index/detail/rtree/quadratic/redistribute_elements.hpp

@@ -74,9 +74,10 @@ inline void pick_seeds(Elements const& elements,
 
             bounded_indexable_view bounded_ind1(ind1, strategy);
             bounded_indexable_view bounded_ind2(ind2, strategy);
-            content_type free_content = ( index::detail::content(enlarged_box)
-                                            - index::detail::content(bounded_ind1) )
-                                                - index::detail::content(bounded_ind2);
+            content_type free_content = (std::max)(content_type(0),
+                ( index::detail::content(enlarged_box)
+                    - index::detail::content(bounded_ind1) )
+                        - index::detail::content(bounded_ind2));
 
             if ( greatest_free_content < free_content )
             {
@@ -273,7 +274,7 @@ struct redistribute_elements<MembersHolder, quadratic_tag>
         typedef typename boost::iterator_value<It>::type element_type;
         typedef typename rtree::element_indexable_type<element_type, translator_type>::type indexable_type;
 
-        content_type greatest_content_incrase_diff = 0;
+        content_type greatest_content_increase_diff = 0;
         It out_it = first;
         out_content_increase1 = 0;
         out_content_increase2 = 0;
@@ -291,18 +292,20 @@ struct redistribute_elements<MembersHolder, quadratic_tag>
             content_type enlarged_content1 = index::detail::content(enlarged_box1);
             content_type enlarged_content2 = index::detail::content(enlarged_box2);
 
-            content_type content_incrase1 = (enlarged_content1 - content1);
-            content_type content_incrase2 = (enlarged_content2 - content2);
+            content_type const content_increase1
+                = (std::max)(content_type(0), enlarged_content1 - content1);
+            content_type const content_increase2
+                = (std::max)(content_type(0), enlarged_content2 - content2);
 
-            content_type content_incrase_diff = content_incrase1 < content_incrase2 ?
-                content_incrase2 - content_incrase1 : content_incrase1 - content_incrase2;
+            content_type content_increase_diff = content_increase1 < content_increase2 ?
+                content_increase2 - content_increase1 : content_increase1 - content_increase2;
 
-            if ( greatest_content_incrase_diff < content_incrase_diff )
+            if ( greatest_content_increase_diff < content_increase_diff )
             {
-                greatest_content_incrase_diff = content_incrase_diff;
+                greatest_content_increase_diff = content_increase_diff;
                 out_it = el_it;
-                out_content_increase1 = content_incrase1;
-                out_content_increase2 = content_incrase2;
+                out_content_increase1 = content_increase1;
+                out_content_increase2 = content_increase2;
             }
         }
 

include/boost/geometry/index/detail/rtree/rstar/choose_next_node.hpp

@@ -117,7 +117,7 @@ class choose_next_node<MembersHolder, choose_by_overlap_diff_tag>
 
             // areas difference
             content_type content = index::detail::content(box_exp);
-            content_type content_diff = content - index::detail::content(ch_i.first);
+            content_type content_diff = index::detail::content_diff(box_exp, ch_i.first);
 
             children_contents[i].set(i, content, content_diff);
 
@@ -246,7 +246,7 @@ class choose_next_node<MembersHolder, choose_by_overlap_diff_tag>
 
             // areas difference
             content_type content = index::detail::content(box_exp);
-            content_type content_diff = content - index::detail::content(ch_i.first);
+            content_type content_diff = index::detail::content_diff(box_exp, ch_i.first);
 
             // update the result
             if ( content_diff < smallest_content_diff ||

include/boost/geometry/index/detail/rtree/visitors/insert.hpp

@@ -89,7 +89,7 @@ class choose_next_node<MembersHolder, choose_by_content_diff_tag>
 
             // areas difference
             content_type content = index::detail::content(box_exp);
-            content_type content_diff = content - index::detail::content(ch_i.first);
+            content_type content_diff = index::detail::content_diff(box_exp, ch_i.first);
 
             // update the result
             if ( content_diff < smallest_content_diff ||