Fix _small_angle_of_2_vectors to evaluate Kahan's formula accurately

test/grid/geometry/test_angles.py

@@ -3,9 +3,13 @@
 """
 
 import numpy as np
+import pytest
 import xarray as xr
 
 import uxarray as ux
+from uxarray.grid.arcs import _angle_of_2_vectors
+from uxarray.grid.coordinates import _lonlat_rad_to_xyz
+from uxarray.grid.utils import _small_angle_of_2_vectors
 
 
 def test_face_node_angles_triangle():
@@ -109,3 +113,91 @@ def test_equiangle_skewness():
     assert skewness.max() < 0.15
     # skewness is well defined for all faces
     assert not np.any(skewness.isnull())
+
+
+def test_angle_of_2_vectors():
+    # Test the angle between two vectors
+    v1 = np.array([1.0, 0.0, 0.0])
+    v2 = np.array([0.0, 1.0, 0.0])
+    assert pytest.approx(_angle_of_2_vectors(v1, v2)) == np.pi / 2.0
+
+    v1 = np.array([1.0, 0.0, 0.0])
+    v2 = np.array([1.0, 0.0, 0.0])
+    assert pytest.approx(_angle_of_2_vectors(v1, v2)) == 0.0
+
+def test_angle_of_2_vectors_180_degree():
+    GCR1_cart = np.array([
+        _lonlat_rad_to_xyz(np.deg2rad(0.0),
+                                np.deg2rad(0.0)),
+        _lonlat_rad_to_xyz(np.deg2rad(181.0),
+                                np.deg2rad(0.0))
+    ])
+
+    res = _angle_of_2_vectors(GCR1_cart[0], GCR1_cart[1])
+
+    # The angle between the two vectors should be 181 degree
+    assert pytest.approx(res, abs=1e-8) == np.deg2rad(181.0)
+
+    GCR1_cart = np.array([
+        _lonlat_rad_to_xyz(np.deg2rad(170.0),
+                                np.deg2rad(89.0)),
+        _lonlat_rad_to_xyz(np.deg2rad(170.0),
+                                np.deg2rad(-10.0))
+    ])
+
+    res = _angle_of_2_vectors(GCR1_cart[0], GCR1_cart[1])
+
+    # The angle between the two vectors should be 99 degrees
+    assert pytest.approx(res, abs=1e-8) == np.deg2rad(89.0+10.0)
+
+
+def _exact_small_angle(u, v):
+    """Angle between two float vectors, for angles well below 1 rad.
+
+    The cross and dot products are formed exactly in rational arithmetic, then
+    atan(|u x v| / (u . v)) is evaluated with a 50-digit Taylor series.
+    """
+    from decimal import Decimal, localcontext
+    from fractions import Fraction
+
+    u = [Fraction(float(c)) for c in u]
+    v = [Fraction(float(c)) for c in v]
+    cross = (u[1] * v[2] - u[2] * v[1], u[2] * v[0] - u[0] * v[2], u[0] * v[1] - u[1] * v[0])
+    cross2 = sum(c * c for c in cross)
+    dot = sum(a * b for a, b in zip(u, v))
+    with localcontext() as ctx:
+        ctx.prec = 50
+        t = (Decimal(cross2.numerator) / Decimal(cross2.denominator)).sqrt() / (
+            Decimal(dot.numerator) / Decimal(dot.denominator))
+        term, total, k = t, Decimal(0), 0
+        while abs(term) > Decimal(10) ** -45:
+            total += term / (2 * k + 1) * (-1) ** k
+            term *= t * t
+            k += 1
+        return float(total)
+
+
+def test_small_angle_of_2_vectors_special_values():
+    assert _small_angle_of_2_vectors((1.0, 0.0, 0.0), (0.0, 3.0, 0.0)) == pytest.approx(np.pi / 2, rel=1e-15)
+    assert _small_angle_of_2_vectors((2.0, 0.0, 0.0), (0.5, 0.0, 0.0)) == 0.0
+    assert _small_angle_of_2_vectors((0.0, 0.0, 1.0), (0.0, 0.0, -7.0)) == pytest.approx(np.pi, rel=1e-15)
+    assert _small_angle_of_2_vectors((0.0, 0.0, 0.0), (1.0, 2.0, 3.0)) == 0.0
+
+
+@pytest.mark.parametrize("theta, rtol", [(1e-3, 1e-14), (1e-5, 1e-14), (1e-7, 1e-14), (1e-9, 1e-13)])
+def test_small_angle_of_2_vectors_accuracy(theta, rtol):
+    """Nearly parallel vectors of different, arbitrary lengths.
+
+    A naive evaluation of v|u| - u|v| loses accuracy as eps / theta (about 1e-7
+    at theta = 1e-9); with compensated products it stays within a few ulps,
+    plus a second-order (eps / theta)^2 term from the rounded norms.
+    """
+    rng = np.random.default_rng(0)
+    for _ in range(50):
+        u = rng.normal(size=3)
+        w = np.cross(u, rng.normal(size=3))
+        v = u / np.linalg.norm(u) * np.cos(theta) + w / np.linalg.norm(w) * np.sin(theta)
+        u, v = u * rng.uniform(0.1, 10.0), v * rng.uniform(0.1, 10.0)
+        expected = _exact_small_angle(u, v)
+        assert _small_angle_of_2_vectors(u, v) == pytest.approx(expected, rel=rtol, abs=0.0)
+        assert _small_angle_of_2_vectors(v, u) == pytest.approx(expected, rel=rtol, abs=0.0)

test/test_helpers.py

@@ -8,7 +8,7 @@
 from uxarray.grid.connectivity import _replace_fill_values
 from uxarray.constants import INT_DTYPE, INT_FILL_VALUE
 from uxarray.grid.coordinates import _lonlat_rad_to_xyz, _normalize_xyz, _xyz_to_lonlat_rad
-from uxarray.grid.arcs import point_within_gca, _angle_of_2_vectors, in_between
+from uxarray.grid.arcs import point_within_gca, in_between
 from uxarray.grid.utils import _get_cartesian_face_edge_nodes_array, _get_lonlat_rad_face_edge_nodes_array
 from uxarray.grid.geometry import pole_point_inside_polygon, _pole_point_inside_polygon_cartesian
 
@@ -210,41 +210,6 @@ def test_in_between():
     assert in_between(0, 1, 2)
     assert in_between(-1, -1.5, -2)
 
-def test_angle_of_2_vectors():
-    # Test the angle between two vectors
-    v1 = np.array([1.0, 0.0, 0.0])
-    v2 = np.array([0.0, 1.0, 0.0])
-    assert pytest.approx(_angle_of_2_vectors(v1, v2)) == np.pi / 2.0
-
-    v1 = np.array([1.0, 0.0, 0.0])
-    v2 = np.array([1.0, 0.0, 0.0])
-    assert pytest.approx(_angle_of_2_vectors(v1, v2)) == 0.0
-
-def test_angle_of_2_vectors_180_degree():
-    GCR1_cart = np.array([
-        _lonlat_rad_to_xyz(np.deg2rad(0.0),
-                                np.deg2rad(0.0)),
-        _lonlat_rad_to_xyz(np.deg2rad(181.0),
-                                np.deg2rad(0.0))
-    ])
-
-    res = _angle_of_2_vectors(GCR1_cart[0], GCR1_cart[1])
-
-    # The angle between the two vectors should be 181 degree
-    assert pytest.approx(res, abs=1e-8) == np.deg2rad(181.0)
-
-    GCR1_cart = np.array([
-        _lonlat_rad_to_xyz(np.deg2rad(170.0),
-                                np.deg2rad(89.0)),
-        _lonlat_rad_to_xyz(np.deg2rad(170.0),
-                                np.deg2rad(-10.0))
-    ])
-
-    res = _angle_of_2_vectors(GCR1_cart[0], GCR1_cart[1])
-
-    # The angle between the two vectors should be 99 degrees
-    assert pytest.approx(res, abs=1e-8) == np.deg2rad(89.0+10.0)
-
 def test_get_cartesian_face_edge_nodes_pipeline():
     vertices = [[0.5, 0.5, 0.5], [-0.5, 0.5, 0.5], [-0.5, -0.5, 0.5], [0.5, -0.5, 0.5]]
     vertices = [x / np.linalg.norm(x) for x in vertices]

test/utils/test_numba_math.py

@@ -16,6 +16,7 @@
     _numba_norm3,
     _numba_dot3,
     _numba_cross3,
+    _numba_diff_of_products3,
 )
 
 def test_numba_add3():
@@ -79,3 +80,16 @@ def test_numba_cross3():
     assert _numba_cross3((1, 2, 30), (4, 0, 6)) == (12, 114, -8)
     assert _numba_cross3([0,1,0],[0,0,10]) == (10,0,0)
     assert _numba_cross3(np.array([1,10,0]), np.array([0,-1,5])) == (50, -5, -1)
+
+def test_numba_diff_of_products3():
+    """ensure _numba_diff_of_products3 gives u*a - v*b accurately when the products nearly cancel."""
+    from fractions import Fraction
+    assert _numba_diff_of_products3((1.0, 2.0, 3.0), 2.0, (0.5, 1.0, 1.0), 4.0) == ((0.0, 0.0, 2.0), (0.0, 0.0, 0.0))
+    u = np.array([0.1, 0.7, -1.3])
+    v = u * (1.0 + 1e-12)       # nearly parallel, so u*a and v*b nearly cancel
+    a, b = 1.0 + 1e-12, 1.0
+    hi, lo = _numba_diff_of_products3(u, a, v, b)
+    for i in range(3):
+        exact = Fraction(u[i]) * Fraction(a) - Fraction(v[i]) * Fraction(b)
+        assert hi[i] == u[i] * a - v[i] * b                     # hi alone is the naive result
+        assert Fraction(hi[i]) + Fraction(lo[i]) == exact       # hi + lo is exact here

uxarray/grid/utils.py

@@ -6,18 +6,30 @@
 from uxarray.utils.numba_math import (
     _numba_add3,
     _numba_cross3,
+    _numba_diff_of_products3,
     _numba_dot3,
     _numba_mul3_scalar,
     _numba_norm3,
-    _numba_sub3,
 )
 from uxarray.utils.parallel import parallel_njit
 
 
 @njit(cache=True)
 def _small_angle_of_2_vectors(u, v):
     """
