EFT-based compensated arithmetic for spherical geometry (GCA intersection, point-in-face, face bounds)

[rajeeja]

Closes #1509

Ports the compensated-arithmetic tier from AccuSphGeom into UXarray's Numba spherical geometry stack, replacing cancellation-prone naive cross-product paths with compensated primitives throughout arc intersection, point-in-face, and face bounds.

What changed

• utils/computing.py — EFT and compensated primitives: two_sum, two_prod, diff_of_products, accucross, accucross_pair, acc_sqrt_re, and fixed-size compensated dot/sum-of-squares helpers. Breaking change: the old cross_fma and dot_fma functions (which depended on the optional pyfma package) are removed. Any code importing those names directly will get an ImportError; use accucross and the _cdp* helpers instead.
• grid/arcs.py — new _orient3d_on_sphere_value, orient3d_on_sphere, on_minor_arc predicates using compensated arithmetic.
• grid/intersections.py — gca_gca_intersection and gca_const_lat_intersection rewritten as three-layer stacks: pure numerical kernel (L1), integer-mask status layer (L2), UXarray dispatcher (L3). Removed dead _gca_gca_intersection_cartesian shim.
• grid/point_in_face.py — _face_contains_point delegates to a new ray-casting SPIP implementation (_point_in_polygon_sphere) using orient3d_on_sphere instead of the old winding-number via arctan2.
• grid/bounds.py — new _construct_face_bounds_array_gca path for pure-GCA grids, computing interior arc z-extrema correctly via the compensated kernel. Old path retained for latlon/mixed-edge grids.

What this is not

Not a full port of AccuSphGeom's robustness stack. Excluded: Shewchuk adaptive predicates, Simulation of Simplicity (requires per-vertex global IDs not available in UXarray's polygon representation), geogram exact fallback. This implements only the compensated-arithmetic tier — roughly twice as accurate as naive floating-point on near-tangent cases, with no measurable runtime overhead on intersection kernels.

Performance

Operation	vs. naive
gca_gca_intersection	~0% overhead (≈1 µs/call)
gca_const_lat_intersection	~0% overhead (≈1 µs/call)
Grid.bounds (face bounds, cached once)	~40% slower per face
Cross-section / zonal mean	~5% overhead

Accuracy

Measured on the AccuSphGeom baseline suite (31 near-tangent GCA-GCA pairs, angles down to 10⁻⁵°):

Arc angle	Naive error	EFT error	Improvement
≈ 0.000001°	1.07 × 10⁻⁷	1.24 × 10⁻¹⁶	861 million×
≈ 0.001°	5.04 × 10⁻¹²	1.14 × 10⁻¹⁶	44 000×
≈ 5°	8.47 × 10⁻¹⁶	1.58 × 10⁻¹⁶	5×

Tests

241 new baseline regression tests ported from the AccuSphGeom C++ suite covering near-tangent GCA-GCA pairs, GCA/constant-latitude cases, and spherical point-in-polygon.

[review-notebook-app]

Check out this pull request on 

See visual diffs & provide feedback on Jupyter Notebooks.

---

Powered by ReviewNB

[hongyuchen1030]

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)

[rajeeja]

Done: merged this into uxarray/utils/computing.py and renamed the docs/API section to compensated arithmetic. I only call two_sum/two_prod EFT now; the cross-product helpers are described as compensated algorithms.

[hongyuchen1030]

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.

[rajeeja]

Done: rewrote this path closer to the AccuSphGeom structure. The local bounds computation now uses mask-selection for extrema/snap decisions instead of branching through separate cases.

[hongyuchen1030]

The implementation for this entire function involves lots of new branching and operations that are not existed in the AccuSphGeom

[rajeeja]

The remaining branching in _generate_lat_lon_bounds_local and _generate_lat_lon_bounds_pole is UXarray-specific post-kernel formatting, not part of the EFT computation path. These functions convert the z-extrema from _face_location_info into degree lat/lon bounds using UXarray's encoding conventions (antimeridian signaled by lon_min > lon_max, pole-enclosing faces returned as 0/360). AccuSphGeom does not have an equivalent layer since it uses different output conventions. Both functions are now documented as UXarray-specific formatting steps, clearly separated from the EFT kernel.

[hongyuchen1030]

See the later comments: Uxarray actually uses the same criteria as the AccuSphGeom, which is also the common best-practice criteria: 1. the input faces are always counterclockwise 2. The GCA will always pick the minor arc that's less than 180 degree. A new issue should be opened to address the convention cleanup. For detailed definition, please refer to the https://egusphere.copernicus.org/preprints/2026/egusphere-2026-636/egusphere-2026-636.pdf, Section 2: Terminology and Problem Definition about node, edge, face

[hongyuchen1030]

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

[rajeeja]

Done: same cleanup here. The polar bounds path now keeps the AccuSphGeom-style local/polar structure and uses mask-selection for the snap logic.

[hongyuchen1030]

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

[rajeeja]

Done: _face_location_info now follows the AccuSphGeom face-location flow. The per-edge extrema use mask-selection; only the final face classification branches remain.

[hongyuchen1030]

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.

[hongyuchen1030]

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

[hongyuchen1030]

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);

