First pass at EFT-based spherical geometry accuracy (port from AccuSphGeom)

docs/userguide.rst

@@ -94,6 +94,9 @@ Supplementary Guides
 
 These user guides provide additional details about specific features in UXarray.
 
+`Accurate Spherical Geometry <user-guide/spherical-geometry-accuracy.ipynb>`_
+ How UXarray uses error-free transformations to avoid catastrophic cancellation in cross-product and point-in-polygon operations
+
 `Working with HEALPix Grids <user-guide/healpix.ipynb>`_
  Use UXarray with HEALPix
 
@@ -127,6 +130,7 @@ These user guides provide additional details about specific features in UXarray.
    user-guide/dual-mesh.ipynb
    user-guide/structured.ipynb
    user-guide/from-points.ipynb
+   user-guide/spherical-geometry-accuracy.ipynb
    user-guide/healpix.ipynb
    user-guide/holoviz.ipynb
    user-guide/from_file.ipynb

uxarray/grid/_eft.py

 def _two_sum(a, b): 
 def _two_sum(a, b): 
@@ -0,0 +1,167 @@
+"""Error-free transformations (EFT) for accurate floating-point arithmetic.
+
+In spherical-geometry computations the critical operations are cross products
+and dot products over unit vectors. When two vectors are nearly parallel, the
+difference of products that forms each cross-product component suffers
+catastrophic cancellation: both products round to the same floating-point
+value and their difference carries no significant bits. This affects
+GCA-GCA intersection of nearly tangent arcs, constant-latitude intersection
+near arc endpoints, and the ray-crossing test in point-in-polygon near polygon
+edges.
+
+The functions here represent each result as an unevaluated sum of two
+``float64`` values ``(hi, lo)`` such that ``hi + lo`` equals the
+mathematically exact result. This effectively doubles the significant bits
+available for cross-product components without resorting to arbitrary-
+precision arithmetic.
+
+These primitives are a Python/Numba port of the error-free transformation
+layer from the AccuSphGeom C++ library:
+
+    Chen, H. (2026). Accurate and Robust Algorithms for Spherical Polygon
+    Operations. EGUsphere preprint.
+    https://egusphere.copernicus.org/preprints/2026/egusphere-2026-636/
+
+    Chen, H. Accurate and Robust Great Circle Arc Intersection and Great
+    Circle Arc Constant Latitude Intersection on the Sphere. SIAM J. Sci.
+    Comput. https://doi.org/10.1137/25M1737614
+
+AccuSphGeom reference implementation (C++):
+    https://github.com/hongyuchen1030/AccuSphGeom
+
+What this module omits: AccuSphGeom's full robustness stack has three
+tiers — an EFT filter (what this module implements), Shewchuk adaptive
+predicates for results that fall inside the filter threshold, and a geogram
+exact-arithmetic fallback. This port implements only the EFT tier. For
+non-degenerate inputs in double precision this is sufficient; callers that
+need to handle geometrically degenerate inputs (coincident arcs, a query
+point exactly on a polygon edge) should add their own perturbation or
+fall-back logic.
+"""
+
+from numba import njit
+
+
+@njit(cache=True, inline="always")
+def two_sum(a, b):
+    """Knuth's TwoSum: return (s, e) with s = fl(a + b) and s + e = a + b exactly.
+
+    Floating-point addition rounds the mathematical result to the nearest
+    representable value. ``two_sum`` captures that rounding error in the
+    companion term ``e`` so that ``s + e`` equals the true sum with no
+    information lost. The cost is four extra floating-point operations beyond
+    the addition itself.
+
+    Parameters
+    ----------
+    a, b : float
+        Input values.
+
+    Returns
+    -------
+    s : float
+        Rounded sum fl(a + b).
+    e : float
+        Rounding error term; s + e = a + b exactly.
+    """
+    s = a + b
+    bp = s - a
+    e = (a - (s - bp)) + (b - bp)
+    return s, e
+
+
+@njit(cache=True, inline="always")
+def two_prod(a, b):
+    """Dekker/Veltkamp TwoProd: return (p, e) with p = fl(a * b) and p + e = a * b exactly.
+
+    Like ``two_sum`` for multiplication. Uses the Veltkamp splitting constant
+    2**27 + 1 to decompose each operand into a high and low half, then
+    reconstructs the exact rounding error from the four partial products.
+    On hardware with a fused multiply-add (FMA) instruction the error term
+    could be obtained in one step as ``fma(a, b, -p)``; the split used here
+    is portable across all Numba targets.
+
+    Parameters
+    ----------
+    a, b : float
+        Input values.
+
+    Returns
+    -------
+    p : float
+        Rounded product fl(a * b).
+    e : float
+        Rounding error term; p + e = a * b exactly.
+    """
+    p = a * b
+    factor = 134217729.0  # 2**27 + 1
+    a_hi = factor * a - (factor * a - a)
+    a_lo = a - a_hi
+    b_hi = factor * b - (factor * b - b)
+    b_lo = b - b_hi
+    e = a_lo * b_lo - (((p - a_hi * b_hi) - a_lo * b_hi) - a_hi * b_lo)
+    return p, e
+
+
+@njit(cache=True, inline="always")
+def diff_of_products(a, b, c, d):
+    """Kahan's accurate a*b - c*d using two_prod and two_sum.
+
+    Naive evaluation of ``a*b - c*d`` loses all significant bits when the two
+    products are nearly equal (catastrophic cancellation). This routine
+    computes each product exactly via ``two_prod``, subtracts the rounded
+    high parts, then folds the residual low parts back in. The result has
+    rounding error bounded by one ulp of the true value regardless of
+    cancellation.
+
+    This is the core operation that makes cross products accurate: every
+    component of ``a x b`` is a difference of two products of exactly this
+    form.
+
+    Parameters
+    ----------
+    a, b, c, d : float
+        Input scalars; computes a*b - c*d.
+
+    Returns
+    -------
+    hi : float
+        High-order part of the accurate result.
+    lo : float
+        Low-order correction term; hi + lo equals the accurate value.
+    """
+    w, e_w = two_prod(c, d)
+    x, e_x = two_prod(a, b)
+    s, e_s = two_sum(x, -w)
+    lo = (e_x - e_w) + e_s
+    return s, lo
+
+
+@njit(cache=True, inline="always")
+def accucross(a0, a1, a2, b0, b1, b2):
+    """Accurate cross product a x b returning (hi[3], lo[3]) component pairs.
+
+    Each component of a cross product is a difference of two products — the
+    exact form that ``diff_of_products`` handles. This function computes all
+    three components that way, returning six scalars such that the
+    mathematically exact cross product satisfies ``result[i] = hi[i] + lo[i]``
+    for each component. Callers that need single-precision accuracy can use
+    the hi parts alone; callers that need the full compensated result add
+    hi and lo before further use.
+
+    Parameters
+    ----------
+    a0, a1, a2 : float
+        Components of vector a.
+    b0, b1, b2 : float
+        Components of vector b.
+
+    Returns
+    -------
+    x_hi, y_hi, z_hi, x_lo, y_lo, z_lo : float
+        High and low parts of each cross-product component.
+    """
+    x_hi, x_lo = diff_of_products(a1, b2, a2, b1)
+    y_hi, y_lo = diff_of_products(a2, b0, a0, b2)
+    z_hi, z_lo = diff_of_products(a0, b1, a1, b0)
+    return x_hi, y_hi, z_hi, x_lo, y_lo, z_lo