-    Compute the smallest angle between two vectors using the new _angle_of_2_vectors.
+    Compute the angle in [0, pi] between two vectors of arbitrary lengths.
+
+    Uses Kahan's formula
+
+        angle = 2 * atan2(| v |u| - u |v| |, | v |u| + u |v| |),
+
+    which avoids the cancellation of the arcsin/arccos forms and of the cross
+    product in atan2(|u x v|, u . v). When the vectors are nearly parallel,
+    every component of ``v |u| - u |v|`` is a difference of two nearly equal
+    products, so it is evaluated with the FMA-based ``_numba_diff_of_products3``
+    (Kahan's 2x2 determinant per component). The result is accurate to a few ulps for angles
+    down to about 1e-7 rad, for any input lengths; below that, the rounding of
+    the norms adds a second-order relative error of about (eps / angle)^2.
 
     Parameters
     ----------
@@ -29,14 +41,31 @@ def _small_angle_of_2_vectors(u, v):
     Returns
     -------
     float
-        The smallest angle between `u` and `v` in radians.
+        The angle between `u` and `v` in radians. Zero if either vector is zero.
+
+    References
+    ----------
+    Kahan, W. (2006). How Futile are Mindless Assessments of Roundoff in
+    Floating-Point Computation? https://people.eecs.berkeley.edu/~wkahan/Mindless.pdf
+
+    Chen, H., Ullrich, P. A., Panetta, J., Marsico, D., Hanke, M., Jain, R.,
+    Zhang, C., and Jacob, R. L. (2026). Accurate and robust geometric
+    algorithms for regridding on the sphere. Geoscientific Model
+    Development, 19(14), 6545-6570. https://doi.org/10.5194/gmd-19-6545-2026
     """
-    u_times_v_norm = _numba_mul3_scalar(v, _numba_norm3(u))
-    v_times_u_norm = _numba_mul3_scalar(u, _numba_norm3(v))
-    vec_minus = _numba_sub3(u_times_v_norm, v_times_u_norm)
-    vec_sum = _numba_add3(u_times_v_norm, v_times_u_norm)
-    angle_u_v_rad = 2 * np.arctan2(_numba_norm3(vec_minus), _numba_norm3(vec_sum))
-    return angle_u_v_rad
+    nu = _numba_norm3(u)
+    nv = _numba_norm3(v)
+
+    # v|u| - u|v| as a compensated difference of products: m_hi alone is the
+    # naive rounded difference, adding m_lo applies the products' rounding errors.
+    m_hi, m_lo = _numba_diff_of_products3(v, nu, u, nv)
+    m = _numba_add3(m_hi, m_lo)
+
+    # v|u| + u|v|: its cancellation near antiparallel vectors does not affect
+    # the returned angle, so it is evaluated directly.
+    s = _numba_add3(_numba_mul3_scalar(v, nu), _numba_mul3_scalar(u, nv))
+
+    return 2.0 * np.arctan2(_numba_norm3(m), _numba_norm3(s))
 
 
 @njit(cache=True)

uxarray/utils/numba_math.py

@@ -15,6 +15,8 @@
 import numpy as np
 from numba import njit
 
+from uxarray.utils.computing import diff_of_products
+
 # ------- basic arithmetic with vectors ------- #
 
 
@@ -114,6 +116,27 @@ def _numba_cross3(u, v):
     return (cx, cy, cz)
 
 
+@njit(cache=True, inline="always")
+def _numba_diff_of_products3(u, a, v, b):
+    """component-wise u * a - v * b for 3-vectors u, v and scalars a, b, without cancellation.
+
+    Each component is evaluated with Kahan's 2x2 determinant
+    (:func:`uxarray.utils.computing.diff_of_products`), so it is accurate even
+    when ``u * a`` and ``v * b`` nearly cancel.
+
+    Returns
+    -------
+    hi, lo : tuple of 3 floats
+        High and low parts of each component; ``hi[i] + lo[i]`` is the
+        accurate value. ``hi`` alone equals the naive ``u * a - v * b``, so the
+        two must be added, e.g. with ``_numba_add3(hi, lo)``.
+    """
+    h0, l0 = diff_of_products(u[0], a, v[0], b)
+    h1, l1 = diff_of_products(u[1], a, v[1], b)
+    h2, l2 = diff_of_products(u[2], a, v[2], b)
+    return (h0, h1, h2), (l0, l1, l2)
+
+
 # ------- convenience functions / helpers ------- #