[rajeeja]

Done: removed _normalize_pair. The intersection code now keeps compensated normals/intersection vectors, and the baseline near-tangent cases pass.

[rajeeja]

Done. _normalize_pair is removed. For _accux_gca, the direction vector is normalized using plain math.sqrt on the collapsed scalars rather than sum_of_squares_c + acc_sqrt_re — we found that passing collapsed scalars as hi/lo pairs to _sum_sq_c3 misrepresents their error structure and actually degrades accuracy on the baseline suite (pair_id=8 goes from <1e-15 to ~7e-10 error). The collapsed norm is sufficient because the unit-sphere error is dominated by the accucross_pair step, not the normalization. This is documented in a comment in _accux_gca.

[hongyuchen1030]

Thanks for checking this. I do not think plain math.sqrt should be treated as more accurate here. In general floating-point computations, the square-root / normalization step is one of the major contributors to rounding error, which is exactly why the EFT-style accurate square-root path is used.

The degradation you described for the EFT-style square-root path sounds more like an implementation bug to me, especially if collapsed scalars are being repackaged as hi/lo pairs, or if sum_of_squares_c and acc_sqrt_re are being combined in a way that does not match the original algorithm.

Please refer to the following resources to double-check whether the EFT accurate square root is implemented correctly:

from https://egusphere.copernicus.org/preprints/2026/egusphere-2026-636/egusphere-2026-636.pdf

and

from https://link.springer.com/article/10.1007/s13160-023-00593-8

The accurate square-root routine should use a Taylor-expansion correction together with the normal floating-point sqrt operation to obtain both the rounded square-root value and its rounding-error term.

One important detail is that, for the Taylor expansion to work correctly, the low part passed into acc_sqrt_re must be much smaller than the hi part. We need to make sure the hi/lo input is normalized first, for example by enforcing:

new_hi, new_lo = FastTwoSum(original_hi, original_low)

before feeding the hi/lo pair into acc_sqrt_re.

If the UXarray version produces worse accuracy than plain math.sqrt, we should first isolate that path and verify the implementation before documenting plain math.sqrt as the intended choice.

[hongyuchen1030]

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

[rajeeja]

Done: rewrote GCA-GCA around the AccuSphGeom structure: accucross normals, accucross_pair for the intersection direction, then candidate filtering with on_minor_arc.

[hongyuchen1030]

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

[hongyuchen1030]

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

[rajeeja]

Done: rewrote const-lat around the AccuSphGeom accux_constlat flow and filter invalid candidates after computing them.

[hongyuchen1030]

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]});

[rajeeja]

Done: denom now comes from _sum_sq_c2(...) as s2_hi + s2_lo.

[hongyuchen1030]

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

[rajeeja]

Done: removed the endpoint early return. The code computes candidates first, then snaps near-endpoint results afterward.

[hongyuchen1030]

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

[rajeeja]

Done: removed the endpoint pre-exits. I kept only invalid denom / negative-discriminant exits.

[hongyuchen1030]

This is not the new EFT-fused algorithmn

[rajeeja]

Done: replaced this with the compensated s2/s3, two_prod, cdp4, and acc_sqrt_re flow from AccuSphGeom.

[rajeeja]

@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

[hongyuchen1030]

We're technically still unable to claim "For non-degenerate inputs in double precision this is sufficient". We can claim "This is twice more accurate than the direct floating point implementation without computation overhead with vectorization and parralization"

[hongyuchen1030]

If we don't use any adaptive arithmetic, we cannot claim for robustness here. We can only claim our point-in-face use the new algorithm that are more accurate (if we use any EFT operations)

[rajeeja]

Fixed in commit 002b9fe. Removed all claims of 'sufficient', 'non-degenerate inputs', and 'robustness'. The module docstring now says the compensated routines are roughly twice as accurate as direct floating-point equivalents; robustness against all degenerate inputs would require adding an adaptive predicate or exact-arithmetic fallback tier. Same language used for the point-in-face docstring.

[hongyuchen1030]

These edge cases should be included in the test cases already. These scenarios don't have more risks than the "normal input" here. (They probably have the same risk since the intersection point is a newly constructed point)

[rajeeja]

Agreed. The denom==0 and planar_sq<0 edge cases propagate as inf/NaN through the kernel and are correctly handled by the isfinite mask in the status layer — no special-casing needed. The guards were removed from the kernel in commit 56c6ce1; only the status layer classifies them.

[hongyuchen1030]

The a = min(max(a_raw, 0.0), 1.0) should be able to prevent the divide by zero here. But we can keep it just in case

[rajeeja]

Kept as-is. The denom==0 guard sets a_raw=-1 so the clamp produces a=0 (use endpoint), which is safe. Matches your comment — kept just in case, and the guard costs nothing.

[hongyuchen1030]

The guard is still a branching inside computation hot path, just because it looks condensed into one line, it's still semantically equivalent to a normal if-else branch

[hongyuchen1030]

Consider utilizing the boolean operation from AccuSphHGeom to reduce branching

  const MaskType<T> north_pole_candidate_mask = (face_z_max >= polar_cap_z);
  const MaskType<T> south_pole_candidate_mask = (face_z_min <= -polar_cap_z);
  const MaskType<T> local_mask = !(north_pole_candidate_mask | south_pole_candidate_mask);

