github-actions[bot] commented on code in PR #69032:
URL: https://github.com/apache/doris/pull/69032#discussion_r4376498668
##########
be/src/exprs/function/geo/functions_geo.cpp:
##########
@@ -826,142 +2103,114 @@ struct StIsClosed {
}
};
-struct StDistance {
- static constexpr auto NAME = "st_distance";
- static const size_t NUM_ARGS = 2;
- using Type = DataTypeFloat64;
+static double geography_distance(const GeoShape& left_shape, const GeoShape&
right_shape,
+ StringRef left_wkb, StringRef right_wkb) {
+ // Keep the common point-pair path allocation-free while canonicalizing
date-line/pole aliases.
+ if (left_shape.type() == GeoShapeType::GEO_SHAPE_POINT &&
+ right_shape.type() == GeoShapeType::GEO_SHAPE_POINT) {
+ return S2Earth::GetDistanceMeters(
+ canonical_geography_point(*assert_cast<const
GeoPoint&>(left_shape).point()),
+ canonical_geography_point(*assert_cast<const
GeoPoint&>(right_shape).point()));
+ }
+ const GeographyEdgeNodes nodes(left_wkb, right_wkb);
+ MutableS2ShapeIndex left;
+ MutableS2ShapeIndex right;
+ geography_shape_index(left_shape, nodes, &left);
+ geography_shape_index(right_shape, nodes, &right);
+ S2ClosestEdgeQuery query(&left);
+ S2ClosestEdgeQuery::ShapeIndexTarget target(&right);
+ return S2Earth::ToMeters(query.GetDistance(&target).ToAngle());
Review Comment:
[P2] Stabilize near-antipodal raw arcs before distance queries
For default Geography WKB representing A=`LINESTRING(-44.99999999999997
-45,134.99999999999997 45)` and P=`POINT(45 0.000001)`, both inputs are valid
and A's endpoints are not antipodal. Its raw-degree midpoint is exactly (45,0),
so the true distance to P is at most 0.112 m and `ST_DWithin(A,P,1)` must be
true. However, P is off A's exact plane, so `collect_circle_vertices` adds no
subdivision and this query receives the original independently rounded S2 edge.
Pinned S2's [ToPoint
conversion](https://github.com/google/s2geometry/blob/v0.10.0/src/s2/s2latlng.cc#L37-L44)
and [edge-distance
calculation](https://github.com/google/s2geometry/blob/v0.10.0/src/s2/s2edge_distances.cc#L154-L187)
measure the distorted rounded plane: a static arithmetic model gives about 237
km, making DWithin false. Adding the exact midpoint to A changes this result
despite preserving its spherical point set.
The existing membership/equality fixes do not cover this nearby off-arc
proximity path. Stabilize the raw-arc representation independently of vertices
contributed by the other operand, and cover Distance/DWithin in both orders and
under equivalent subdivision. Evidence is source/arithmetic analysis; no Doris
runtime test was run.
##########
be/src/exprs/function/geo/functions_geo.cpp:
##########
@@ -18,32 +18,1030 @@
#include "exprs/function/geo/functions_geo.h"
#include <glog/logging.h>
+#include <s2/mutable_s2shape_index.h>
+#include <s2/s2boolean_operation.h>
+#include <s2/s2builderutil_s2polyline_vector_layer.h>
+#include <s2/s2closest_edge_query.h>
+#include <s2/s2closest_point_query.h>
+#include <s2/s2earth.h>
+#include <s2/s2lax_polygon_shape.h>
+#include <s2/s2loop.h>
#include <s2/s2point.h>
+#include <s2/s2point_index.h>
+#include <s2/s2point_vector_shape.h>
+#include <s2/s2pointutil.h>
+#include <s2/s2polygon.h>
+#include <s2/s2polyline.h>
+#include <s2/s2predicates.h>
#include <algorithm>
+#include <array>
+#include <bit>
#include <boost/iterator/iterator_facade.hpp>
+#include <boost/multiprecision/cpp_int.hpp>
+#include <cmath>
+#include <cstdint>
+#include <limits>
+#include <map>
+#include <optional>
+#include <set>
#include <utility>
+#include "common/cast_set.h"
#include "common/compiler_util.h"
#include "core/assert_cast.h"
#include "core/block/block.h"
#include "core/block/column_with_type_and_name.h"
#include "core/column/column.h"
#include "core/column/column_array.h"
#include "core/column/column_execute_util.h"
+#include "core/column/column_map.h"
#include "core/column/column_nullable.h"
+#include "core/column/column_struct.h"
+#include "core/column/column_varbinary.h"
+#include "core/custom_allocator.h"
#include "core/data_type/data_type_array.h"
+#include "core/data_type/data_type_map.h"
#include "core/data_type/data_type_nullable.h"
#include "core/data_type/data_type_number.h"
+#include "core/data_type/data_type_spatial.h"
#include "core/data_type/data_type_string.h"
+#include "core/data_type/data_type_struct.h"
+#include "core/data_type/data_type_varbinary.h"
#include "core/data_type/define_primitive_type.h"
#include "core/string_ref.h"
#include "exprs/function/geo/geo_common.h"
#include "exprs/function/geo/geo_types.h"
+#include "exprs/function/geo/wkb_parse.h"
#include "exprs/function/simple_function_factory.h"
+#include "exprs/function/string_hex_util.h"
namespace doris {
+static bool is_spatial_type(const DataTypePtr& type) {
+ const auto primitive_type = remove_nullable(type)->get_primitive_type();
+ return primitive_type == TYPE_GEOMETRY || primitive_type == TYPE_GEOGRAPHY;
+}
+
+static Status validate_geography_semantics(const DataTypePtr& type, const
char* function_name) {
+ const auto& nested_type = remove_nullable(type);
+ if (!is_spatial_type(nested_type)) {
+ return Status::OK();
+ }
+
+ const auto* spatial_type = dynamic_cast<const
DataTypeSpatial*>(nested_type.get());
+ DCHECK(spatial_type != nullptr);
+ if (spatial_type->get_primitive_type() == TYPE_GEOGRAPHY &&
+ spatial_type->crs() == "OGC:CRS84" && spatial_type->algorithm() ==
"spherical") {
+ return Status::OK();
+ }
+ return Status::NotSupported(
+ "Function {} requires GEOGRAPHY(OGC:CRS84, spherical) for spatial
inputs",
+ function_name);
+}
+
+static std::unique_ptr<GeoShape> decode_geo_shape(StringRef value, const
DataTypePtr& type) {
+ if (!is_spatial_type(type)) {
+ return GeoShape::from_encoded(value.data, value.size);
+ }
+
+ GeoParseStatus status;
+ auto shape = GeoShape::from_wkb_bytes(value.data, value.size, status);
+ return status == GeoParseStatus::GEO_PARSE_OK ? std::move(shape) : nullptr;
+}
+
+static bool any_geo_argument_is_null(const Block& block, const ColumnNumbers&
arguments,
+ size_t row) {
+ return std::ranges::any_of(arguments, [&](auto argument) {
+ return block.get_by_position(argument).column->is_null_at(row);
+ });
+}
+
+static Status decode_calculation_shape(StringRef value, const DataTypePtr&
type,
+ std::unique_ptr<GeoShape>* shape,
+ const char* function_name) {
+ *shape = decode_geo_shape(value, type);
+ // POLYGON EMPTY decodes successfully, but distance helpers retain
infinity/max sentinels.
+ // Reject it for typed calculations just like empty points/lines, leaving
raw WKB untouched.
+ if (is_spatial_type(type) &&
+ (*shape == nullptr ||
+ ((*shape)->type() == GeoShapeType::GEO_SHAPE_POLYGON &&
+ assert_cast<const GeoPolygon&>(**shape).polygon()->num_loops() ==
0))) {
+ return Status::NotSupported("Function {} does not support this
Geography WKB shape",
+ function_name);
+ }
+ return Status::OK();
+}
+
+// Calculations use canonical positions; raw spatial accessors still return
the original WKB.
+static S2Point canonical_geography_point(const S2Point& point) {
+ const auto coordinate = S2LatLng(point);
+ const double latitude = coordinate.lat().radians();
+ // All longitudes at a pole and +/-180 at the date line denote the same
position.
+ if (std::abs(latitude) == M_PI_2) {
+ return {0, 0, latitude > 0 ? 1.0 : -1.0};
+ }
+ if (coordinate.lng().radians() == M_PI) {
+ // Flip only sin(longitude): rebuilding latitude can change its last
bit.
+ return {point.x(), -point.y(), point.z()};
+ }
+ return point;
+}
+
+// Error-free TwoSum preserves tiny differences from a half-turn (e.g. -1e-15
+ 180).
+static std::pair<double, double> geography_degree_sum(double a, double b) {
+ const double high = a + b;
+ const double offset = high - a;
+ return {high, (a - (high - offset)) + (b - offset)};
+}
+
+#ifdef BE_TEST
+static thread_local size_t geography_exact_proofs = 0;
+#endif
+
+// Test coplanarity of the original degree ordinates, not independently
rounded S2 vectors.
+// For example (45,-45), (90,0), (135,45) share an arc although ToPoint's
determinant is nonzero.
+static bool geography_common_great_circle(const GeoCoordinate& a, const
GeoCoordinate& b,
+ const GeoCoordinate& c) {
+#ifdef BE_TEST
+ ++geography_exact_proofs;
+#endif
+ using boost::multiprecision::cpp_int;
+ const std::array<double, 6> ordinates {a.y, a.x, b.y, b.x, c.y, c.x};
+ std::array<uint64_t, 6> significands;
+ std::array<int, 6> exponents;
+ int denominator_bits = 0;
+ for (size_t i = 0; i < ordinates.size(); ++i) {
+ const auto bits = std::bit_cast<uint64_t>(ordinates[i]);
+ const int biased = (bits >> 52) & 0x7ff;
+ uint64_t significand = bits & ((uint64_t {1} << 52) - 1);
+ int exponent = biased ? biased - 1075 : -1074;
+ if (biased) {
+ significand |= uint64_t {1} << 52;
+ }
+ if (significand != 0) {
+ const int zeros = std::countr_zero(significand);
+ significand >>= zeros;
+ exponent += zeros;
+ denominator_bits = std::max(denominator_bits, -exponent);
+ } else {
+ exponent = 0;
+ }
+ significands[i] = significand;
+ exponents[i] = exponent;
+ }
+ std::array<cpp_int, 6> angles;
+ for (size_t i = 0; i < angles.size(); ++i) {
+ angles[i] = cpp_int(significands[i]) << (denominator_bits +
exponents[i]);
+ if (std::signbit(ordinates[i])) {
+ angles[i] = -angles[i];
+ }
+ }
+ // Binary-rational degrees have denominator 2^q. Exp(i*angle) is a power
of a
+ // root of unity of order 360*2^q. Its minimal polynomial is Phi_90(t),
with
+ // t = z^(2^(q+2)), and Phi_90(t) = t^24+t^21-t^15-t^12-t^9+t^3+1.
+ // Expand cos(lat_i)cos(lat_j)sin(lat_k)sin(lon_j-lon_i), cyclically, into
+ // 48 integer-coefficient powers. Equal residues modulo 2^(q+2) share one
+ // polynomial; reduction by Phi_90 proves zero exactly, including tiny
nonzero controls.
+ const int shift = denominator_bits + 2;
+ const cpp_int basis = cpp_int(1) << shift;
+ const cpp_int period = 90 * basis;
+ std::map<cpp_int, std::array<int, 90>> polynomials;
+ for (int i = 0; i < 3; ++i) {
+ const int j = (i + 1) % 3;
+ const int k = (i + 2) % 3;
+ for (int mask = 0; mask < 16; ++mask) {
+ std::array<int, 4> signs;
+ for (int bit = 0; bit < 4; ++bit) {
+ signs[bit] = (mask & (1 << bit)) ? 1 : -1;
+ }
+ cpp_int power = signs[0] * angles[2 * i] + signs[1] * angles[2 *
j] +
+ signs[2] * angles[2 * k] +
+ signs[3] * (angles[2 * j + 1] - angles[2 * i + 1]);
+ power %= period;
+ if (power < 0) {
+ power += period;
+ }
+ const int quotient = (power >> shift).convert_to<int>();
+ const cpp_int residue = power - (cpp_int(quotient) << shift);
+ polynomials[residue][quotient] += signs[2] * signs[3];
+ }
+ }
+ for (auto& [residue, coefficients] : polynomials) {
+ for (int i = 89; i >= 24; --i) {
+ const int coefficient = coefficients[i];
+ coefficients[i] = 0;
+ coefficients[i - 3] -= coefficient;
+ coefficients[i - 9] += coefficient;
+ coefficients[i - 12] += coefficient;
+ coefficients[i - 15] += coefficient;
+ coefficients[i - 21] -= coefficient;
+ coefficients[i - 24] -= coefficient;
+ }
+ if (std::ranges::any_of(coefficients, [](int value) { return value !=
0; })) {
+ return false;
+ }
+ }
+ return true;
+}
+
+// Keep each atomic segment in its own shape. Pinned S2 ProcessEdge1 uses a
chain's
+// global start as a local offset in one endpoint path, so multiple chains per
line
+// shape are unsafe. Closed AA shapes restore internal vertices in the OPEN
model.
+class GeographyLineShape final : public S2Shape {
+public:
+ GeographyLineShape(const S2Point& a, const S2Point& b) {
_edge.emplace_back(a, b); }
+ int num_edges() const override { return 1; }
+ Edge edge(int i) const override {
+ DCHECK_EQ(i, 0);
+ return _edge.front();
+ }
+ int dimension() const override { return 1; }
+ ReferencePoint GetReferencePoint() const override { return
ReferencePoint::Contained(false); }
+ int num_chains() const override { return 1; }
+ Chain chain(int i) const override {
+ DCHECK_EQ(i, 0);
+ return {0, 1};
+ }
+ Edge chain_edge(int i, int j) const override {
+ DCHECK_EQ(i, 0);
+ return edge(j);
+ }
+ ChainPosition chain_position(int i) const override {
+ DCHECK_EQ(i, 0);
+ return {0, 0};
+ }
+
+private:
+ // Owning storage is checked/accounted by the current query's allocator.
+ DorisVector<Edge> _edge;
+};
+
+class GeographyEdgeNodes {
+public:
+ GeographyEdgeNodes(StringRef left, StringRef right) {
+ GeoCoordinateList coordinates;
+ const auto left_status = WkbParse::xy_coordinates(left.data,
left.size, &coordinates);
+ const auto right_status = WkbParse::xy_coordinates(right.data,
right.size, &coordinates);
+ DCHECK_EQ(left_status, GeoParseStatus::GEO_PARSE_OK);
+ DCHECK_EQ(right_status, GeoParseStatus::GEO_PARSE_OK);
+ for (auto coordinate : coordinates.list) {
+ if (coordinate.x == 180) {
+ coordinate.x = -180;
+ }
+ const auto point = canonical_geography_point(
+ S2LatLng::FromDegrees(coordinate.y,
coordinate.x).ToPoint());
+ if (!_coordinates.emplace(point, coordinate).second) {
+ continue;
+ }
+ _vertex_index.Add(point, 0);
+ if (std::abs(coordinate.y) != 90) {
+ _meridians[meridian(coordinate.x)][coordinate.x < 0].push_back(
+ {coordinate.y, point});
+ }
+ }
+ for (auto& entry : _meridians) {
+ for (auto& side : entry.second) {
+ side.push_back({-90, {0, 0, -1}});
+ side.push_back({90, {0, 0, 1}});
+ std::ranges::sort(side, {}, &Vertex::latitude);
+ }
+ }
+ }
+
+ bool contains_point(const S2Point& raw_start, const S2Point& raw_end,
+ const S2Point& raw_point) const {
+ const auto start = canonical_geography_point(raw_start);
+ const auto end = canonical_geography_point(raw_end);
+ const auto point = canonical_geography_point(raw_point);
+ // Preserve exact endpoints, including date-line and pole aliases.
+ if (point == start || point == end) {
+ return true;
+ }
+ if (start == end) {
+ return false;
+ }
+ const auto& a = _coordinates.at(start);
+ const auto& b = _coordinates.at(end);
+ const auto& c = _coordinates.at(point);
+ const auto key = meridian(std::abs(a.y) == 90 ? b.x : a.x);
+ bool member = false;
+ if (meridian_edge(start, end, [&](const auto& vertices, double from,
double to) {
+ // Opposite-longitude edges can pass through a pole; the pole
belongs
+ // to both legs, whereas other positions must be on the
selected side.
+ const bool same_side =
+ std::abs(c.y) == 90 ||
+ (meridian(c.x) == key && &_meridians.at(key)[c.x < 0]
== &vertices);
+ member |= same_side && c.y >= std::min(from, to) && c.y <=
std::max(from, to);
+ })) {
+ return member;
+ }
+ const bool wraps = geography_degree_sum(std::max(a.x, b.x),
-std::min(a.x, b.x)) >
+ std::pair {180.0, 0.0};
+ const bool within = wraps ? (c.x < std::min(a.x, b.x) || c.x >
std::max(a.x, b.x))
+ : (c.x > std::min(a.x, b.x) && c.x <
std::max(a.x, b.x));
+ // Near-antipodal endpoint conversion can rotate the S2 plane
substantially.
+ // POINT(45 0) is on (-44.99999999999997 -45,134.99999999999997 45);
+ // decide from raw degrees rather than accepting a larger distance
tolerance.
+ return within && geography_common_great_circle(a, b, c);
+ }
+
+ void append_edge(const S2Point& raw_start, const S2Point& raw_end,
+ DorisVector<S2Point>* output) const {
+ const auto start = canonical_geography_point(raw_start);
+ const auto end = canonical_geography_point(raw_end);
+ output->push_back(start);
+ if (start == end) {
+ return;
+ }
+ if (meridian_edge(start, end, [&](const auto& vertices, double from,
double to) {
+ append_range(vertices, from, to, output);
+ if (std::abs(to) == 90) {
+ const S2Point pole(0, 0, to > 0 ? 1 : -1);
+ if (end != pole) {
+ output->push_back(pole);
+ }
+ }
+ })) {
+ return;
+ }
+ const auto& a = _coordinates.at(start);
+ const auto& b = _coordinates.at(end);
+ const bool wraps = geography_degree_sum(std::max(a.x, b.x),
-std::min(a.x, b.x)) >
+ std::pair {180.0, 0.0};
+ const bool east = (a.x < b.x) != wraps;
+ S2ClosestPointQuery<int> query(&_vertex_index);
+ const double conversion_error = 8 *
std::numeric_limits<double>::epsilon();
+ // Near antipodes, endpoint rounding rotates the edge plane by
O(epsilon/sin(length)).
+ // Widen candidate search only; the exact raw-angle proof still
decides membership.
+ // Every point on a minor arc is at most pi/2 from one endpoint,
including when
+ // rounded antipodes make the computed normal vanish.
+ const double search_radius =
+ start.DotProd(end) >= 0
+ ? conversion_error
+ : std::min(M_PI_2 + conversion_error,
+ 4 * conversion_error /
start.CrossProd(end).Norm());
+
query.mutable_options()->set_conservative_max_distance(S1Angle::Radians(search_radius));
+ S2ClosestPointQuery<int>::EdgeTarget target(start, end);
+ DorisVector<S2Point> on_edge;
+ for (const auto& candidate : query.FindClosestPoints(&target)) {
+ const auto& point = candidate.point();
+ const auto& c = _coordinates.at(point);
+ // Longitude is strictly monotone along every non-meridian minor
arc.
+ // Check the original interval (including date-line wrap), not a
distance
+ // tolerance that could append a point just beyond an endpoint.
+ const bool within = wraps ? (c.x < std::min(a.x, b.x) || c.x >
std::max(a.x, b.x))
+ : (c.x > std::min(a.x, b.x) && c.x <
std::max(a.x, b.x));
+ if (within && geography_common_great_circle(a, b, c)) {
+ on_edge.push_back(point);
+ }
+ }
+ std::ranges::sort(on_edge, [&](const S2Point& x, const S2Point& y) {
+ const double x_lon = _coordinates.at(x).x;
+ const double y_lon = _coordinates.at(y).x;
+ if ((x_lon < a.x) != (y_lon < a.x)) {
+ return east ? x_lon >= a.x : x_lon < a.x;
+ }
+ return east ? x_lon < y_lon : x_lon > y_lon;
+ });
+ output->insert(output->end(), on_edge.begin(), on_edge.end());
+ }
+
+ void add_line(const GeoLine& shape, MutableS2ShapeIndex* index) const {
+ using Edge = std::pair<S2Point, S2Point>;
+ using EdgeSet = std::set<Edge, std::less<>, CustomStdAllocator<Edge>>;
+ using Interval = std::pair<double, double>;
+ DorisMap<const DorisVector<Vertex>*, DorisVector<Interval>> intervals;
+ EdgeSet raw_edges;
+ const auto* line = shape.polyline();
+ const auto first = canonical_geography_point(line->vertex(0));
+ const auto last =
canonical_geography_point(line->vertex(line->num_vertices() - 1));
+ for (int i = 0; i + 1 < line->num_vertices(); ++i) {
+ const auto a = canonical_geography_point(line->vertex(i));
+ const auto b = canonical_geography_point(line->vertex(i + 1));
+ if (a == b) {
+ continue;
+ }
+ if (!meridian_edge(a, b, [&](const auto& vertices, double from,
double to) {
+ intervals[&vertices].emplace_back(std::min(from, to),
std::max(from, to));
+ })) {
+ raw_edges.emplace(std::min(a, b), std::max(a, b));
+ }
+ }
+ EdgeSet edges;
+ auto add_edge = [&](const S2Point& a, const S2Point& b) {
+ edges.emplace(std::min(a, b), std::max(a, b));
+ };
+ // Merge latitude coverage before splitting: N repeated or overlapping
edges
+ // and M collinear vertices produce O(N+M) scratch, never an N*M
vertex walk.
+ for (auto& [vertices, ranges] : intervals) {
+ std::ranges::sort(ranges);
+ auto emit = [&](double from, double to) {
+ auto begin = std::ranges::lower_bound(*vertices, from, {},
&Vertex::latitude);
+ const auto end = std::ranges::upper_bound(*vertices, to, {},
&Vertex::latitude);
+ for (auto next = begin + 1; next != end; ++begin, ++next) {
+ add_edge(begin->point, next->point);
+ }
+ };
+ double from = ranges.front().first;
+ double to = ranges.front().second;
+ for (const auto& range : ranges) {
+ if (range.first > to) {
+ emit(from, to);
+ from = range.first;
+ }
+ to = std::max(to, range.second);
+ }
+ emit(from, to);
+ }
+ add_circle_edges(raw_edges, add_edge);
+ auto add = [&](const S2Point& a, const S2Point& b) {
+ index->Add(std::make_unique<GeographyLineShape>(a, b));
+ };
+ if (edges.empty()) {
+ // A raw pole-alias line may collapse to one canonical calculation
position.
+ add(first, first);
+ }
+ std::set<S2Point, std::less<>, CustomStdAllocator<S2Point>>
interior_vertices;
+ for (const auto& [a, b] : edges) {
+ add(a, b);
+ interior_vertices.insert(a);
+ interior_vertices.insert(b);
+ }
+ for (const auto& point : interior_vertices) {
+ if (first == last || (point != first && point != last)) {
+ add(point, point);
+ }
+ }
+ }
+
+private:
+ struct Vertex {
+ double latitude;
+ S2Point point;
+ };
+ struct Circle {
+ S2Point a;
+ S2Point b;
+ S2Point normal;
+ double error;
+ DorisVector<std::pair<double, double>> ranges;
+ DorisMap<double, S2Point> vertices;
+ };
+
+ template <typename Edges, typename AddEdge>
+ void add_circle_edges(const Edges& edges, AddEdge&& add_edge) const {
+ DorisVector<Circle> circles;
+ DorisMap<int, S2PointIndex<size_t>> normal_indexes;
+ const double rounding = 32 * std::numeric_limits<double>::epsilon();
+ for (const auto& [a, b] : edges) {
+ const auto cross = a.CrossProd(b);
+ const double norm = cross.Norm();
+ // A rounded near-antipodal normal may be arbitrary: error then
covers
+ // every unoriented plane. Normals are only candidate keys, never
identity.
+ const auto normal = norm == 0 ? S2::Ortho(a) : cross / norm;
+ const double error = std::min(M_PI_2, rounding / norm);
+ size_t group = circles.size();
+ // Isolate uncertainty by exponent: one near-antipodal
representative must
+ // not widen every ordinary-circle search and recreate quadratic
work.
+ for (const auto& [exponent, index] : normal_indexes) {
+ S2ClosestPointQuery<size_t> query(&index);
+
query.mutable_options()->set_conservative_max_distance(S1Angle::Radians(
+ std::min(M_PI, error + std::ldexp(1.0, exponent) +
rounding)));
Review Comment:
[P2] Bound candidate work for distinct uncertain circles
For N spokes, use default Geography A=`LINESTRING(0 0,179.99999999999997
e1,0 0,179.99999999999997 e2,0 0,...,179.99999999999997 eN,0 0)`, where
`ei=1e-14*i/(N+1)`, and evaluate `ST_Distance(A,POINT(20 20))`. These are valid
non-antipodal retraces ([pinned S2
validation](https://github.com/google/s2geometry/blob/v0.10.0/src/s2/s2polyline.cc#L160-L182));
reversed-edge deduplication leaves N distinct planes. Every normal has
`error=pi/2` and exponent 1, so this radius becomes pi. The [closest-point
query](https://github.com/google/s2geometry/blob/v0.10.0/src/s2/s2closest_point_query.h#L32-L34)
consequently returns all prior circle keys. None can merge, forcing at least
N*(N-1)/2 full 48-term bigint/map coplanarity proofs before distance runs.
At N=10,000, only 320,025 bytes of WKB cause at least 49,995,000 proofs,
although final edges remain O(N); the later broad vertex searches add quadratic
candidate visits. This is a different instance from 4375245664's overlapping
edges on one shared circle: these genuinely different circles hit the newly
added uncertain-bucket grouping fallback. Use robust raw-coordinate circle
identities or otherwise bound this work, and add increasing-size coverage with
many uncertain circles. Counts follow statically; no benchmark was run.
--
This is an automated message from the Apache Git Service.
To respond to the message, please log on to GitHub and use the
URL above to go to the specific comment.
To unsubscribe, e-mail: [email protected]
For queries about this service, please contact Infrastructure at:
[email protected]
---------------------------------------------------------------------
To unsubscribe, e-mail: [email protected]
For additional commands, e-mail: [email protected]
