First pass at EFT-based spherical geometry accuracy (port from AccuSphGeom) by rajeeja · Pull Request #1513 · UXARRAY/uxarray

rajeeja · 2026-05-22T12:00:23Z

Ports EFT primitives from AccuSphGeom (Chen 2026) into uxarray/grid/_eft.py and replaces naive cross-product arithmetic in arcs.py, intersections.py, bounds.py, and point_in_face.py with EFT-backed orient3d_on_sphere and on_minor_arc. Adds a user guide notebook covering the problem, API, and a live point-in-polygon demo.

Closes #1509

review-notebook-app · 2026-05-22T12:00:29Z

Check out this pull request on

See visual diffs & provide feedback on Jupyter Notebooks.

Powered by ReviewNB

hongyuchen1030 · 2026-05-22T19:00:36Z

We have an existed implementation for these functions in

uxarray/uxarray/utils/computing.py

Line 189 in 735a1dd

def _two_sum(a, b):

. Consider either merge them or only keep one of them

And I would suggest use other name instead of eft here, the term "EFT" specifically refers to "Error-Free" floating point operations, but the AccuCross and diff_of_productthemself is not completely Error-Free, is just more accurate than normal floating point cross (doubling the precision)

hongyuchen1030 · 2026-05-22T19:23:44Z

+
+
+@njit(cache=True)
+def _generate_lat_lon_bounds_local(face_vertices, z_min, z_max, snap_tol_deg):


For a better vectorization, consider keeping the current design in https://github.com/hongyuchen1030/AccuSphGeom/blob/56cbd6e30270ec5845a853ec72ba5b19c9128017/include/accusphgeom/algorithms/lat_lon_bounds.hpp#L107:

It uses lots of masking and avoid branching for computation steps.

hongyuchen1030 · 2026-05-22T19:25:29Z

+
+
+@njit(cache=True)
+def _generate_lat_lon_bounds_pole(face_vertices, label, z_min, z_max, snap_tol_deg):


Same here, consider keeping the structure here: https://github.com/hongyuchen1030/AccuSphGeom/blob/56cbd6e30270ec5845a853ec72ba5b19c9128017/include/accusphgeom/algorithms/lat_lon_bounds.hpp#L159

use masks instead of branching as much as possible for the vectorization purpose

hongyuchen1030 · 2026-05-22T19:26:56Z

+
+
+@njit(cache=True)
+def _face_location_info(face_vertices, polar_cap_z):


The branching here will probably slow down the vectorization, consider keeping the structure here:
https://github.com/hongyuchen1030/AccuSphGeom/blob/56cbd6e30270ec5845a853ec72ba5b19c9128017/include/accusphgeom/algorithms/lat_lon_bounds.hpp#L49

hongyuchen1030 · 2026-05-22T19:29:02Z

+@njit(cache=True)
+def _normalize_pair(x_hi, y_hi, z_hi, x_lo, y_lo, z_lo):
+    """Normalize an (hi, lo) compensated vector, returning the unit vector and magnitude."""
+    x = x_hi + x_lo


This _normalize_pair is not utilizing the EFT and it will have the same precision as the direct floating point precision. And the normalization is one of the big killer for the precision.

Also if the input are normalized, then both the gca_gca_intersection and the gca_constLat_intersection will return the normalized results already

And just in case you want an "accurately calculated norm", it is const auto sum = numeric::sum_of_squares_c<T, 3>(v.hi, v.lo); const auto norm = numeric::acc_sqrt_re(sum.hi, sum.lo);

hongyuchen1030 · 2026-05-22T19:35:16Z

@@ -301,139 +314,198 @@ def _gca_gca_intersection_cartesian(gca_a_xyz, gca_b_xyz):

 @njit(cache=True)


consider keeping the entire structure from here : https://github.com/hongyuchen1030/AccuSphGeom/blob/e4e13215dd4771b7ef01b7edf81eaf58dd6e6995/include/accusphgeom/constructions/gca_gca_intersection.hpp#L31

These EFT-fused operations are extremely sensitive for each operations and order. The accuracy can only be guaranteed iff following the exact algorithm

hongyuchen1030 · 2026-05-22T19:41:30Z



 @njit(cache=True)
 def gca_const_lat_intersection(gca_cart, const_z):


The algorithm is almost very stable around floating point precision. Within the current Uxarray Error tolerance, you almost don't have to use any other safety pre-check other than if the input are valid (The relative error bound is 3*\sqrt{1-constz^2}*machine_epsilon as long as it's constZ is at least 10^15 from the equator) . And probably the only check needed Then only check you can make is if the const latitude is at the equator . Again, make sure to follow the entire algorithm structure from below. Such precision is only guaranteed
iff it follows the exact algorithm here
https://github.com/hongyuchen1030/AccuSphGeom/blob/e4e13215dd4771b7ef01b7edf81eaf58dd6e6995/include/accusphgeom/constructions/gca_constlat_intersection.hpp#L33