[rajeeja]

Done in commit 002b9fe. The label computation now uses pure integer multiplication instead of boolean branching: label = north_pole_candidate * _FACE_LOC_NORTH_POLAR + (1 - north_pole_candidate) * south_pole_candidate * _FACE_LOC_SOUTH_POLAR. The unused local variable was also removed.

[hongyuchen1030]

We probably don't need this work around anymore with the current latlon bounds implementation

[rajeeja]

This function is still needed and cannot be removed. UXarray uses a lon_min > lon_max encoding to signal antimeridian-crossing faces throughout the bounds and cross-section APIs — AccuSphGeom uses a union-of-intervals convention instead. The largest-gap algorithm translates between them. Removed the snap logic (which you flagged separately) but the antimeridian detection itself must stay. Added a docstring explaining this.

[hongyuchen1030]

Thank you for the clarification. I understand why this function is still being retained for the current UXarray lon_min > lon_max antimeridian convention; however, I believe we should be cautious about relying on it in this PR.

This function was originally referenced from TempestExtremes, which was designed around vertex point clouds rather than fully connected face geometry. It was introduced early in UXarray’s development as a placeholder when the framework was first implemented.

For spherical face geometry, UXarray has generally relied on two assumptions to define the intended face and edge interpretation:

1. Face vertices follow a consistent orientation, typically counterclockwise.
2. Each edge is interpreted as a minor arc, meaning the arc interval is less than 180 degrees, with near-180-degree ambiguous cases either rejected or handled explicitly.

Under this convention, separate longitude or antimeridian inference should not be necessary for the AccuX-style geometric path, and relying on it may introduce inconsistencies with the intended spherical geometry interpretation.

Additionally, I do not believe AccuSphGeom represents this using a union-of-intervals convention as described here. Instead, it assumes consistent face orientation and interprets edges as minor arcs.

Therefore, while this function can remain for now if it is still required by existing bounds or cross-section APIs, we should avoid making this PR dependent on it for the AccuX path. I recommend opening a separate issue to more thoroughly address the longitude-bound and antimeridian convention, as it has broader implications for UXarray’s geometry assumptions.

[hongyuchen1030]

The degeneracy check like this will greatly impact the vectorization performance, any check like this should be isolated from the intersection computation kernel. And the AccuXGCA itself is able to handle extremely short arcs

[rajeeja]

Done. Degeneracy checks (denom==0, planar_sq<0) are handled by letting NaN/inf propagate through the kernel — the status layer uses isfinite masks to classify them without any branching in the hot path. The _accux_gca and _accux_constlat kernels contain no if/else guards; only the 1.0/vn if vn!=0.0 else np.inf guard remains, which produces inf so the isfinite mask rejects the candidate without branching.

[hongyuchen1030]

the remaining expression

1.0 / vn if vn != 0.0 else np.inf

It is semantically equivalent to:

if vn != 0.0:
    inv_vn = 1.0 / vn
else:
    inv_vn = np.inf

is still a conditional guard in the hot path. Even though it is written as a Python conditional expression rather than a multiline if/else, it is semantically still equivalent to an if vn != 0.0 branch. Numba may lower it to a branch or a conditional select depending on the generated code, but either way it is still a predicate inside the kernel.

We should be able to handle it like the Degeneracy check and let error/inf be caught by the later try-catch

[hongyuchen1030]

Again, all these if-else branch will impact the vectorization behavior, that's why the original try_gca_gca_intersection only use mask/boolean operations. The point is, we do not try to prevent the NaN or Divide by zero in the intersection kernel, these errors will be recorded in the "status" so an outside dispatch function will know how to handle each case.

[rajeeja]

Done in commits 51c00c2 and 002b9fe. Both _try_gca_gca_intersection and _try_gca_const_lat_intersection now use integer mask arithmetic throughout — no if/else in the hot path. pos_fin = int(isfinite(...)), validity via pos_valid = pos_fin * pos_on_a * pos_on_b, point selection via pos_mask * pos + neg_mask * neg, status via both + none * 2. The only remaining guards wrap on_minor_arc calls to prevent calling it with inf inputs, which is unavoidable.

[hongyuchen1030]

Why it's still unavoidable?

the potential `inf` propagation from the `on_minor_arc` should be handled through the `_finite` mask already

[hongyuchen1030]

Isolate these if-else branch outside of the computation kernel

[rajeeja]

Done. All branching is now in L3 (gca_gca_intersection / gca_const_lat_intersection) — the dispatcher layer. L1 kernels (_accux_gca, _accux_constlat) are branch-free. L2 (_try_gca_gca_intersection, _try_gca_const_lat_intersection) uses integer mask arithmetic with no if/else in the computational path.

[hongyuchen1030]

Snapping involving branching should be isolated outside of the computation kernels. The intersection computation kernels should not include any branching and always stay in the SIMD-packed vectorization friendly form

[rajeeja]

Done. _snap_const_lat_endpoint lives entirely in L3 (gca_const_lat_intersection) — it is called after the kernel and status layer have completed, never inside _accux_constlat or _try_gca_const_lat_intersection. The L1 and L2 layers remain branch-free and kernel-pure.