uxarray/grid/arcs.py

@@ -4,11 +4,19 @@
 from numba import njit
 
 from uxarray.constants import ERROR_TOLERANCE, MACHINE_EPSILON
+from uxarray.grid._eft import diff_of_products, two_sum
 from uxarray.grid.coordinates import (
     _normalize_xyz_scalar,
 )
 from uxarray.grid.utils import _angle_of_2_vectors
 
+# Tolerance used to classify orient3d results as zero.  For double-precision
+# unit-vector inputs this covers rounding error in the EFT cross product.
+_PREDICATE_ZERO_TOL = 1e-15
+
+# Default tolerance for the on_minor_arc collinearity and interval tests.
+_ON_MINOR_ARC_TOL = 1e-10
+
 
 def _to_list(obj):
     if not isinstance(obj, list):
@@ -364,3 +372,103 @@ def compute_arc_length(pt_a, pt_b):
     delta_theta = np.arctan2(cross_2d, dot_2d)
 
     return rho * abs(delta_theta)
+
+
+@njit(cache=True)
+def _orient3d_on_sphere_value(a, b, q):
+    """Return the EFT-accurate value of the orient3d-on-sphere predicate.
+
+    Computes the scalar (a x b) . q using ``diff_of_products`` for the
+    cross-product components and ``two_sum`` for the final accumulation.
+    For unit vectors all coordinates are in [-1, 1], so the EFT cross product
+    provides roughly double the effective precision of a naive evaluation.
+    The result is positive when q lies to the left of the directed arc a->b,
+    negative when to the right, and near zero when q is on the great circle
+    through a and b.
+
+    Parameters
+    ----------
+    a, b, q : np.ndarray, shape (3,)
+        Unit vectors on the unit sphere.
+
+    Returns
+    -------
+    float
+        Signed determinant value.
+    """
+    x_hi, x_lo = diff_of_products(a[1], b[2], a[2], b[1])
+    y_hi, y_lo = diff_of_products(a[2], b[0], a[0], b[2])
+    z_hi, z_lo = diff_of_products(a[0], b[1], a[1], b[0])
+    p0 = (x_hi + x_lo) * q[0]
+    p1 = (y_hi + y_lo) * q[1]
+    p2 = (z_hi + z_lo) * q[2]
+    s, e = two_sum(p0, p1)
+    s, e2 = two_sum(s, p2)
+    return s + (e + e2)
+
+
+@njit(cache=True)
+def orient3d_on_sphere(a, b, q, tol=_PREDICATE_ZERO_TOL):
+    """Sign of the orient3d predicate on the unit sphere: -1, 0, or +1.
+
+    Evaluates the sign of ``(a x b) . q`` using error-free transformations to
+    avoid false zero results from floating-point cancellation near great-circle
+    boundaries. The sign determines which side of the great circle through a
+    and b the point q lies on.
+
+    Parameters
+    ----------
+    a, b, q : np.ndarray, shape (3,)
+        Unit vectors on the unit sphere.
+    tol : float, optional
+        Magnitude below which the result is classified as zero.
+
+    Returns
+    -------
+    int
+        +1 if q is to the left of a->b, -1 if to the right, 0 if collinear
+        within ``tol``.
+    """
+    v = _orient3d_on_sphere_value(a, b, q)
+    if v > tol:
+        return 1
+    if v < -tol:
+        return -1
+    return 0
+
+
+@njit(cache=True)
+def on_minor_arc(q, a, b, tol=_ON_MINOR_ARC_TOL):
+    """Return True if q lies on the minor great-circle arc from a to b.
+
+    Uses ``_orient3d_on_sphere_value`` (a compensated cross product) for the
+    collinearity test and dot products for the interval check. Compared to
+    ``point_within_gca``, this avoids the ``arctan2`` call that guards against
+    180-degree arcs and avoids the separate plane-membership check via
+    ``np.cross`` + ``np.dot``.
+
+    Parameters
+    ----------
+    q : np.ndarray, shape (3,)
+        Query point (unit vector).
+    a, b : np.ndarray, shape (3,)
+        Endpoints of the great-circle arc (unit vectors).
+    tol : float, optional
+        Tolerance for the collinearity and interval checks.
+
+    Returns
+    -------
+    bool
+        True if q lies on the minor arc ab, False otherwise.
+    """
+    # Coincident endpoints: degenerate arc, no interior.
+    if a[0] == b[0] and a[1] == b[1] and a[2] == b[2]:
+        return False
+    # Collinearity check: q must lie on the great circle through a and b.
+    if abs(_orient3d_on_sphere_value(a, b, q)) > tol:
+        return False
+    # Interval check: q must lie on the minor-arc side of both endpoints.
+    qa = a[0] * q[0] + a[1] * q[1] + a[2] * q[2]
+    qb = b[0] * q[0] + b[1] * q[1] + b[2] * q[2]
+    ab = a[0] * b[0] + a[1] * b[1] + a[2] * b[2]
+    return (qb - ab * qa) >= -tol and (qa - qb * ab) >= -tol