But you can use the same check as in https://github.com/hongyuchen1030/AccuSphGeom/blob/e4e13215dd4771b7ef01b7edf81eaf58dd6e6995/include/accusphgeom/constructions/gca_gca_intersection.hpp#L51

to remove any invalid points

hongyuchen1030 · 2026-05-22T19:44:55Z

+    nx = nx_hi + nx_lo
+    ny = ny_hi + ny_lo
+    nz = nz_hi + nz_lo
+    denom = nx * nx + ny * ny


Denome should be calculated from
const T denom = s2.hi + s2.lo;
where
const auto s2 = numeric::sum_of_squares_c<T, 2>( {normal.hi[0], normal.hi[1]}, {normal.lo[0], normal.lo[1]});

hongyuchen1030 · 2026-05-22T19:49:10Z

@@ -301,139 +314,198 @@ def _gca_gca_intersection_cartesian(gca_a_xyz, gca_b_xyz):

 @njit(cache=True)
 def gca_gca_intersection(gca_a_xyz, gca_b_xyz):


Consider keeping the entire structure from here https://github.com/hongyuchen1030/AccuSphGeom/blob/e4e13215dd4771b7ef01b7edf81eaf58dd6e6995/include/accusphgeom/constructions/gca_gca_intersection.hpp#L31.

These algorithms are extremely sensitive to the operations and order, the accuracy is only guaranteed iff we follow the exact same algorithm described

hongyuchen1030 · 2026-05-22T19:55:30Z

-    x2_at_const_z = np.isclose(
-        x2[2], const_z, rtol=ERROR_TOLERANCE, atol=ERROR_TOLERANCE
-    )
+    # 1. Endpoint coincidence with the latitude line.


This case is still valid and calculable with the new algorithm. And we can improve the vectorization by only use mask-selection after the intersection is calculated to snap the intersection point with the endpoints

hongyuchen1030 · 2026-05-22T19:56:12Z

    z_max = extreme_gca_z(gca_cart, extreme_type="max")
-
-    # Check if the constant latitude is within the GCA range
    if not in_between(z_min, const_z, z_max):


I am not sure if this early exist will counter-effect the branching it brings to the vectorization

hongyuchen1030 · 2026-05-22T19:56:41Z

-    elif p2_intersects_gca:
-        res[0] = p2
+    # 4. Solve for the two candidate points on the latitude circle.
+    r2 = 1.0 - const_z * const_z


This is not the new EFT-fused algorithmn

rajeeja · 2026-05-22T21:19:32Z

@hongyuchen1030 Thanks for the review. The point-in-polygon predicate and orient3d sign check look good, but I can see now that the intersection functions are basically the old code with a few error-free transformation calls sprinkled in rather than a real port. I'm planning to rewrite gca_gca_intersection and gca_const_lat_intersection from scratch following your C++ implementation exactly. I will be adding the missing pieces:

compensated_dot_product
sum_of_squares_c
acc_sqrt_re
The second accucross overload that takes the hi/lo pairs

Port AccuSphGeom EFT algorithms and add spherical geometry user guide

7563aba

Closes #1509

rajeeja added 3 commits May 22, 2026 07:54

Fix pre-commit: remove unused imports and variable

16cc263

Fix pre-commit: split semicolons in notebook cells

ab3d514

Merge branch 'main' into rajeeja/accusphere

b4624fe

rajeeja requested a review from hongyuchen1030 May 22, 2026 13:24

hongyuchen1030 reviewed May 22, 2026

View reviewed changes

hongyuchen1030 requested changes May 22, 2026

View reviewed changes

rajeeja marked this pull request as draft May 22, 2026 21:18



		@njit(cache=True)
		def _generate_lat_lon_bounds_local(face_vertices, z_min, z_max, snap_tol_deg):



		@njit(cache=True)
		def _generate_lat_lon_bounds_pole(face_vertices, label, z_min, z_max, snap_tol_deg):



		@njit(cache=True)
		def _face_location_info(face_vertices, polar_cap_z):

		@@ -301,139 +314,198 @@ def _gca_gca_intersection_cartesian(gca_a_xyz, gca_b_xyz):

		@njit(cache=True)



		@njit(cache=True)
		def gca_const_lat_intersection(gca_cart, const_z):

		@@ -301,139 +314,198 @@ def _gca_gca_intersection_cartesian(gca_a_xyz, gca_b_xyz):

		@njit(cache=True)
		def gca_gca_intersection(gca_a_xyz, gca_b_xyz):

Conversation

rajeeja commented May 22, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

review-notebook-app Bot commented May 22, 2026

Uh oh!

hongyuchen1030 May 22, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

hongyuchen1030 May 22, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

hongyuchen1030 May 22, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

rajeeja commented May 22, 2026

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

2 participants

rajeeja commented May 22, 2026 •

edited

Loading

hongyuchen1030 May 22, 2026 •

edited

Loading

hongyuchen1030 May 22, 2026 •

edited

Loading

hongyuchen1030 May 22, 2026 •

edited

Loading