[hongyuchen1030]

The implementation in

include/accusphgeom/constructions/gca_constlat_intersection.hpp

was intentionally designed with three separate layers, and I think it is important to keep these layers conceptually and structurally separated.

1. Core computation kernel: accux_constlat

   This is the numerical kernel. It contains the carefully designed AccuSphere algorithm for computing the great-circle-arc and constant-latitude intersection.

   This layer should stay isolated and intact. It is designed to be SIMD-packing friendly, branch-free in the hot path, and easy to reason about for accuracy, performance, and future optimization. This function should not be mixed with high-level filtering, status interpretation, or UXarray-specific logic.

2. SIMD-friendly batch API: try_gca_constlat_intersection

   This is the API intended for heavy computation.

   It is still SIMD-packing friendly and is designed to compute the full matrix of possible intersection results from the input faces or edge lists. It does not immediately discard invalid intersections. Instead, it returns the computed candidate points together with a status flag indicating whether each result is valid.

   This is the layer we want to use for large-scale vectorized computation, because it keeps the computation uniform. Even invalid candidates are part of the output matrix, and validity is represented by status instead of control flow.

3. Dispatcher / convenience API: gca_constlat_intersection

   This is the higher-level user-facing dispatcher.

   This layer is allowed to branch. It reads the status returned by try_gca_constlat_intersection, filters out invalid results, and returns only the valid intersection points. This is useful as a lightweight convenience API, especially when the caller wants clean geometry results instead of the full vectorized computation matrix.

The important point is that these three layers serve different purposes.

The core kernel should focus only on accurate numerical computation. The batch API should focus on uniform, SIMD-friendly execution. The dispatcher should handle branching, filtering, and convenience behavior.

Mixing these layers together is not ideal because it makes the heavy computation path harder to vectorize, harder to optimize, and harder to verify numerically. Once filtering, branching, and application-specific logic are pushed into the core kernel or SIMD batch layer, the implementation becomes less predictable and less suitable for large-scale computation.

For this kind of heavy numerical geometry kernel, the best practice is usually:

• keep the low-level numerical kernel pure, isolated, and branch-minimized;
• keep the batch/vectorized API uniform and status-based;
• move branching, filtering, and user-facing convenience behavior to a separate dispatcher layer;
• avoid mixing UXarray-specific data handling with the core computational algorithm.

So for UXarray integration, I think the preferred workflow should be:

1. Use try_gca_constlat_intersection for the main large-scale computation.
2. Preserve the full output matrix and status information during the vectorized computation stage.
3. Apply filtering only afterward, either through gca_constlat_intersection or through UXarray-side post-processing logic.

That way, we preserve the original AccuSphere design: accurate computation first, SIMD-friendly batch execution second, and lightweight branching/filtering only at the outer layer.

[rajeeja]

Understood — the key point is the design separation, not the C++ specifics. We've now split both intersection functions into the three layers you described:

• _accux_constlat / _accux_gca — pure numerical kernels, no branching, no validity filtering, compensated operations in the exact AccuSphGeom order
• _try_gca_const_lat_intersection / _try_gca_gca_intersection — batch/status layer using integer mask arithmetic (pos_valid * (1 - neg_valid)) to select candidates without branching in the hot path; returns (point, status, pos, neg)
• gca_const_lat_intersection / gca_gca_intersection — outer dispatcher that reads status, applies endpoint snapping, and packages UXarray's NaN-filled result format; all UXarray-specific branching lives here

[rajeeja]

@hongyuchen1030 At the bare-kernel level — the analog of your benchmark_gca_constlat_SIMDPack.cpp (kernels L232/L247, thread loop omp_set_num_threads L361, median L431) — the AccuX/FP64 ratio shows the decreasing trend you expected, matching the shape of your Fig 6:

threads	fp64 ns	accux ns	ratio
1	2.67	12.86	4.82
2	1.37	6.54	4.76
4	0.88	3.25	3.67
8	0.64	1.84	2.85

At equal thread count this tracks your numbers closely — your Fig 6 is 2.38 at 8 threads, ours is 2.85. Your ratio only reaches ~1.4 at 16 threads, which this machine can't run: it is an M1 Max with 8 performance cores (the 2 efficiency cores distort the loop, so I cap at 8). N=32M to saturate, median (not best-of-N).

Could you run our Python bench on the same architecture you used for Fig 6? Then we would have a true same-hardware comparison out to 16 threads and can confirm the ratio converges the same way. Code: python/bench_threads_py.py (you have access).

The earlier flat plot was the full gca_const_lat_intersection dispatcher (benchmarks/thread_scaling_constlat.py), where shared plumbing — array allocation, validity masks, endpoint snapping — dilutes the ratio. The kernel comparison above is the apples-to-apples match to your figure.

[rajeeja]

Re-ran at your 2M working set (kDefaultDataSize, L33); code: python/bench_threads_py.py.

threads	fp64 ns	accux ns	ratio
1	2.70	12.87	4.77
2	1.42	6.47	4.55
4	0.82	3.30	4.04
8	0.80	2.86	3.58

