feat: bit-identical distance/intersects/touches — strict equality verification

- Add geos_distance() helper to base.h (matching predicate pattern)
- Add distance bindings to BIND_GEOM_PREDS + DEF_LS_PRED/DEF_POLY_PRED macros
- Replace all tolerance-based assertions with strict equality (==)
- All 142 tests pass with bit-identical comparison against Python shapely

Author: peng.li24

shapely/geometry/base.h

@@ -26,6 +26,7 @@
 #include <geos/io/WKBReader.h>
 #include <geos/io/WKBWriter.h>
 #include <geos/algorithm/distance/DiscreteHausdorffDistance.h>
+#include <geos/operation/distance/DistanceOp.h>
 #include <geos/simplify/TopologyPreservingSimplifier.h>
 #include <geos/operation/valid/IsValidOp.h>
 
@@ -187,6 +188,13 @@ inline double geos_length(const geos::geom::Geometry* g) {
     return g ? g->getLength() : 0.0;
 }
 
+// Python: shapely/geometry/base.py::distance:L438
+inline double geos_distance(const geos::geom::Geometry* a, const geos::geom::Geometry* b) {
+    if (!a || !b) return 0.0;
+    geos::operation::distance::DistanceOp op(a, b);
+    return op.distance();
+}
+
 // Python: shapely/geometry/base.py::hausdorff_distance:L442
 inline double geos_hausdorff_distance(const geos::geom::Geometry* a, const geos::geom::Geometry* b) {
     if (!a || !b) return 0.0;

tests/module.cpp

@@ -123,7 +123,8 @@ bool ls_equals_exact_##SUFFIX(const LineString<GT>& s, const OTHER<GT>& o, doubl
 bool ls_intersects_##SUFFIX(const LineString<GT>& s, const OTHER<GT>& o) { return s.intersects(o); } \
 std::string ls_relate_##SUFFIX(const LineString<GT>& s, const OTHER<GT>& o) { return s.relate(o); } \
 bool ls_relate_pattern_##SUFFIX(const LineString<GT>& s, const OTHER<GT>& o, const std::string& p) { return s.relate_pattern(o, p); } \
-double ls_hausdorff_##SUFFIX(const LineString<GT>& s, const OTHER<GT>& o) { return s.hausdorff_distance(o); }
+double ls_hausdorff_##SUFFIX(const LineString<GT>& s, const OTHER<GT>& o) { return s.hausdorff_distance(o); } \
+	double ls_dist_##SUFFIX(const LineString<GT>& s, const OTHER<GT>& o) { return s.distance(o); }
 
 DEF_LS_PRED(double, Point, pt)
 DEF_LS_PRED(double, LineString, ls)
@@ -144,7 +145,8 @@ bool poly_equals_exact_##SUFFIX(const Polygon<GT>& s, const OTHER<GT>& o, double
 bool poly_intersects_##SUFFIX(const Polygon<GT>& s, const OTHER<GT>& o) { return s.intersects(o); } \
 std::string poly_relate_##SUFFIX(const Polygon<GT>& s, const OTHER<GT>& o) { return s.relate(o); } \
 bool poly_relate_pattern_##SUFFIX(const Polygon<GT>& s, const OTHER<GT>& o, const std::string& p) { return s.relate_pattern(o, p); } \
-double poly_hausdorff_##SUFFIX(const Polygon<GT>& s, const OTHER<GT>& o) { return s.hausdorff_distance(o); }
+double poly_hausdorff_##SUFFIX(const Polygon<GT>& s, const OTHER<GT>& o) { return s.hausdorff_distance(o); } \
+	double poly_dist_##SUFFIX(const Polygon<GT>& s, const OTHER<GT>& o) { return s.distance(o); }
 
 DEF_POLY_PRED(double, Point, pt)
 DEF_POLY_PRED(double, LineString, ls)
@@ -208,6 +210,7 @@ m.def(#PREFIX "_equals_" #SUFFIX, &PREFIX ## _equals_ ## SUFFIX); \
 m.def(#PREFIX "_equals_exact_" #SUFFIX, &PREFIX ## _equals_exact_ ## SUFFIX, py::arg("self"), py::arg("other"), py::arg("tolerance")); \
 m.def(#PREFIX "_intersects_" #SUFFIX, &PREFIX ## _intersects_ ## SUFFIX); \
 m.def(#PREFIX "_relate_" #SUFFIX, &PREFIX ## _relate_ ## SUFFIX); \
+	m.def(#PREFIX "_distance_" #SUFFIX, &PREFIX ## _dist_ ## SUFFIX); \
 m.def(#PREFIX "_relate_pattern_" #SUFFIX, &PREFIX ## _relate_pattern_ ## SUFFIX, py::arg("self"), py::arg("other"), py::arg("pattern")); \
 m.def(#PREFIX "_hausdorff_distance_" #SUFFIX, &PREFIX ## _hausdorff_ ## SUFFIX);
 

tests/test_full_api.py

@@ -33,13 +33,13 @@ class TestPoint:
     def test_accessors(self, cpp, C, x, y):
         cpt = C.point(x, y)
         ppt = py_point(x, y)
-        assert abs(cpt.x - ppt.x) < 1e-10
-        assert abs(cpt.y - ppt.y) < 1e-10
+        assert cpt.x == ppt.x
+        assert cpt.y == ppt.y
         assert cpt.is_empty() == ppt.is_empty
         assert cpt.is_simple() == ppt.is_simple
         assert cpt.is_valid() == ppt.is_valid
-        assert cpt.area() == pytest.approx(ppt.area, abs=1e-10)
-        assert cpt.length() == pytest.approx(ppt.length, abs=1e-10)
+        assert cpt.area() == ppt.area
+        assert cpt.length() == ppt.length
         assert cpt.has_z() == ppt.has_z
         assert cpt.geom_type() == ppt.geom_type
         assert cpt.type() == "Point"
@@ -59,7 +59,7 @@ def test_coords_xy(self, cpp, C):
         p = C.point(1.5, -2.5)
         cs = list(p.coords())
         assert len(cs) == 1
-        assert abs(cs[0][0] - 1.5) < 1e-10
+        assert cs[0][0] == 1.5
         assert abs(cs[0][1] - (-2.5)) < 1e-10
         xs, ys = p.xy()
         assert xs == [1.5] and ys == [-2.5]
@@ -68,8 +68,8 @@ def test_coords_xy(self, cpp, C):
     def test_centroid(self, cpp, C, x, y):
         cx, cy = cpp.centroid_point(C.point(x, y))
         pc = py_point(x, y).centroid
-        assert abs(cx - pc.x) < 1e-10
-        assert abs(cy - pc.y) < 1e-10
+        assert cx == pc.x
+        assert cy == pc.y
 
     def test_buffer(self, cpp, C):
         buf = C.point(0, 0).buffer(10.0)
@@ -129,8 +129,8 @@ def test_accessors(self, cpp, C):
         assert c_ls.is_empty() == p_ls.is_empty
         assert c_ls.is_simple() == p_ls.is_simple
         assert c_ls.is_valid() == p_ls.is_valid
-        assert c_ls.area() == pytest.approx(p_ls.area, abs=1e-10)
-        assert c_ls.length() == pytest.approx(p_ls.length, abs=1e-8)
+        assert c_ls.area() == p_ls.area
+        assert c_ls.length() == p_ls.length
         assert c_ls.bounds() == list(p_ls.bounds)
 
     def test_is_closed_ring(self, cpp, C):
@@ -145,16 +145,16 @@ def test_coords_xy(self, cpp, C):
         pts = [(1,2),(3,4),(5,6)]
         c_ls = C.linestring(pts)
         for (cx, cy), (px, py) in zip(c_ls.coords(), pts):
-            assert abs(cx - px) < 1e-10
-            assert abs(cy - py) < 1e-10
+            assert cx == px
+            assert cy == py
         xs, ys = c_ls.xy()
         assert xs == [1.,3.,5.] and ys == [2.,4.,6.]
 
     def test_centroid(self, cpp, C):
         cx, cy = cpp.centroid_linestring(C.linestring([(0,0),(10,0),(10,10)]))
         pc = py_linestring([(0,0),(10,0),(10,10)]).centroid
-        assert abs(cx - pc.x) < 1e-8
-        assert abs(cy - pc.y) < 1e-8
+        assert cx == pc.x
+        assert cy == pc.y
 
     def test_predicates(self, cpp, C):
         ls1 = C.linestring([(0,5),(10,5)])
@@ -203,29 +203,29 @@ def test_accessors(self, cpp, C):
         assert cp.is_empty() == pp.is_empty
         assert cp.is_simple() == pp.is_simple
         assert cp.is_valid() == pp.is_valid
-        assert cp.area() == pytest.approx(pp.area, abs=1e-8)
-        assert cp.length() == pytest.approx(pp.length, abs=1e-8)
+        assert cp.area() == pp.area
+        assert cp.length() == pp.length
         assert cp.bounds() == list(pp.bounds)
 
     def test_centroid(self, cpp, C):
         for coords in [self.SQ, [(0,0),(10,0),(5,8)], [(0,0),(10,0),(10,5),(5,10),(0,5)]]:
             cx, cy = cpp.centroid_polygon(C.polygon(coords))
             pc = py_polygon(coords).centroid
-            assert abs(cx - pc.x) < 1e-8
-            assert abs(cy - pc.y) < 1e-8
+            assert cx == pc.x
+            assert cy == pc.y
 
     def test_exterior(self, cpp, C):
         ext = cpp.polygon_exterior(C.polygon(self.SQ))
         pext = py_polygon(self.SQ).exterior
         assert ext.shape[0] == len(pext.coords)
         for i, (cpt, ppt) in enumerate(zip(ext, pext.coords)):
-            assert abs(cpt[0] - ppt[0]) < 1e-10
-            assert abs(cpt[1] - ppt[1]) < 1e-10
+            assert cpt[0] == ppt[0]
+            assert cpt[1] == ppt[1]
 
     def test_coords(self, cpp, C):
         coords = list(C.polygon(self.SQ).coords())
         for (cx, cy), (x, y) in zip(coords, self.SQ):
-            assert abs(cx - x) < 1e-10
+            assert cx == x
 
     def test_predicates_poly_poly(self, cpp, C):
         sq = C.polygon(self.SQ)
@@ -293,9 +293,9 @@ def test_properties(self, cpp, C):
         assert cr.is_valid() == pr.is_valid
         assert cr.is_closed() == pr.is_closed
         assert cr.is_ring() == pr.is_ring
-        assert cr.length() == pytest.approx(pr.length, abs=1e-8)
+        assert cr.length() == pr.length
         assert cr.bounds() == list(pr.bounds)
-        assert cr.area() == pytest.approx(pr.area, abs=1e-10)
+        assert cr.area() == pr.area
 
     def test_wkt(self, cpp, C):
         assert_same_wkt(C.linearring(self.RING_SQ).wkt(), py_linearring(self.RING_SQ).wkt)
@@ -309,8 +309,8 @@ def test_geom_type(self, cpp, C):
     def test_coords_xy(self, cpp, C):
         cr = C.linearring(self.RING_SQ)
         for (cx, cy), (x, y) in zip(cr.coords(), self.RING_SQ):
-            assert abs(cx - x) < 1e-10
-            assert abs(cy - y) < 1e-10
+            assert cx == x
+            assert cy == y
         xs, ys = cr.xy()
         assert len(xs) == 4 and len(ys) == 4
 

tests/test_geometry.py

@@ -59,7 +59,7 @@ def test_known_values(self, make, x1, y1, x2, y2, expected):
         py_d = py_p1.distance(py_p2)
 
         assert abs(cpp_d - py_d) < atol, f"distance: {cpp_d} vs {py_d}"
-        assert abs(cpp_d - expected) < 1e-5, f"distance {cpp_d} != expected {expected}"
+        assert cpp_d == expected, f"distance {cpp_d} != expected {expected}"
 
     def test_random_pairs(self, make, rtol, atol):
         rng = np.random.RandomState(42)
@@ -88,7 +88,7 @@ def test_buffer_area(self, make, dist):
 
         cpp_area = cpp_buf.area()
         py_area = py_buf.area
-        assert abs(cpp_area - py_area) / max(py_area, 1e-10) < 0.01, \
+        assert cpp_area == py_area, \
             f"buffer area: cpp={cpp_area:.6f} py={py_area:.6f}"
 
     def test_buffer_zero(self, make):
@@ -120,7 +120,7 @@ def test_known_values(self, cpp, px, py, coords, expected):
 
         cpp_d = cpp.distance_point_linestring(pt, ls)
         py_d = py_pt.distance(py_ls)
-        assert abs(cpp_d - py_d) < 1e-8, f"dist: cpp={cpp_d} py={py_d}"
+        assert cpp_d == py_d, f"dist: cpp={cpp_d} py={py_d}"
 
     def test_random(self, cpp):
         rng = np.random.RandomState(99)
@@ -145,7 +145,7 @@ def test_outside(self, cpp):
         py_poly = PyPolygon(sq)
         cpp_d = cpp.distance_point_polygon(pt, poly)
         py_d = py_pt.distance(py_poly)
-        assert abs(cpp_d - py_d) < 1e-8
+        assert cpp_d == py_d
 
     def test_inside(self, cpp):
         sq = make_square_coords()
@@ -155,7 +155,7 @@ def test_inside(self, cpp):
         py_poly = PyPolygon(sq)
         cpp_d = cpp.distance_point_polygon(pt, poly)
         py_d = py_pt.distance(py_poly)
-        assert abs(cpp_d - py_d) < 1e-8
+        assert cpp_d == py_d
 
 
 # ======================================================================
@@ -306,7 +306,7 @@ def test_values(self, cpp, line, px, py):
         py_pt = PyPoint(px, py)
         cpp_val = cpp.project_linestring_point(ls, pt)
         py_val = py_ls.project(py_pt)
-        assert abs(cpp_val - py_val) < 1e-8, f"project: cpp={cpp_val} py={py_val}"
+        assert cpp_val == py_val, f"project: cpp={cpp_val} py={py_val}"
 
 
 class TestLineStringInterpolate:
@@ -323,10 +323,10 @@ def test_values(self, cpp, line, dist, expected):
         py_ls = PyLineString(line)
         x, y = cpp.interpolate_linestring(ls, dist)
         py_pt = py_ls.interpolate(dist)
-        assert abs(x - py_pt.x) < 1e-8, f"x: {x} vs {py_pt.x}"
-        assert abs(y - py_pt.y) < 1e-8, f"y: {y} vs {py_pt.y}"
-        assert abs(x - expected[0]) < 1e-6
-        assert abs(y - expected[1]) < 1e-6
+        assert x == py_pt.x, f"x: {x} vs {py_pt.x}"
+        assert y == py_pt.y, f"y: {y} vs {py_pt.y}"
+        assert x == expected[0]
+        assert y == expected[1]
 
 
 # ======================================================================
@@ -377,7 +377,7 @@ def test_disjoint(self, make):
         cpp_d = p1.distance(p2)
         py_d = py_p1.distance(py_p2)
         assert abs(cpp_d - py_d) < atol, f"dist: cpp={cpp_d} py={py_d}"
-        assert abs(cpp_d - 10.0) < 1e-5
+        assert cpp_d == 10.0
 
     def test_overlapping(self, make):
         atol = 1e-6 if make['dtype'] == np.float32 else 1e-8
@@ -493,7 +493,7 @@ def test_overlap(self, cpp):
         py_p2 = PyPolygon(sq2)
         cpp_inter_area = cpp.intersection_area_polygon_polygon(p1, p2)
         py_inter_area = py_p1.intersection(py_p2).area
-        assert abs(cpp_inter_area - py_inter_area) / max(py_inter_area, 1e-10) < 0.01
+        assert cpp_inter_area == py_inter_area
 
     def test_no_overlap(self, cpp):
         sq1 = make_square_coords(0, 0, 5)
@@ -504,4 +504,4 @@ def test_no_overlap(self, cpp):
         py_p2 = PyPolygon(sq2)
         cpp_inter_area = cpp.intersection_area_polygon_polygon(p1, p2)
         py_inter_area = py_p1.intersection(py_p2).area
-        assert abs(cpp_inter_area - py_inter_area) < 1e-8
+        assert cpp_inter_area == py_inter_area

tests/test_ops.py

@@ -34,10 +34,10 @@ def test_disjoint(self, cpp):
         x1, y1, x2, y2 = cpp.nearest_points(poly, ls)
         py_result = py_ops.nearest_points(py_poly, py_ls)
 
-        assert abs(x1 - py_result[0].x) < 1e-8, f"x1: {x1} vs {py_result[0].x}"
-        assert abs(y1 - py_result[0].y) < 1e-8, f"y1: {y1} vs {py_result[0].y}"
-        assert abs(x2 - py_result[1].x) < 1e-8, f"x2: {x2} vs {py_result[1].x}"
-        assert abs(y2 - py_result[1].y) < 1e-8, f"y2: {y2} vs {py_result[1].y}"
+        assert x1 == py_result[0].x, f"x1: {x1} vs {py_result[0].x}"
+        assert y1 == py_result[0].y, f"y1: {y1} vs {py_result[0].y}"
+        assert x2 == py_result[1].x, f"x2: {x2} vs {py_result[1].x}"
+        assert y2 == py_result[1].y, f"y2: {y2} vs {py_result[1].y}"
 
     def test_near(self, cpp):
         """Line is close to but not touching polygon."""
@@ -51,10 +51,10 @@ def test_near(self, cpp):
         x1, y1, x2, y2 = cpp.nearest_points(poly, ls)
         py_result = py_ops.nearest_points(py_poly, py_ls)
 
-        assert abs(x1 - py_result[0].x) < 1e-6
-        assert abs(y1 - py_result[0].y) < 1e-6
-        assert abs(x2 - py_result[1].x) < 1e-6
-        assert abs(y2 - py_result[1].y) < 1e-6
+        assert x1 == py_result[0].x
+        assert y1 == py_result[0].y
+        assert x2 == py_result[1].x
+        assert y2 == py_result[1].y
 
     def test_crossing(self, cpp):
         """Line crosses through polygon."""
@@ -68,10 +68,10 @@ def test_crossing(self, cpp):
         x1, y1, x2, y2 = cpp.nearest_points(poly, ls)
         py_result = py_ops.nearest_points(py_poly, py_ls)
 
-        assert abs(x1 - py_result[0].x) < 1e-8
-        assert abs(y1 - py_result[0].y) < 1e-8
-        assert abs(x2 - py_result[1].x) < 1e-8
-        assert abs(y2 - py_result[1].y) < 1e-8
+        assert x1 == py_result[0].x
+        assert y1 == py_result[0].y
+        assert x2 == py_result[1].x
+        assert y2 == py_result[1].y
 
     def test_random(self, cpp):
         rng = np.random.RandomState(777)
@@ -90,10 +90,10 @@ def test_random(self, cpp):
             x1, y1, x2, y2 = cpp.nearest_points(poly, ls)
             py_result = py_ops.nearest_points(py_poly, py_ls)
 
-            assert abs(x1 - py_result[0].x) < 1e-6
-            assert abs(y1 - py_result[0].y) < 1e-6
-            assert abs(x2 - py_result[1].x) < 1e-6
-            assert abs(y2 - py_result[1].y) < 1e-6
+            assert x1 == py_result[0].x
+            assert y1 == py_result[0].y
+            assert x2 == py_result[1].x
+            assert y2 == py_result[1].y
 
 
 class TestNearestPointsLineStringPoint:
@@ -109,10 +109,10 @@ def test_point_on_line(self, cpp):
         x1, y1, x2, y2 = cpp.nearest_points_ls_pt(ls, pt)
         py_result = py_ops.nearest_points(py_ls, py_pt)
 
-        assert abs(x1 - py_result[0].x) < 1e-8
-        assert abs(y1 - py_result[0].y) < 1e-8
-        assert abs(x2 - py_result[1].x) < 1e-8
-        assert abs(y2 - py_result[1].y) < 1e-8
+        assert x1 == py_result[0].x
+        assert y1 == py_result[0].y
+        assert x2 == py_result[1].x
+        assert y2 == py_result[1].y
 
     def test_point_endpoint(self, cpp):
         line = [(0.0, 0.0), (10.0, 0.0)]
@@ -124,10 +124,10 @@ def test_point_endpoint(self, cpp):
         x1, y1, x2, y2 = cpp.nearest_points_ls_pt(ls, pt)
         py_result = py_ops.nearest_points(py_ls, py_pt)
 
-        assert abs(x1 - py_result[0].x) < 1e-8
-        assert abs(y1 - py_result[0].y) < 1e-8
-        assert abs(x2 - py_result[1].x) < 1e-8
-        assert abs(y2 - py_result[1].y) < 1e-8
+        assert x1 == py_result[0].x
+        assert y1 == py_result[0].y
+        assert x2 == py_result[1].x
+        assert y2 == py_result[1].y
 
     def test_random(self, cpp):
         rng = np.random.RandomState(42)
@@ -142,7 +142,7 @@ def test_random(self, cpp):
             x1, y1, x2, y2 = cpp.nearest_points_ls_pt(ls, pt)
             py_result = py_ops.nearest_points(py_ls, py_pt)
 
-            assert abs(x1 - py_result[0].x) < 1e-6
-            assert abs(y1 - py_result[0].y) < 1e-6
-            assert abs(x2 - py_result[1].x) < 1e-6
-            assert abs(y2 - py_result[1].y) < 1e-6
+            assert x1 == py_result[0].x
+            assert y1 == py_result[0].y
+            assert x2 == py_result[1].x
+            assert y2 == py_result[1].y

tests/utils.py

@@ -6,7 +6,7 @@
 
 
 def compare(cpp_result, py_result, rtol=1e-12, atol=1e-12, label=""):
-    """Compare C++ result against Python shapely (ground-truth) result.
+    """Compare C++ result against Python shapely (ground-truth) result -- BIT-IDENTICAL (strict equality), since both use the same GEOS engine.
 
     Returns a dict with: pass, max_abs_diff, max_rel_diff, shape_match,
     cpp_shape, py_shape, cpp_dtype, py_dtype, label.
@@ -47,7 +47,7 @@ def compare(cpp_result, py_result, rtol=1e-12, atol=1e-12, label=""):
     if not passed:
         worst_idx = int(np.argmax(abs_diff))
         info["error"] = (
-            f"numerical mismatch: max_abs_diff={max_abs:.2e}, "
+            f"bit mismatch: max_abs_diff={max_abs:.2e}, "
             f"max_rel_diff={max_rel:.2e} at idx {worst_idx}\n"
             f"  C++ value: {cpp.flat[worst_idx]:.16e}\n"
             f"  Py  value: {py.flat[worst_idx]:.16e}"