At 2M our FP64 flatlines 4→8 (0.82 → 0.80), same as your Fig 6 (1.03/1.04) — both hit the bandwidth wall. Ratio decreases; caps at 8 (M1 Max, 8 P-cores). My earlier 32M numbers hid this since the larger set kept FP64 scaling.

[cmdupuis3]

Hmm, I think it's certainly worth exploring more. It seems like the scaling behavior just isn't as aggressive in Python yet.

I think though that we should be looking to wind this PR down if there are no further accuracy issues, and we can continue working on the scaling in a new issue.

@hongyuchen1030 Does that sound okay with you?

[hongyuchen1030]

Hmm, I think it's certainly worth exploring more. It seems like the scaling behavior just isn't as aggressive in Python yet.

I think though that we should be looking to wind this PR down if there are no further accuracy issues, and we can continue working on the scaling in a new issue.

@hongyuchen1030 Does that sound okay with you?

@cmdupuis3,

Moving the Python scaling issue to a separate PR sounds very reasonable.

For the current work, I suggest splitting the remaining changes into two additional PRs:

1. A PR for the new AccuXConstLat and AccuXGCA implementation.
2. A PR for the latitude–longitude bounds related generation.

For the new AccuXConstLat and AccuXGCA implementation, the current benchmarks are still not sufficient to determine whether the port is fully correct. Given the complexity of the floating-point operations, even the existing accuracy test cases should be treated only as a preliminary validation.

We have also already identified several mismatches between the Python implementation and the AccuSphGeom implementation that caused problems. Since the computational kernels themselves are relatively small, I think a manual line-by-line comparison is appropriate here.

@rajeeja and @cmdupuis3, could both of you please perform an independent manual line-by-line review of the following implementations? Thank you.

I would prefer that this review be performed manually rather than delegated to an LLM. These kernels are explicitly structured around floating-point evaluation order, and LLMs are generally unreliable when reviewing this type of carefully designed numerical code.

1. Low-level EFT kernels

UXarray implementation:

https://github.com/UXARRAY/uxarray/blob/rajeeja/accusphere/uxarray/utils/computing.py

AccuSphGeom reference implementation:

https://github.com/hongyuchen1030/AccuSphGeom/blob/main/include/accusphgeom/numeric/eft.hpp

The non-FMA versions of two_sum and two_prod are provided here:

two_sum:

https://github.com/hongyuchen1030/regridding-geom-benchmark/blob/4081171248040c60ea72484d92ed3f037a08df52/regrid_benchmark_utils.h#L75

two_prod:

https://www.tuhh.de/ti3/paper/rump/OgRuOi05.pdf

2. GCA × ConstLat computation

2.1 Low-level computation kernel

No branching should be introduced in this kernel.

UXarray implementation:

uxarray/uxarray/grid/intersections.py

Line 474 in a7d8e88

def _accux_constlat_scalar(a0, a1, a2, b0, b1, b2, const_z):

AccuSphGeom reference implementation:

https://github.com/hongyuchen1030/AccuSphGeom/blob/90b76b2e68220c15dbe1c85c00b92b385da19624/include/accusphgeom/constructions/gca_constlat_intersection.hpp#L33

2.2 Batch-processing API

The batch-processing path should also avoid branching. Boolean masking should be sufficient.

UXarray implementation:

uxarray/uxarray/grid/intersections.py

Line 557 in a7d8e88

def _try_gca_const_lat_intersection(gca_cart, const_z):

AccuSphGeom reference implementation:

https://github.com/hongyuchen1030/AccuSphGeom/blob/90b76b2e68220c15dbe1c85c00b92b385da19624/include/accusphgeom/constructions/gca_constlat_intersection.hpp#L79

3. GCA × GCA computation

3.1 Low-level AccuXGCA kernel

Similarly, no branching should be introduced in this kernel.

UXarray implementation:

uxarray/uxarray/grid/intersections.py

Line 305 in a7d8e88

def _accux_gca(w0, w1, v0, v1):

AccuSphGeom reference implementation:

https://github.com/hongyuchen1030/AccuSphGeom/blob/90b76b2e68220c15dbe1c85c00b92b385da19624/include/accusphgeom/constructions/gca_gca_intersection.hpp#L31

3.2 Batch-processing API

UXarray implementation:

uxarray/uxarray/grid/intersections.py

Line 364 in a7d8e88

def _try_gca_gca_intersection(w0, w1, v0, v1):

AccuSphGeom reference implementation:

https://github.com/hongyuchen1030/AccuSphGeom/blob/90b76b2e68220c15dbe1c85c00b92b385da19624/include/accusphgeom/constructions/gca_gca_intersection.hpp#L51

Once the manual comparisons are complete and the implementations are confirmed to match, I will be happy to perform a final review and merge the PR.

For the point-in-face and latitude–longitude bounds generation work, I noticed several major algorithmic mismatches. Many of them appear to be related to legacy code that should be cleaned up first.

For that reason, I think moving the point-in-face and bounds-generation changes into a separate PR would be more manageable and would keep the scope of the current PR focused.

[rajeeja]

Agreed — the split is the high-value move. I'll keep this PR focused on the AccuXConstLat / AccuXGCA kernels and move point-in-face and lat-lon bounds into a separate PR. Since those depend on the EFT kernels here, I'll base the new PR on this branch for now and retarget it to main once this lands. I'll go through the kernel comparisons against eft.hpp and the AccuSphGeom constructions, and hopefully @cmdupuis3 can verify independently too. Will open a separate issue for the Python thread-scaling.

[hongyuchen1030]

Agreed — the split is the high-value move. I'll keep this PR focused on the AccuXConstLat / AccuXGCA kernels and move point-in-face and lat-lon bounds into a separate PR. Since those depend on the EFT kernels here, I'll base the new PR on this branch for now and retarget it to main once this lands. I'll go through the kernel comparisons against eft.hpp and the AccuSphGeom constructions, and hopefully @cmdupuis3 can verify independently too. Will open a separate issue for the Python thread-scaling.

Thanks, and for the manual verification, if you or @cmdupuis3 spot ANY mismatch between its C++ implementation, please label it out so we can discuss it together. In my quick checking, I already spot several new degeneracy handling rules were added in this python implementation and introduced branching in the computation hot-path.

In this way we can have a clear understanding for the potential scaling problems reasoning for the future implementation.

[rajeeja]

I fixed the degeneracy issues you flagged. This should be ready for another look whenever you get a chance — thanks!

[hongyuchen1030]

I fixed the degeneracy issues you flagged. This should be ready for another look whenever you get a chance — thanks!

Thanks for addressing those issues. Just to double-check, have you completed a line-by-line comparison between the C++ and Python implementations to confirm that the computational logic is identical throughout and to explicitly identify any remaining mismatches?

[rajeeja]

Ok, went again through both kernels line by line. Fixed one real mismatch, found two intentional deviations I don't think we should change.

Fixed: diff_of_products (computing.py:208-212) had a different addition grouping than eft.hpp's accurate_difference_of_products (56-63). Reordered to match. Pushed.

Deviation 1 — reciprocal-multiply vs division (intersections.py:354-358, 514-518 vs gca_gca_intersection.hpp:42-44, gca_constlat_intersection.hpp:70-73): we precompute 1/denom once and multiply; C++ divides per component. These can differ by 1 ULP (~25% of the time), but that's 4-20x smaller than our existing test tolerances, and switching costs ~15-17% in the hot path. Keeping it as is.

Deviation 2 — degeneracy check (arcs.py:512-516 vs on_minor_arc.hpp:53-56 / mask.hpp): we use |a×b|² < 1e-30 instead of exact equality. Already landed (793f351, tests in f78f5ff) — exact equality misses near-duplicate endpoints from independent trig paths (e.g. lon vs lon+2π), and when that happens the predicate degrades to always-true for any query point on the sphere. C++ has a Shewchuk fallback tier to catch this that we don't port, so the wider threshold is standing in for that. False-positive floor is ~1e-13deg, way below any real mesh edge.

@hongyuchen1030 let me know if you want to dig into either of these further. Also fine to hold off — @cmdupuis3 said he'd take a look tomorrow, might make sense to wait and review both passes together.

[hongyuchen1030]

Ok, went again through both kernels line by line. Fixed one real mismatch, found two intentional deviations I don't think we should change.

Fixed: diff_of_products (computing.py:208-212) had a different addition grouping than eft.hpp's accurate_difference_of_products (56-63). Reordered to match. Pushed.

Deviation 1 — reciprocal-multiply vs division (intersections.py:354-358, 514-518 vs gca_gca_intersection.hpp:42-44, gca_constlat_intersection.hpp:70-73): we precompute 1/denom once and multiply; C++ divides per component. These can differ by 1 ULP (~25% of the time), but that's 4-20x smaller than our existing test tolerances, and switching costs ~15-17% in the hot path. Keeping it as is.

Deviation 2 — degeneracy check (arcs.py:512-516 vs on_minor_arc.hpp:53-56 / mask.hpp): we use |a×b|² < 1e-30 instead of exact equality. Already landed (793f351, tests in f78f5ff) — exact equality misses near-duplicate endpoints from independent trig paths (e.g. lon vs lon+2π), and when that happens the predicate degrades to always-true for any query point on the sphere. C++ has a Shewchuk fallback tier to catch this that we don't port, so the wider threshold is standing in for that. False-positive floor is ~1e-13deg, way below any real mesh edge.

@hongyuchen1030 let me know if you want to dig into either of these further. Also fine to hold off — @cmdupuis3 said he'd take a look tomorrow, might make sense to wait and review both passes together.

I took a quick look and found an unaddressed implementation mismatch:

uxarray/uxarray/grid/intersections.py

Line 351 in a7d8e88

inv = 1.0 / vn if vn != 0.0 else np.inf

This branching does not exist in the C++ implementation, and it should not be present here. Introducing this type of branching blocks vectorization and changes the computational structure of the implementation. It should have been identified and explicitly flagged during the requested line-by-line comparison.

This is precisely why I was concerned about silent implementation mismatches and requested a manual line-by-line review without relying on LLM-generated comparisons.

[rajeeja]

Hey, not all LLM — used tooling to cross-check against the C++ reference. Both spots you'd find are from the pinned commit; already fixed in 4af6ed1a on 7/18 (_accux_gca's inv, _accux_constlat_scalar's inv_denom, plus the mask logic around them — all branch-free now). Please pull latest (1d00adbd) before re-checking.

[cmdupuis3]

I can't find anything else that would be substantially wrong.

[hongyuchen1030]

Thanks for the great work, @rajeeja. I've flagged a few minor changes — mostly comment fixes, plus one small implementation note. Once those are addressed, I'm happy to approve the merge.

[hongyuchen1030]

Thanks for the work.

Before reviewing, I want to share an API-level benchmark setup I recently implemented for the AccuSphGeom package. I think it can serve as a useful reference for structuring this kind of benchmark in either Python or C++.

The benchmark covers both GCA x GCA and GCA x ConstLat intersection APIs. It compares the naive floating-point API against the corresponding AccuX API under scalar, SIMD-vectorized, and OpenMP-parallel execution. I also wrote a detailed README explaining the motivation, scientific reasoning, benchmark design, and expected interpretation of the results.

The README includes a checklist for verifying that the benchmark is moving in the correct direction. I am copying it here for convenience. Once you go through this checklist and verify the current implementation against it, I will proceed with reviewing the changes.

Checklist for reproducing the benchmarks

This checklist is intended to help users reproduce the benchmark on their own machines. This performance test is non-trivial, delicate, and often confusing for users who are not already familiar with SIMD and HPC performance benchmarking. The goal is to make the setup easier to follow, verify, and reproduce. The bold items are the most important sanity checks.

• Build the project in Release mode.
• Confirm that the benchmark is compiled with -O3 -march=native -ffp-contract=off -fno-fast-math.
• Confirm that OpenMP is enabled for the benchmark target.
• Set up the naive floating-point APIs for the API level you want to benchmark.
• Set up the corresponding accux_* API.
• Make sure both API families return the same full API output shape.
• Make sure the scalar path uses direct double, not EigenPack<1>.
• Make sure the SIMD path uses the same pack type for both the naive floating-point implementation and the AccuX implementation.
• Before using the real AccuX arithmetic, temporarily replace only the implementation body inside the accux_* API with the same implementation used in the naive floating-point API. Keep the API entry point, return type, and output storage unchanged.
• Run the benchmark between the naive floating-point API and this temporary same-implementation accux_* API.
• Confirm that the two API entry points produce highly similar, and ideally nearly identical, performance.
• If the performance is not similar at this stage, inspect the engineering design carefully before continuing.
• Common things to inspect include inlining visibility, translation-unit separation, output storage, scalar-vs-pack dispatch, OpenMP run order, and repeated-trial stability.
• Restore the real AccuX implementation.
• Measure the saturation data size for your system.
• Update either DATA_SIZE in gca_constLat/run_gca_constlat_SIMDPack_benchmark.sh or kDefaultDataSize in gca_constLat/benchmark_gca_constlat_SIMDPack.cpp to match your measured saturation size.
• Rerun the benchmark using the saturation data size.
• Confirm that the naive floating-point API performance is stable at the chosen data size.
• Confirm that the AccuX API performance is also stable at the chosen data size.
• Confirm that, with the optimized setup and saturation data size, the accux_* APIs show no meaningful computational overhead relative to the naive floating-point APIs.

[hongyuchen1030]

Thank you for the clarification. I understand why this function is still being retained for the current UXarray lon_min > lon_max antimeridian convention; however, I believe we should be cautious about relying on it in this PR.

This function was originally referenced from TempestExtremes, which was designed around vertex point clouds rather than fully connected face geometry. It was introduced early in UXarray’s development as a placeholder when the framework was first implemented.

For spherical face geometry, UXarray has generally relied on two assumptions to define the intended face and edge interpretation:

1. Face vertices follow a consistent orientation, typically counterclockwise.
2. Each edge is interpreted as a minor arc, meaning the arc interval is less than 180 degrees, with near-180-degree ambiguous cases either rejected or handled explicitly.

Under this convention, separate longitude or antimeridian inference should not be necessary for the AccuX-style geometric path, and relying on it may introduce inconsistencies with the intended spherical geometry interpretation.

Additionally, I do not believe AccuSphGeom represents this using a union-of-intervals convention as described here. Instead, it assumes consistent face orientation and interprets edges as minor arcs.

Therefore, while this function can remain for now if it is still required by existing bounds or cross-section APIs, we should avoid making this PR dependent on it for the AccuX path. I recommend opening a separate issue to more thoroughly address the longitude-bound and antimeridian convention, as it has broader implications for UXarray’s geometry assumptions.

[hongyuchen1030]

the remaining expression

1.0 / vn if vn != 0.0 else np.inf

It is semantically equivalent to:

if vn != 0.0:
    inv_vn = 1.0 / vn
else:
    inv_vn = np.inf

is still a conditional guard in the hot path. Even though it is written as a Python conditional expression rather than a multiline if/else, it is semantically still equivalent to an if vn != 0.0 branch. Numba may lower it to a branch or a conditional select depending on the generated code, but either way it is still a predicate inside the kernel.

We should be able to handle it like the Degeneracy check and let error/inf be caught by the later try-catch

[hongyuchen1030]

See the later comments: Uxarray actually uses the same criteria as the AccuSphGeom, which is also the common best-practice criteria: 1. the input faces are always counterclockwise 2. The GCA will always pick the minor arc that's less than 180 degree. A new issue should be opened to address the convention cleanup. For detailed definition, please refer to the https://egusphere.copernicus.org/preprints/2026/egusphere-2026-636/egusphere-2026-636.pdf, Section 2: Terminology and Problem Definition about node, edge, face

[hongyuchen1030]

If nx_hi and nx_lo are the output of accucross, then it's highly likely that fl(nx_hi + nx_lo) = nx_hi. To fully utilize the EFT kernels we have, you can try using _cdp2(a0, b0, a1, b1) and the the existed def _cdp8( a0, a1, a2, a3, a4, a5, 0, 0, b0, b1, b2, b3, b4, b5, 0, 0 ) to simulate the _cdp6(You can drop in 0 in the last two elements): you can refer to my old C++ implementation below. I've since dropped this EFT support from AccuSphGeom, but you should still be able to view it here: https://github.com/hongyuchen1030/AccuSphGeom/blob/8e2e35c7b1838cf86209f1d4a992fd2d5cabc251/include/spip/predicates/eft/basic.hpp#L282-L313. Hope it helps!

// -----------------------------------------------------------------------------
// Compensated scalar triple product.
// Returns an accurate approximation to
//
//   a · (b x c).
//
// The cross product b x c is first represented componentwise in two-term form,
//
//   (b x c)_k = h_k + l_k,   k = 0,1,2,
//
// using three compensated differences of products. The final triple product is
// then evaluated as the compensated dot product
//
//   a_0 (h_0 + l_0) + a_1 (h_1 + l_1) + a_2 (h_2 + l_2).
//
// This quantity is the determinant det[a, b, c], i.e., the orient3d sign on the
// sphere up to the usual geometric interpretation.
// -----------------------------------------------------------------------------
template <typename T>
inline T compensated_triple_product(const V3_T<T>& a,
                                    const V3_T<T>& b,
                                    const V3_T<T>& c) {
  const TwoTerm<T> x = _cdp2(b[1], c[2], b[2], c[1]);
  const TwoTerm<T> y = _cdp2(b[2], c[0], b[0], c[2]);
  const TwoTerm<T> z = _cdp2(b[0], c[1], b[1], c[0]);
  const TwoTerm<T> dot = _cdp6(
      a[0], x.hi, a[0], x.lo,
      a[1], y.hi, a[1], y.lo,
      a[2], z.hi, a[2], z.lo);
  return dot.hi + dot.lo;
}

[hongyuchen1030]

See the comment above about using my previous C++ implementation

[hongyuchen1030]

Since this still uses if-else branching, we shouldn't call it a "Branch-free mask form" — the comment doesn't match the implementation. We don't need to fix the actual branching in this PR; removing the incorrect claim from the comment is enough for now, since we want to keep moving fast. I'll open a separate issue documenting this mismatch for a future fix.

A true branch-free mask form is what's implemented in AccuSphGeom's on_minor_arc.hpp, [line 43](https://github.com/hongyuchen1030/AccuSphGeom/blob/90b76b2e68220c15dbe1c85c00b92b385da19624/include/accusphgeom/predicates/on_minor_arc.hpp#L43). It relies on the fact that a boolean maps to 1/0, so multiplication, XOR, and OR can select values without any if-else.

For background, see [Avoiding Shader Conditionals](https://theorangeduck.com/page/avoiding-shader-conditionals). Worth noting for our case specifically: this trick pays off most when the condition isn't uniform/coherent across threads — and that's exactly our situation, since every edge and every intersection point has a different shape. There's no shared constant to branch on, so this is inherently a case-by-case, per-element condition rather than something a GPU could cheaply branch on anyway.

[hongyuchen1030]

The so-called "co-planar" or "co-linear" degeneracy case is an extremely sensitive problem that sits inside a genuinely active research area — the general study of how to detect these degeneracies and prove global consistency around them, grounded in mesh topology and math, is still ongoing work in computational geometry. So a claim to "resolve" or "detect" such cases isn't something you can casually assert: it needs to guarantee global consistency, and that guarantee requires substantial proof grounded in mesh topology and math. This is a predicate problem, it requires robustness, and as I mentioned here in another PR, UXarray won't claim any "robustness" due to our domain nature. #1566

The current best-practice, proven solution is "Simulation of Simplicity," which both AccuSphGeom and S2Geometry use. That's why AccuSphGeom's README doesn't recommend the approach at on_minor_arc.hpp, line 43, and instead consistently advocates for Simulation of Simplicity.

We don't need any actual implementation change here — just documenting our strategy around |a x b|^2 and removing the claim that

"|a x b|^2 cleanly separates true degeneracy (~1e-32, the squared rounding-noise floor for unit vectors) from even sub-millimeter real edges (~1e-28 at 1e-12 degree separation)"

will be enough. (All these numbers claims are normally used for describing accuracy for construction problem instead of the predicate problem)