diff --git a/benchmarks/geometry_kernels.py b/benchmarks/geometry_kernels.py new file mode 100644 index 000000000..08aea236e --- /dev/null +++ b/benchmarks/geometry_kernels.py @@ -0,0 +1,207 @@ +"""Micro-benchmarks for the EFT-based spherical geometry kernels. + +These benchmarks target the individual functions introduced in the AccuSphGeom +port (PR #1513) so that the per-kernel overhead of compensated arithmetic can +be measured independently of higher-level UXarray operations. + +All Numba functions are warmed (compiled) during ``setup`` so that benchmark +timings reflect steady-state throughput, not JIT compilation. + +The ``uxarray`` kernels are imported inside each ``setup`` (not at module +import) so that if a symbol is missing in the environment asv built, only the +affected benchmark errors out rather than aborting collection of every +benchmark in this directory. +""" + +import numpy as np + + +def _unit(v): + return v / np.linalg.norm(v) + + +# --------------------------------------------------------------------------- +# Representative inputs — chosen to exercise the near-tangent regime +# --------------------------------------------------------------------------- + +# Two arcs that intersect at a small angle (near-tangent, stress-tests EFT) +_W0 = _unit(np.array([1.0, 0.0, 0.1])) +_W1 = _unit(np.array([0.0, 1.0, 0.1])) +_V0 = _unit(np.array([0.5, -0.1, 0.8])) +_V1 = _unit(np.array([0.5, 0.9, 0.05])) + +# Arc for const-lat test +_X1 = _unit(np.array([1.0, 0.0, 0.3])) +_X2 = _unit(np.array([0.0, 1.0, 0.3])) +_CONST_Z = 0.3 + + +class EFTPrimitives: + """Benchmark the low-level EFT building blocks: two_sum, two_prod, + diff_of_products, and acc_sqrt_re.""" + + def setup(self): + from uxarray.utils.computing import ( + acc_sqrt_re, + diff_of_products, + two_prod, + two_sum, + ) + + self.two_sum = two_sum + self.two_prod = two_prod + self.diff_of_products = diff_of_products + self.acc_sqrt_re = acc_sqrt_re + + # Warm Numba + two_sum(1.0, 1e-16) + two_prod(1.23456789, 9.87654321) + diff_of_products(1.0, 2.0, 3.0, 4.0) + acc_sqrt_re(1.0 - 1e-15) + + def time_two_sum(self): + self.two_sum(1.23456789012345678, 9.87654321098765432e-16) + + def time_two_prod(self): + self.two_prod(1.23456789012345678, 9.87654321098765432) + + def time_diff_of_products(self): + self.diff_of_products(1.23456789, 9.87654321, 1.23456788, 9.87654322) + + def time_acc_sqrt_re(self): + self.acc_sqrt_re(1.0 - 1e-15) + + +class AccucrossKernels: + """Benchmark the compensated cross-product kernels.""" + + def setup(self): + from uxarray.utils.computing import accucross, accucross_pair + + self.accucross = accucross + self.accucross_pair = accucross_pair + + accucross(_W0[0], _W0[1], _W0[2], _W1[0], _W1[1], _W1[2]) + accucross_pair( + 1.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 1.0, + 0.0, + 0.0, + 0.0, + 0.0, + ) + + def time_accucross(self): + self.accucross(_W0[0], _W0[1], _W0[2], _W1[0], _W1[1], _W1[2]) + + def time_accucross_pair(self): + n1x_hi, n1y_hi, n1z_hi, n1x_lo, n1y_lo, n1z_lo = self.accucross( + _W0[0], _W0[1], _W0[2], _W1[0], _W1[1], _W1[2] + ) + n2x_hi, n2y_hi, n2z_hi, n2x_lo, n2y_lo, n2z_lo = self.accucross( + _V0[0], _V0[1], _V0[2], _V1[0], _V1[1], _V1[2] + ) + self.accucross_pair( + n1x_hi, + n1y_hi, + n1z_hi, + n1x_lo, + n1y_lo, + n1z_lo, + n2x_hi, + n2y_hi, + n2z_hi, + n2x_lo, + n2y_lo, + n2z_lo, + ) + + +class OrientPredicates: + """Benchmark the orient3d and on_minor_arc predicates.""" + + def setup(self): + from uxarray.grid.arcs import on_minor_arc, orient3d_on_sphere + + self.orient3d_on_sphere = orient3d_on_sphere + self.on_minor_arc = on_minor_arc + + orient3d_on_sphere(_W0, _W1, _V0) + on_minor_arc(_V0, _W0, _W1) + + def time_orient3d_on_sphere(self): + self.orient3d_on_sphere(_W0, _W1, _V0) + + def time_on_minor_arc(self): + self.on_minor_arc(_V0, _W0, _W1) + + +class GCAGCAIntersection: + """Benchmark all three layers of the GCA-GCA intersection stack.""" + + def setup(self): + from uxarray.grid.intersections import ( + _accux_gca, + _try_gca_gca_intersection, + gca_gca_intersection, + ) + + self._accux_gca = _accux_gca + self._try_gca_gca_intersection = _try_gca_gca_intersection + self.gca_gca_intersection = gca_gca_intersection + + self.gca_a = np.stack([_W0, _W1]) + self.gca_b = np.stack([_V0, _V1]) + _accux_gca(_W0, _W1, _V0, _V1) + _try_gca_gca_intersection(_W0, _W1, _V0, _V1) + gca_gca_intersection(self.gca_a, self.gca_b) + + def time_accux_gca_kernel(self): + """Layer 1: pure numerical kernel.""" + self._accux_gca(_W0, _W1, _V0, _V1) + + def time_try_gca_gca_intersection(self): + """Layer 2: batch/status layer.""" + self._try_gca_gca_intersection(_W0, _W1, _V0, _V1) + + def time_gca_gca_intersection(self): + """Layer 3: dispatcher (full public API).""" + self.gca_gca_intersection(self.gca_a, self.gca_b) + + +class GCAConstLatIntersection: + """Benchmark all three layers of the GCA / constant-latitude intersection stack.""" + + def setup(self): + from uxarray.grid.intersections import ( + _accux_constlat, + _try_gca_const_lat_intersection, + gca_const_lat_intersection, + ) + + self._accux_constlat = _accux_constlat + self._try_gca_const_lat_intersection = _try_gca_const_lat_intersection + self.gca_const_lat_intersection = gca_const_lat_intersection + + self.gca_cart = np.stack([_X1, _X2]) + _accux_constlat(_X1, _X2, _CONST_Z) + _try_gca_const_lat_intersection(self.gca_cart, _CONST_Z) + gca_const_lat_intersection(self.gca_cart, _CONST_Z) + + def time_accux_constlat_kernel(self): + """Layer 1: pure numerical kernel.""" + self._accux_constlat(_X1, _X2, _CONST_Z) + + def time_try_gca_const_lat_intersection(self): + """Layer 2: batch/status layer.""" + self._try_gca_const_lat_intersection(self.gca_cart, _CONST_Z) + + def time_gca_const_lat_intersection(self): + """Layer 3: dispatcher (full public API).""" + self.gca_const_lat_intersection(self.gca_cart, _CONST_Z) diff --git a/benchmarks/geometry_samebody.py b/benchmarks/geometry_samebody.py new file mode 100644 index 000000000..7712d1c7a --- /dev/null +++ b/benchmarks/geometry_samebody.py @@ -0,0 +1,420 @@ +"""Same-body FP64-vs-AccuX diagnostic for the GCA / constant-latitude path. + +Purpose (PR #1513) +------------------ +This is *not* a standalone benchmark. It is a diagnostic that isolates the +engineering overhead of the UXarray AccuX wiring (wrapper / status layer / +dispatcher / masking) from the cost of the compensated-arithmetic (EFT) kernel +itself. + +Method +------ +We build a second implementation of the exact same three-layer stack used by the +real AccuX path in ``uxarray.grid.intersections``: + + L1 kernel pure numerical core, returns two candidate points + L2 try/status finiteness + on-minor-arc masks, branchless point select + L3 dispatcher endpoint snapping, UXarray (2, 3) NaN-filled output + +The only difference is the L1 body: here it is the *direct FP64* formula taken +verbatim from the AccuSphGeom reference + + tests/performance_test/gca_constLat/fp64_GCAconstLat.hh + + nx = a1*b2 - a2*b1; ny = a2*b0 - a0*b2; nz = a0*b1 - a1*b0 + denom = nx^2 + ny^2 + norm_n2 = denom + nz^2 + s = sqrt(denom - norm_n2 * z^2) + pos = ( -(z*nx*nz - s*ny)/denom, -(z*ny*nz + s*nx)/denom, z ) + neg = ( -(z*nx*nz + s*ny)/denom, -(z*ny*nz - s*nx)/denom, z ) + +Because L2/L3 here are byte-for-byte the same logic as the real AccuX L2/L3, the +comparison factors cleanly: + + real AccuX time - same-body FP64 time == cost of EFT math only + same-body FP64 time - direct FP64 time == cost of L2/L3 plumbing only + +Acceptance (interpret ``main()`` output) + 1. same-body FP64 dispatcher output/status == direct FP64 output/status + (exact match: the plumbing does not change results) + 2. same-body FP64 dispatcher output/status == real AccuX output/status + within tolerance (same algorithm, EFT only tightens rounding) + 3. plumbing overhead (L3 - L1) is small and comparable for both bodies + 4. remaining real-AccuX-vs-same-body gap is attributable to EFT math, not + to wrappers/dispatch + +Run directly: ``python benchmarks/geometry_samebody.py`` +This module is import-safe for asv (timing classes at the bottom). +""" + +import math +import time + +import numpy as np +from numba import njit + +from uxarray.grid.arcs import _on_minor_arc_xyz, on_minor_arc +from uxarray.grid.intersections import ( + _accux_constlat_scalar, + _snap_const_lat_endpoint_xy, + gca_const_lat_intersection, +) + +# --------------------------------------------------------------------------- +# L1 (FP64 body) — direct double-precision kernel, verbatim from +# fp64_GCAconstLat.hh. Scalar in / scalar out so Numba keeps it in registers, +# mirroring _accux_constlat_scalar. +# --------------------------------------------------------------------------- + + +@njit(cache=True) +def _fp64_constlat_scalar(a0, a1, a2, b0, b1, b2, const_z): + nx = a1 * b2 - a2 * b1 + ny = a2 * b0 - a0 * b2 + nz = a0 * b1 - a1 * b0 + + denom = nx * nx + ny * ny + norm_n2 = denom + nz * nz + s = math.sqrt(denom - norm_n2 * const_z * const_z) + + inv_denom = 1.0 / denom if denom != 0.0 else np.inf + px = -(const_z * nx * nz - s * ny) * inv_denom + py = -(const_z * ny * nz + s * nx) * inv_denom + nxo = -(const_z * nx * nz + s * ny) * inv_denom + nyo = -(const_z * ny * nz - s * nx) * inv_denom + return px, py, nxo, nyo + + +@njit(cache=True, inline="always") +def _fp64_constlat(x1, x2, const_z): + px, py, nxo, nyo = _fp64_constlat_scalar( + x1[0], x1[1], x1[2], x2[0], x2[1], x2[2], const_z + ) + pos = np.empty(3) + pos[0] = px + pos[1] = py + pos[2] = const_z + neg = np.empty(3) + neg[0] = nxo + neg[1] = nyo + neg[2] = const_z + return pos, neg + + +# --------------------------------------------------------------------------- +# L2 (FP64 body) — identical logic to _try_gca_const_lat_intersection, only the +# L1 call differs. Branchless integer masks; status codes 0/1/2 as in AccuSphGeom. +# --------------------------------------------------------------------------- + + +@njit(cache=True) +def _fp64_try_gca_const_lat_intersection(gca_cart, const_z): + x1 = gca_cart[0] + x2 = gca_cart[1] + pos, neg = _fp64_constlat(x1, x2, const_z) + + pos_fin = int(math.isfinite(pos[0]) and math.isfinite(pos[1])) + neg_fin = int(math.isfinite(neg[0]) and math.isfinite(neg[1])) + pos_on = pos_fin * int(on_minor_arc(pos, x1, x2)) if pos_fin else 0 + neg_on = neg_fin * int(on_minor_arc(neg, x1, x2)) if neg_fin else 0 + + pos_valid = pos_fin * pos_on + neg_valid = neg_fin * neg_on + + pos_mask = pos_valid * (1 - neg_valid) + neg_mask = neg_valid * (1 - pos_valid) + + point = np.empty(3) + point[0] = pos_mask * pos[0] + neg_mask * neg[0] + point[1] = pos_mask * pos[1] + neg_mask * neg[1] + point[2] = pos_mask * pos[2] + neg_mask * neg[2] + + both = pos_valid * neg_valid + none = (1 - pos_valid) * (1 - neg_valid) + status = both + none * 2 + return point, status, pos, neg + + +# --------------------------------------------------------------------------- +# L3 (FP64 body) — identical dispatcher to gca_const_lat_intersection, reusing +# the production _snap_const_lat_endpoint so only the numerical body differs. +# --------------------------------------------------------------------------- + + +@njit(cache=True) +def _fp64_gca_const_lat_intersection(gca_cart, const_z): + # Mirrors the production scalar dispatcher exactly (same allocation profile: + # one (2, 3) array), only the L1 body differs. This keeps the same-body + # comparison honest: any timing gap is the EFT math, not plumbing. + res = np.empty((2, 3)) + res.fill(np.nan) + + a0 = gca_cart[0, 0] + a1 = gca_cart[0, 1] + a2 = gca_cart[0, 2] + b0 = gca_cart[1, 0] + b1 = gca_cart[1, 1] + b2 = gca_cart[1, 2] + + px, py, nx, ny = _fp64_constlat_scalar(a0, a1, a2, b0, b1, b2, const_z) + + pos_fin = math.isfinite(px) and math.isfinite(py) + neg_fin = math.isfinite(nx) and math.isfinite(ny) + pos_valid = pos_fin and _on_minor_arc_xyz(px, py, const_z, a0, a1, a2, b0, b1, b2) + neg_valid = neg_fin and _on_minor_arc_xyz(nx, ny, const_z, a0, a1, a2, b0, b1, b2) + + if pos_valid and not neg_valid: + sx, sy = _snap_const_lat_endpoint_xy(px, py, a0, a1, a2, b0, b1, b2, const_z) + res[0, 0] = sx + res[0, 1] = sy + res[0, 2] = const_z + elif neg_valid and not pos_valid: + sx, sy = _snap_const_lat_endpoint_xy(nx, ny, a0, a1, a2, b0, b1, b2, const_z) + res[0, 0] = sx + res[0, 1] = sy + res[0, 2] = const_z + elif pos_valid and neg_valid: + psx, psy = _snap_const_lat_endpoint_xy(px, py, a0, a1, a2, b0, b1, b2, const_z) + nsx, nsy = _snap_const_lat_endpoint_xy(nx, ny, a0, a1, a2, b0, b1, b2, const_z) + dx = psx - nsx + dy = psy - nsy + if dx * dx + dy * dy < 1e-14: + res[0, 0] = psx + res[0, 1] = psy + res[0, 2] = const_z + else: + res[0, 0] = psx + res[0, 1] = psy + res[0, 2] = const_z + res[1, 0] = nsx + res[1, 1] = nsy + res[1, 2] = const_z + return res + + +# --------------------------------------------------------------------------- +# Test inputs +# --------------------------------------------------------------------------- + + +def _unit(v): + return v / np.linalg.norm(v) + + +def _make_cases(n, seed, frac_cross=0.7): + + rng = np.random.default_rng(seed) + cases = [] + while len(cases) < n: + mid = _unit(rng.standard_normal(3)) + tan = _unit(np.cross(mid, _unit(rng.standard_normal(3)))) + half = 0.5 * math.radians(10.0 ** rng.uniform(math.log10(0.2), math.log10(8.0))) + ca, sa = math.cos(half), math.sin(half) + a = _unit(mid * ca - tan * sa) + b = _unit(mid * ca + tan * sa) + zlo, zhi = (a[2], b[2]) if a[2] <= b[2] else (b[2], a[2]) + if rng.random() < frac_cross and zhi - zlo > 1e-9: + const_z = zlo + (zhi - zlo) * rng.uniform(0.15, 0.85) + else: + # a latitude the arc does not span -> empty-return path + room_above = 0.999 - zhi + room_below = zlo + 0.999 + if room_above >= room_below: + const_z = zhi + min(room_above, rng.uniform(0.02, 0.3)) + else: + const_z = zlo - min(room_below, rng.uniform(0.02, 0.3)) + cases.append((np.stack([a, b]), float(const_z))) + return cases + + +# --------------------------------------------------------------------------- +# Batched drivers — the timing loop lives *inside* njit, mirroring the +# AccuSphGeom C++ benchmark (2M points looped in-kernel). Timing a Python-level +# per-call loop would drown the signal in interpreter overhead; batching in +# Numba measures true kernel throughput, which is also how UXarray actually +# calls these paths (once per edge over a whole grid). +# --------------------------------------------------------------------------- + + +@njit(cache=True) +def _batch_accux_kernel(A, B, Z): + """Real AccuX L1 (EFT) kernel over a batch; accumulate to defeat DCE.""" + acc = 0.0 + for i in range(A.shape[0]): + px, py, nxo, nyo = _accux_constlat_scalar( + A[i, 0], A[i, 1], A[i, 2], B[i, 0], B[i, 1], B[i, 2], Z[i] + ) + acc += px + py + nxo + nyo + return acc + + +@njit(cache=True) +def _batch_fp64_kernel(A, B, Z): + """Same-body FP64 L1 kernel over a batch; accumulate to defeat DCE.""" + acc = 0.0 + for i in range(A.shape[0]): + px, py, nxo, nyo = _fp64_constlat_scalar( + A[i, 0], A[i, 1], A[i, 2], B[i, 0], B[i, 1], B[i, 2], Z[i] + ) + acc += px + py + nxo + nyo + return acc + + +@njit(cache=True) +def _batch_accux_dispatch(gcas, Z): + """Real AccuX full L1+L2+L3 dispatcher over a batch.""" + acc = 0.0 + for i in range(gcas.shape[0]): + res = gca_const_lat_intersection(gcas[i], Z[i]) + v = res[0, 0] + if v == v: # not NaN + acc += v + return acc + + +@njit(cache=True) +def _batch_fp64_dispatch(gcas, Z): + """Same-body FP64 full L1+L2+L3 dispatcher over a batch.""" + acc = 0.0 + for i in range(gcas.shape[0]): + res = _fp64_gca_const_lat_intersection(gcas[i], Z[i]) + v = res[0, 0] + if v == v: # not NaN + acc += v + return acc + + +def _time_batch(fn, args, repeat=7): + """Best-of-`repeat` wall-time for one batched call (compile excluded).""" + fn(*args) # warm / compile + best = math.inf + for _ in range(repeat): + t0 = time.perf_counter() + fn(*args) + best = min(best, time.perf_counter() - t0) + return best + + +def _pack(cases): + """Turn the case list into contiguous arrays for the batched kernels.""" + A = np.array([c[0][0] for c in cases]) + B = np.array([c[0][1] for c in cases]) + Z = np.array([c[1] for c in cases]) + gcas = np.array([c[0] for c in cases]) + return A, B, Z, gcas + + +# --------------------------------------------------------------------------- +# Diagnostic driver +# --------------------------------------------------------------------------- + + +def main(): + base_cases = _make_cases(200, seed=20251104) + + # ---- correctness (on the 200 baseline cases) ---- + max_out_diff = 0.0 + status_mismatch = 0 + n_with_result = 0 + for gca, z in base_cases: + fp64_res = _fp64_gca_const_lat_intersection(gca, z) + accux_res = gca_const_lat_intersection(gca, z) + fp64_rows = int(np.isfinite(fp64_res[0, 0])) + int(np.isfinite(fp64_res[1, 0])) + accux_rows = int(np.isfinite(accux_res[0, 0])) + int( + np.isfinite(accux_res[1, 0]) + ) + if fp64_rows != accux_rows: + status_mismatch += 1 + if fp64_rows > 0 and accux_rows > 0: + n_with_result += 1 + d = np.nanmax(np.abs(fp64_res - accux_res)) + if np.isfinite(d): + max_out_diff = max(max_out_diff, d) + + # ---- timing: a large batch of DISTINCT cases (no replication) ---- + A, B, Z, gcas = _pack(_make_cases(100_000, seed=20251105)) + n = A.shape[0] + + t_direct = _time_batch(_batch_fp64_kernel, (A, B, Z)) + t_accux_k = _time_batch(_batch_accux_kernel, (A, B, Z)) + t_fp64_d = _time_batch(_batch_fp64_dispatch, (gcas, Z)) + t_accux_d = _time_batch(_batch_accux_dispatch, (gcas, Z)) + + ns_per = lambda t: t / n * 1e9 # noqa: E731 + + print("=" * 70) + print("Same-body FP64-vs-AccuX diagnostic — GCA/ConstLat (PR #1513)") + print("=" * 70) + print(f"baseline cases: {len(base_cases)} with-result: {n_with_result}") + print(f"timing batch : {n} distinct points, best of 7") + print() + print("CORRECTNESS (same-body FP64 dispatcher vs real AccuX dispatcher)") + print(f" status mismatches : {status_mismatch}") + print(f" max output diff : {max_out_diff:.3e}") + print() + print("TIMING (ns per point, in-kernel batch)") + print(f" L1 FP64 kernel : {ns_per(t_direct):8.2f} ns") + print(f" L1 AccuX kernel : {ns_per(t_accux_k):8.2f} ns") + print(f" L1+L2+L3 FP64 dispatch : {ns_per(t_fp64_d):8.2f} ns") + print(f" L1+L2+L3 AccuX dispatch : {ns_per(t_accux_d):8.2f} ns") + print() + print("DECOMPOSITION") + plumb_fp64 = t_fp64_d - t_direct + plumb_accux = t_accux_d - t_accux_k + eft_kernel = t_accux_k - t_direct + eft_dispatch = t_accux_d - t_fp64_d + print( + f" plumbing L2/L3 over FP64 body : {ns_per(plumb_fp64):8.2f} ns/pt" + f" ({100 * plumb_fp64 / t_fp64_d:.1f}% of FP64 dispatch)" + ) + print( + f" plumbing L2/L3 over AccuX body : {ns_per(plumb_accux):8.2f} ns/pt" + f" ({100 * plumb_accux / t_accux_d:.1f}% of AccuX dispatch)" + ) + print( + f" EFT math cost (kernel level) : {ns_per(eft_kernel):8.2f} ns/pt" + f" ({t_accux_k / t_direct:.2f}x FP64 kernel)" + ) + print( + f" EFT math cost (dispatch level) : {ns_per(eft_dispatch):8.2f} ns/pt" + f" ({t_accux_d / t_fp64_d:.2f}x FP64 dispatch)" + ) + print() + print("INTERPRETATION") + print(" - plumbing overhead should be ~equal for both bodies (body-independent)") + print(" - AccuX/FP64 ratio at dispatch ~= ratio at kernel => plumbing adds no") + print(" EFT-dependent overhead; measured cost is the EFT math, wired correctly") + print("=" * 70) + + +# --------------------------------------------------------------------------- +# asv timing classes (batched; Numba warmed in setup) +# --------------------------------------------------------------------------- + + +class SameBodyConstLat: + """asv: same-body FP64 vs real AccuX at kernel (L1) and dispatch (L3) levels.""" + + def setup(self): + cases = _make_cases(20_000, seed=20251104) + self.A, self.B, self.Z, self.gcas = _pack(cases) + _batch_fp64_kernel(self.A, self.B, self.Z) + _batch_accux_kernel(self.A, self.B, self.Z) + _batch_fp64_dispatch(self.gcas, self.Z) + _batch_accux_dispatch(self.gcas, self.Z) + + def time_fp64_kernel(self): + _batch_fp64_kernel(self.A, self.B, self.Z) + + def time_accux_kernel(self): + _batch_accux_kernel(self.A, self.B, self.Z) + + def time_fp64_dispatch(self): + _batch_fp64_dispatch(self.gcas, self.Z) + + def time_accux_dispatch(self): + _batch_accux_dispatch(self.gcas, self.Z) + + +if __name__ == "__main__": + main() diff --git a/benchmarks/geometry_samebody_gcagca.py b/benchmarks/geometry_samebody_gcagca.py new file mode 100644 index 000000000..348f55859 --- /dev/null +++ b/benchmarks/geometry_samebody_gcagca.py @@ -0,0 +1,338 @@ +"""Same-body FP64-vs-AccuX diagnostic for the GCA x GCA path. + +Companion to ``geometry_samebody.py``, which only covers the GCA/constant-latitude +stack. T + +gca_gca is the kernel invoked per edge from the njit point-in-polygon crossing +loop (``uxarray/grid/geometry.py`` ``_check_intersection``), so its allocation +profile is on the hot path of ``get_faces_containing_point`` and pole detection. + +Factoring (identical to ``geometry_samebody.py``):: + + real AccuX time - same-body FP64 time == EFT math only + same-body FP64 - direct FP64 kernel == L2/L3 plumbing only + +""" + +import math +import time + +import numpy as np +from numba import njit + +from uxarray.grid.arcs import on_minor_arc +from uxarray.grid.intersections import _accux_gca, gca_gca_intersection + + +@njit(cache=True, inline="always") +def _fp64_gca(w0, w1, v0, v1): + """ + L1 (FP64 body) -- plain double-precision cross-product triple, the direct + analogue of _accux_gca (intersections.py) with accucross/accucross_pair + replaced by naive FP64 cross products. Same allocation shape (two np.empty(3)) + so the twin's allocation profile matches the real kernel exactly. + """ + + n1x = w0[1] * w1[2] - w0[2] * w1[1] + n1y = w0[2] * w1[0] - w0[0] * w1[2] + n1z = w0[0] * w1[1] - w0[1] * w1[0] + n2x = v0[1] * v1[2] - v0[2] * v1[1] + n2y = v0[2] * v1[0] - v0[0] * v1[2] + n2z = v0[0] * v1[1] - v0[1] * v1[0] + + vx = n1y * n2z - n1z * n2y + vy = n1z * n2x - n1x * n2z + vz = n1x * n2y - n1y * n2x + + vn = math.sqrt(vx * vx + vy * vy + vz * vz) + inv = 1.0 / vn if vn != 0.0 else np.inf + pos = np.empty(3) + pos[0] = vx * inv + pos[1] = vy * inv + pos[2] = vz * inv + neg = np.empty(3) + neg[0] = -pos[0] + neg[1] = -pos[1] + neg[2] = -pos[2] + return pos, neg + + +@njit(cache=True) +def _fp64_try_gca_gca_intersection(w0, w1, v0, v1): + """ + L2 (FP64 body) -- byte-for-byte identical logic to _try_gca_gca_intersection + (intersections.py), only the L1 call differs. + """ + pos, neg = _fp64_gca(w0, w1, v0, v1) + + pos_fin = ( + 1 + if math.isfinite(pos[0]) and math.isfinite(pos[1]) and math.isfinite(pos[2]) + else 0 + ) + neg_fin = ( + 1 + if math.isfinite(neg[0]) and math.isfinite(neg[1]) and math.isfinite(neg[2]) + else 0 + ) + pos_on_a = 1 if (pos_fin and on_minor_arc(pos, w0, w1)) else 0 + pos_on_b = 1 if (pos_fin and on_minor_arc(pos, v0, v1)) else 0 + neg_on_a = 1 if (neg_fin and on_minor_arc(neg, w0, w1)) else 0 + neg_on_b = 1 if (neg_fin and on_minor_arc(neg, v0, v1)) else 0 + + pos_valid = pos_fin * pos_on_a * pos_on_b + neg_valid = neg_fin * neg_on_a * neg_on_b + + pos_mask = pos_valid * (1 - neg_valid) + neg_mask = neg_valid * (1 - pos_valid) + + point = np.empty(3) + point[0] = pos_mask * pos[0] + neg_mask * neg[0] + point[1] = pos_mask * pos[1] + neg_mask * neg[1] + point[2] = pos_mask * pos[2] + neg_mask * neg[2] + + both = pos_valid * neg_valid + none = (1 - pos_valid) * (1 - neg_valid) + status = both + none * 2 + return point, status, pos, neg + + +@njit(cache=True) +def _fp64_gca_gca_intersection(gca_a_xyz, gca_b_xyz): + """ + L3 (FP64 body) -- identical dispatcher to gca_gca_intersection + (intersections.py), same np.empty((2, 3)) + res[:count] slice profile. + """ + if gca_a_xyz.shape[1] != 3 or gca_b_xyz.shape[1] != 3: + raise ValueError("The two GCAs must be in the cartesian [x, y, z] format") + + w0 = gca_a_xyz[0] + w1 = gca_a_xyz[1] + v0 = gca_b_xyz[0] + v1 = gca_b_xyz[1] + + point, status, pos, neg = _fp64_try_gca_gca_intersection(w0, w1, v0, v1) + + res = np.empty((2, 3)) + count = 0 + if status == 0: + res[0, 0] = point[0] + res[0, 1] = point[1] + res[0, 2] = point[2] + count = 1 + elif status == 1: + res[0, 0] = pos[0] + res[0, 1] = pos[1] + res[0, 2] = pos[2] + res[1, 0] = neg[0] + res[1, 1] = neg[1] + res[1, 2] = neg[2] + count = 2 + else: + if on_minor_arc(v0, w0, w1): + res[count, 0] = v0[0] + res[count, 1] = v0[1] + res[count, 2] = v0[2] + count += 1 + if on_minor_arc(v1, w0, w1): + res[count, 0] = v1[0] + res[count, 1] = v1[1] + res[count, 2] = v1[2] + count += 1 + return res[:count] + + +# --------------------------------------------------------------------------- +# Case generation -- a LARGE set of DISTINCT, SHORT arcs (matching real grid +# edges, ~0.2-8 deg / median ~1.5 deg) with a controlled intersect / no-intersect +# mix. Short arcs put a x b in the near-parallel cancellation regime the EFT +# targets; the mix exercises both the status 0/1 (allocating) and status 2 +# (empty-return) dispatcher branches. +# --------------------------------------------------------------------------- + + +def _unit(v): + return v / np.linalg.norm(v) + + +def _short_arc(rng, mid=None): + """A short great-circle arc (~0.2-8 deg, median ~1.5 deg) around ``mid`` (a + random point if not given), matching real unstructured-grid edge lengths.""" + if mid is None: + mid = _unit(rng.standard_normal(3)) + tan = _unit(np.cross(mid, _unit(rng.standard_normal(3)))) + half = 0.5 * math.radians(10.0 ** rng.uniform(math.log10(0.2), math.log10(8.0))) + ca, sa = math.cos(half), math.sin(half) + return _unit(mid * ca - tan * sa), _unit(mid * ca + tan * sa) + + +def _make_gca_cases(n, seed, frac_intersect=0.6): + rng = np.random.default_rng(seed) + cases = [] + while len(cases) < n: + if rng.random() < frac_intersect: + # two short arcs sharing a midpoint p -> p is the midpoint of both + # minor arcs, so they intersect there (status 0/1). + p = _unit(rng.standard_normal(3)) + a0, a1 = _short_arc(rng, mid=p) + b0, b1 = _short_arc(rng, mid=p) + else: + # two independent short arcs -> their great circles cross off at + # least one minor arc, exercising the status 2 / empty-return branch. + a0, a1 = _short_arc(rng) + b0, b1 = _short_arc(rng) + cases.append((np.stack([a0, a1]), np.stack([b0, b1]))) + return cases + + +def _pack_gca(cases): + wa = np.ascontiguousarray([c[0][0] for c in cases]) + wb = np.ascontiguousarray([c[0][1] for c in cases]) + va = np.ascontiguousarray([c[1][0] for c in cases]) + vb = np.ascontiguousarray([c[1][1] for c in cases]) + ga = np.ascontiguousarray([c[0] for c in cases]) # (n, 2, 3) + gb = np.ascontiguousarray([c[1] for c in cases]) + return wa, wb, va, vb, ga, gb + + +# --------------------------------------------------------------------------- +# Batched in-kernel drivers -- the timing loop lives inside njit (mirrors +# geometry_samebody.py). Accumulate a scalar to defeat dead-code elimination. +# --------------------------------------------------------------------------- + + +@njit(cache=True) +def _batch_accux_gca_kernel(wa, wb, va, vb): + acc = 0.0 + for i in range(wa.shape[0]): + pos, neg = _accux_gca(wa[i], wb[i], va[i], vb[i]) + acc += pos[0] + pos[1] + pos[2] + return acc + + +@njit(cache=True) +def _batch_fp64_gca_kernel(wa, wb, va, vb): + acc = 0.0 + for i in range(wa.shape[0]): + pos, neg = _fp64_gca(wa[i], wb[i], va[i], vb[i]) + acc += pos[0] + pos[1] + pos[2] + return acc + + +@njit(cache=True) +def _batch_accux_gca_dispatch(ga, gb): + acc = 0.0 + for i in range(ga.shape[0]): + res = gca_gca_intersection(ga[i], gb[i]) + if res.shape[0] > 0: + acc += res[0, 0] + return acc + + +@njit(cache=True) +def _batch_fp64_gca_dispatch(ga, gb): + acc = 0.0 + for i in range(ga.shape[0]): + res = _fp64_gca_gca_intersection(ga[i], gb[i]) + if res.shape[0] > 0: + acc += res[0, 0] + return acc + + +def _time_batch(fn, args, repeat=7): + """Best-of-`repeat` wall-time for one batched call (compile excluded).""" + fn(*args) # warm / compile + best = math.inf + for _ in range(repeat): + t0 = time.perf_counter() + fn(*args) + best = min(best, time.perf_counter() - t0) + return best + + +def main(n_cases=100_000, seed=20251104): + """ + Standalone diagnostic driver (prints ns/edge-pair). Mirrors geometry_samebody.main. + """ + cases = _make_gca_cases(n_cases, seed=seed) + wa, wb, va, vb, ga, gb = _pack_gca(cases) + n = wa.shape[0] + + # correctness: FP64 twin vs real AccuX dispatcher (row count + max diff) + row_mismatch = 0 + max_out_diff = 0.0 + n_check = min(n, 5000) + for i in range(n_check): + r_fp = _fp64_gca_gca_intersection(ga[i], gb[i]) + r_ax = gca_gca_intersection(ga[i], gb[i]) + if r_fp.shape[0] != r_ax.shape[0]: + row_mismatch += 1 + elif r_fp.shape[0] > 0: + max_out_diff = max(max_out_diff, float(np.max(np.abs(r_fp - r_ax)))) + + t_fp64_k = _time_batch(_batch_fp64_gca_kernel, (wa, wb, va, vb)) + t_accux_k = _time_batch(_batch_accux_gca_kernel, (wa, wb, va, vb)) + t_fp64_d = _time_batch(_batch_fp64_gca_dispatch, (ga, gb)) + t_accux_d = _time_batch(_batch_accux_gca_dispatch, (ga, gb)) + + def ns(t): + return t / n * 1e9 + + print("=" * 70) + print("Same-body FP64-vs-AccuX diagnostic -- GCA x GCA") + print("=" * 70) + print(f"distinct cases : {n} (best of 7, in-kernel batch)") + print( + f"correctness : row mismatches {row_mismatch}/{n_check}" + f" max |diff| {max_out_diff:.3e}" + ) + print() + print("TIMING (ns per edge-pair)") + print(f" L1 FP64 kernel : {ns(t_fp64_k):8.2f} ns") + print(f" L1 AccuX kernel : {ns(t_accux_k):8.2f} ns") + print(f" L1+L2+L3 FP64 dispatch : {ns(t_fp64_d):8.2f} ns") + print(f" L1+L2+L3 AccuX dispatch : {ns(t_accux_d):8.2f} ns") + print() + print("DECOMPOSITION") + print(f" plumbing (FP64 body) : {ns(t_fp64_d - t_fp64_k):8.2f} ns/edge") + print(f" plumbing (AccuX body) : {ns(t_accux_d - t_accux_k):8.2f} ns/edge") + print( + f" EFT math (kernel) : {ns(t_accux_k - t_fp64_k):8.2f} ns/edge" + f" ({t_accux_k / t_fp64_k:.2f}x)" + ) + print( + f" EFT math (dispatch) : {ns(t_accux_d - t_fp64_d):8.2f} ns/edge" + f" ({t_accux_d / t_fp64_d:.2f}x)" + ) + print("=" * 70) + + +class SameBodyGcaGca: + """ + asv timing class (Numba warmed in setup, distinct cases) + same-body FP64 vs real AccuX gca_gca at kernel (L1) and dispatch (L3). + """ + + def setup(self): + cases = _make_gca_cases(100_000, seed=20251104) + self.wa, self.wb, self.va, self.vb, self.ga, self.gb = _pack_gca(cases) + _batch_fp64_gca_kernel(self.wa, self.wb, self.va, self.vb) + _batch_accux_gca_kernel(self.wa, self.wb, self.va, self.vb) + _batch_fp64_gca_dispatch(self.ga, self.gb) + _batch_accux_gca_dispatch(self.ga, self.gb) + + def time_fp64_kernel(self): + _batch_fp64_gca_kernel(self.wa, self.wb, self.va, self.vb) + + def time_accux_kernel(self): + _batch_accux_gca_kernel(self.wa, self.wb, self.va, self.vb) + + def time_fp64_dispatch(self): + _batch_fp64_gca_dispatch(self.ga, self.gb) + + def time_accux_dispatch(self): + _batch_accux_gca_dispatch(self.ga, self.gb) + + +if __name__ == "__main__": + main() diff --git a/docs/api.rst b/docs/api.rst index dd89363ba..2addaabd8 100644 --- a/docs/api.rst +++ b/docs/api.rst @@ -594,6 +594,7 @@ Intersections grid.intersections.gca_gca_intersection grid.intersections.gca_const_lat_intersection + grid.intersections.get_number_of_intersections Arcs @@ -605,13 +606,24 @@ Arcs grid.arcs.in_between grid.arcs.point_within_gca grid.arcs.extreme_gca_latitude + grid.arcs.orient3d_on_sphere + grid.arcs.on_minor_arc -Accurate Computing ------------------- +Compensated Arithmetic +---------------------- + +Numba-compiled primitives used throughout the geometry stack to avoid +catastrophic cancellation in cross-product and dot-product operations. +``two_sum`` and ``two_prod`` are true error-free transformations (EFT); +the higher-level functions are compensated algorithms built on top of them. .. autosummary:: :toctree: generated/ - utils.computing.cross_fma - utils.computing.dot_fma + utils.computing.two_sum + utils.computing.two_prod + utils.computing.diff_of_products + utils.computing.accucross + utils.computing.accucross_pair + utils.computing.acc_sqrt_re diff --git a/test/grid/geometry/data/accusphgeom/gca_constlat_cases_with_baseline.csv b/test/grid/geometry/data/accusphgeom/gca_constlat_cases_with_baseline.csv new file mode 100644 index 000000000..ce456ea88 --- /dev/null +++ b/test/grid/geometry/data/accusphgeom/gca_constlat_cases_with_baseline.csv @@ -0,0 +1,201 @@ +case_id,a0_x_sig,a0_x_exp,a0_y_sig,a0_y_exp,a0_z_sig,a0_z_exp,a1_x_sig,a1_x_exp,a1_y_sig,a1_y_exp,a1_z_sig,a1_z_exp,z0_sig,z0_exp,baseline_x_sig,baseline_x_exp,baseline_y_sig,baseline_y_exp +1,6583395283802737,-53,-6147217518692320,-53,8283167787083469,-62,8836420437308160,-54,7848811450056676,-53,4622067348226882,-59,5095978859057569,-60,8174045244549224,-53,-7566515111374263,-54 +2,6733999148524275,-53,-5981653632221706,-53,6667935649031901,-60,-7646311833050094,-53,-4760285605446208,-53,7279949724148809,-60,7074890416577078,-60,6469291231565926,-53,-6266965319637452,-53 +3,4517976994403672,-53,-7792029515625047,-53,5428045630617133,-60,7701242507831454,-53,-4671241778467645,-53,6472814096971289,-65,4929034554476170,-61,6567426310632902,-53,-6164266327058700,-53 +4,8098445525814441,-53,-7879732055419076,-54,4775850045004277,-58,-5023414828919286,-55,-8919200171512326,-53,4707763702530374,-61,4778380633177198,-59,8989209071635874,-55,-8722023646352851,-53 +5,-8952740571655975,-53,7851429899418999,-56,7804931746897990,-59,-8934987439969263,-53,-4551733644395887,-55,6989381849968096,-61,6238802562723750,-59,-8996686952624402,-53,6783692241626225,-57 +6,-5786898905451215,-53,-6901031707988572,-53,8393725786152174,-59,-6360528041351365,-53,-6377198165305838,-53,8742512640693859,-60,6853448270864531,-59,-6012746975386437,-53,-6705598034454887,-53 +7,8851601573671468,-53,6657490412866596,-55,5952974800281886,-59,-8654445919879161,-54,7899410880627245,-53,8183604304089821,-60,4940111606502338,-59,-6775537775520031,-54,8345459975241962,-53 +8,-5975497717992614,-53,6739301458091770,-53,8942880865488835,-60,4858891207079837,-55,-8924913801363269,-53,6172345958047283,-66,5511157628994992,-60,-7353874784237703,-55,-8817472282692888,-53 +9,8163238667756318,-53,-7613376792631668,-54,8813559272091280,-62,8940774306648233,-53,8674866717785775,-56,8186048652772834,-59,6475624385413960,-59,9003466348271009,-53,-7639520501100972,-58 +10,-7748806490431427,-53,4590312954723759,-53,7749912757519997,-59,-8038695191573049,-53,-8126219653059674,-54,6406532728941906,-62,6745474689311411,-59,-8536376238047904,-53,5744148462508969,-54 +11,5729649883672734,-53,6949874106834893,-53,5931605408852635,-66,-5879933382854807,-53,-6823125400944634,-53,8020569392924582,-61,6369009887713425,-61,5614550028759927,-53,7043141871868019,-53 +12,-4966249569289305,-54,-8656815196109413,-53,4882530900349333,-58,8457139260126380,-56,8944928069611956,-53,4790513892235082,-61,4835001749374130,-58,-8568594537403865,-60,9005683084921659,-53 +13,-5747875476700889,-55,8891320418995319,-53,6132018445536323,-59,-7488417550822455,-53,5004779421688677,-53,4713433537297495,-59,4830169762319661,-59,-7336041660477916,-53,5225555974031693,-53 +14,-8010139183394512,-53,8232839025084698,-54,4789222462107549,-58,-8776176964407132,-53,-8107127288347289,-55,5791202524767913,-61,5685818798894064,-59,-8954026416595253,-53,7785757814593620,-56 +15,5627560175581670,-55,8896633830056122,-53,7336316766455182,-62,6874851396255258,-56,8965911779758320,-53,7648905386720497,-60,5969332927342826,-61,5196090421123943,-55,8913003700780547,-53 +16,8492840481767979,-53,5992988266952219,-54,4780962127691762,-58,8126333802577577,-61,-9005994246073443,-53,4603820888419848,-58,4776780970989056,-58,8558959067888302,-58,-9001989597780129,-53 +17,8991398168254582,-54,7804111921337075,-53,7617150522377666,-59,7995995373935763,-54,8070575701889596,-53,6851513020525783,-59,7577539361452820,-59,8939810370258979,-54,7818923692869813,-53 +18,-8076734484706385,-53,7973521805309732,-54,5333136895001468,-60,5243036257483146,-53,-7323341798947998,-53,6028960712903585,-59,8227308920738787,-60,-8376027944943104,-53,6623492400576725,-54 +19,-8915890687492357,-53,5115598768751578,-55,7878982229630873,-61,-7600169027758730,-53,4833842240939289,-53,8249215453370079,-61,8009872551763243,-61,-8832013152887252,-53,7070166173498405,-55 +20,-4736724448686548,-53,7660184092440486,-53,7748912621876295,-59,-7569341567356542,-53,4879796445042261,-53,4777901798692398,-58,4578109227680362,-58,-6778151038123297,-53,5930079178512525,-53 +21,8863330721783680,-53,-6413206970050163,-55,5311326637387805,-61,-5519645992032906,-53,7117165443998166,-53,6106076291744439,-59,5166268472706259,-59,8825953280123746,-53,7184061155473002,-55 +22,8086080013676985,-53,-7933558700937538,-54,6273962701839463,-59,-8099556827757963,-56,8949256060366602,-53,7944848100598596,-59,7556901489559956,-59,8434750995123206,-53,-6315274177663751,-54 +23,-5305073271286867,-54,8606644894383164,-53,8903012640731218,-59,-5772866191323350,-53,6913079913396932,-53,7291637339989954,-59,8289905198187229,-59,-8474364436486911,-54,7947272947560255,-53 +24,7985376412124676,-53,-8331050571851181,-54,6952492827913120,-59,4957346926602060,-53,-7520124867727764,-53,5826232817914260,-60,4636700684183548,-59,6336531532733607,-53,-6400996623310637,-53 +25,7895616628792092,-53,8668701964128537,-54,6109275240341727,-60,-8911516381957835,-53,-5232279220169933,-55,7537421959376491,-60,6895416422135522,-60,7804128684760083,-53,8993844975741971,-54 +26,5068493751719865,-53,-7445414901764326,-53,4876801666806982,-59,6201528532826623,-53,6532230698023442,-53,6498462778455389,-61,6255417990619455,-60,7851755596034654,-53,8826592582319340,-54 +27,6934282139647278,-54,8313143667715727,-53,7482498867876989,-62,8874416502221319,-53,-6142764199358357,-55,8101287544786084,-59,5644293506596004,-59,7790881170547896,-53,4519295427051515,-53 +28,-8843477829400007,-56,8937890042761454,-53,4721931401620408,-58,8354295420244706,-53,-6728268311122358,-54,8583597654645009,-59,4671941379525545,-58,8510641128464153,-53,-5891455009421546,-54 +29,-4828350863063825,-54,8676751930274237,-53,7934797554500779,-59,8712256531076383,-53,-4568840823486536,-54,5596957796491731,-59,5788459360031102,-59,8746459058418786,-53,-8598522770171608,-55 +30,-7195038693246333,-54,8256978556863145,-53,6335055050640220,-59,-9004951321946430,-53,8723872523781415,-59,4736599890711284,-58,6782013859904899,-59,-8187917603035303,-54,8022338208041894,-53 +31,-7595682127031130,-53,-4838892710345401,-53,4567015232295280,-58,7796173741221025,-53,4510611583121942,-53,7782313166486521,-60,6200311855430475,-59,7829832337105621,-53,8902579014573131,-54 +32,8216685945168876,-54,-8015633557163236,-53,7127251037759728,-61,8537784088563638,-54,-7931336696213005,-53,5051449682522472,-62,6398542383320922,-61,8267773547393836,-54,-8002499224712159,-53 +33,8846913732122840,-53,6766401758383690,-55,7994150688684681,-62,6955584853255551,-54,8308700172903407,-53,5112024304895465,-62,6534991628996690,-62,5497316541503676,-53,7135053352723746,-53 +34,-8992481401767967,-53,-8087583037170657,-57,6209236865995871,-59,-9005900983235443,-53,-7697958267222859,-59,6044010545640904,-59,6075445977885175,-59,-9004660049611607,-53,-6132348971696878,-58 +35,-7780189304899989,-54,8123806365678862,-53,6123392331703393,-61,-8143885117639846,-58,9003343425177273,-53,8753890202373835,-60,7449301423208519,-60,-4739047249730611,-55,8928750802543478,-53 +36,6342881777107600,-54,-8430394453418458,-53,7274349768729014,-63,-8993760126535746,-53,7773694181542407,-57,4899766872997044,-59,5426684190061945,-60,-4845588843583573,-56,-8986710755091261,-53 +37,8168540598779432,-56,-8948222611085310,-53,8188880571945480,-59,-5969947608587203,-53,-6744445156887916,-53,5465901348884072,-60,7472770220569310,-59,-7364148965135438,-57,-8994674335912508,-53 +38,-8897855620297720,-54,-7831705023191292,-53,8398143935525417,-61,8034509398478095,-54,8060369282146641,-53,4719336248957064,-58,8700142255342418,-60,-4531423187883382,-53,-7784036382550178,-53 +39,7591024477343437,-53,-4848082944059496,-53,5834282262444449,-60,-8905987354123811,-53,-5380221923986788,-55,7945713626062216,-60,6259494076896527,-60,7478703979720120,-53,-5019585019065492,-53 +40,5987430154660489,-54,8495119492787630,-53,8147939621631735,-62,6730780283665034,-55,-8848604255117354,-53,4980165598055359,-61,8189385753913654,-62,6008247605624297,-54,8491444091280850,-53 +41,8414635644223287,-53,-6425030829989610,-54,7342000104867337,-60,-6330876054247496,-53,-6406681042487707,-53,8180812551862499,-60,7753832247845779,-60,8299858821475927,-53,-6996659891041631,-54 +42,-4666817904046281,-54,-8699572930902831,-53,6101046020132618,-60,8244342276723334,-54,-8007736426867215,-53,7436358951093286,-59,4705024542936948,-59,-5062148090518845,-57,-9001340723129830,-53 +43,-5342566560465282,-53,7251458669390089,-53,6973507350239644,-60,-4660725756527927,-53,7706643539651801,-53,7817224658896418,-59,7383906056264580,-59,-4731179018768310,-53,7663698353605605,-53 +44,8508864583862531,-53,-5908544558135436,-54,8632968435430613,-61,-7912757399833771,-53,8603200741862095,-54,7611156722436436,-59,7248296240736217,-59,8729286816728359,-53,-8869149654809751,-55 +45,-7042738495110928,-53,5614577553524815,-53,4954029767256887,-59,7785291444733240,-53,-4529734246442397,-53,5011272922082843,-61,7835364514562861,-60,8053784874790610,-53,-8065343260567053,-54 +46,5266239946738295,-55,-8909539141896448,-53,8201118695612796,-59,7594612885315268,-55,-8804246341378506,-53,6400230417475745,-59,7140850557748310,-59,6645362014883432,-55,-8851956675695884,-53 +47,-8899426359762398,-53,-5530315556905475,-55,8663552334368522,-59,-6094296256847150,-54,8475080167432525,-53,8476764036048262,-59,8486821294002839,-59,-6118040144470656,-54,8470799850456136,-53 +48,-5796945494158226,-54,-8528097039594337,-53,7596417449760120,-63,-6773719215073852,-53,5936765965120373,-53,8780036199831016,-61,7955711299114477,-62,-4593552531020411,-53,-7747817248388280,-53 +49,-8493522704481427,-53,-5992935204109977,-54,6679464149051980,-59,6112679509756143,-53,-6613744176162359,-53,4871530790417254,-58,4547026701995195,-58,-7778902563157163,-53,-4538515446331232,-53 +50,7591051118782235,-53,4847357869320683,-53,5970556321664841,-59,5591052951616009,-53,7061215514994128,-53,6071828661000874,-59,5986956807403314,-59,7528029005073837,-53,4944660435441206,-53 +51,-5631721150506652,-57,8998964833983582,-53,4996367362121957,-58,6183729103504663,-53,-6547984785204512,-53,7845690488330289,-59,4915146982706384,-58,6706819159161916,-53,-6010376255266272,-53 +52,7787473465745020,-53,4525002001414263,-53,6156095467663277,-59,-8152742403613017,-55,8773187807162609,-53,5204727404114578,-59,5275473451098260,-59,-7634871759060686,-55,8802251965038936,-53 +53,8973229237392732,-53,6189005579608570,-56,7097582124491932,-59,8221117175620286,-53,-7355481451067712,-54,8367606587161904,-59,7755887071489037,-59,8972949060955748,-53,-6202645018316924,-56 +54,4517956231675190,-53,7790974461994384,-53,8687710834303690,-59,5280467243176870,-55,-8908981405073111,-53,8339753664315080,-59,8435522168029806,-59,5420322063572168,-55,-8903708364670241,-53 +55,8787841269304898,-53,7881409313226341,-55,4658824429387542,-58,-8833479736977020,-53,-7039274037156953,-55,6155942241895291,-60,6400691366595966,-60,-8833884802165582,-53,-7030922757664524,-55 +56,8954450879136440,-53,-7741926445293201,-56,6689638320802598,-59,-8888108314554782,-53,-5830638360427896,-55,5119008991295557,-59,5434749104984182,-59,-8877254643136688,-53,-6088385475934662,-55 +57,7575611155136828,-53,-4871745283215318,-53,4895868753245753,-59,6897538650760414,-54,-8320554725576498,-53,8093015592689108,-60,8798096098969420,-60,4603097796038946,-53,-7741860537183101,-53 +58,-5725524486417363,-53,-6953184274251557,-53,9000933579801708,-61,8158865687525795,-53,7625790524377863,-54,4996390592064760,-58,7805829346334243,-60,-5303253293546639,-53,-7280207688501071,-53 +59,8382290900001888,-54,7971632531047732,-53,8347774314207672,-59,-8504980784743272,-54,7939045243000498,-53,8472330293942098,-59,8468590908555150,-59,7947955725976436,-54,8082056008100686,-53 +60,-8027918250125391,-53,-8168742171393760,-54,4574261802921349,-62,8919974043476002,-54,7825473393667678,-53,5616748152480541,-62,5257953496014380,-62,-8185033051993014,-53,-7519246475806476,-54 +61,5820591167772094,-54,8523871450462148,-53,7500324917516266,-60,4660784803532222,-54,8700432636966279,-53,4757741248429687,-60,5706282688901421,-60,5063269700014454,-54,8643984960346702,-53 +62,-7505964713467515,-53,4978862451633784,-53,8337989161484090,-61,-6673826359195730,-53,-6048726627442250,-53,6509562520760380,-60,5095289897333927,-60,-8261879612283795,-53,7174788942758614,-54 +63,-8998604286439707,-53,-5899569976376138,-57,8775837397578392,-59,5600095770999866,-58,9004945414967184,-53,6390533980014037,-59,8420396799816551,-59,-8032206506964754,-55,8779524859360546,-53 +64,-7501788238216111,-53,-4984957163068521,-53,7026868848678285,-60,-7310146700739283,-54,-8232226851924744,-53,5828015697525790,-61,6577401302011267,-60,-7099348841094434,-53,-5543125826773162,-53 +65,9006985090910802,-53,-5556606590468726,-60,5686244593179748,-60,-8495311728146021,-53,-5986358495943380,-54,7150529438114063,-62,8916086530183363,-62,-8433135015945781,-53,-6328212689844411,-54 +66,4671402075118450,-53,7700263470968168,-53,7459101805340417,-59,-5253149313494590,-53,7316589499994064,-53,5085956121176996,-60,8875252138095244,-60,-5905500328091089,-54,8509177230375622,-53 +67,7443884097249668,-53,-5069748559932126,-53,8064381041486251,-59,6125786040133485,-53,-6601977530466425,-53,8652196979899008,-59,8377596355059241,-59,6901693314983620,-53,-5786115536819773,-53 +68,-4742741215889345,-54,-8689430596551785,-53,6123664907318164,-63,-5089170977407037,-55,8916826994804424,-53,8513422331333061,-61,6067083507956773,-61,-7835534064125506,-54,-8110498025932681,-53 +69,6668015275618314,-53,6054955100455302,-53,8802669161176390,-60,-6577908973239637,-53,-6151753826001359,-53,8264880472354099,-59,4598941919297419,-59,6669384313798469,-53,6053411240789422,-53 +70,-7172352705414992,-53,5447716452539175,-53,6198602856169418,-59,-4964150137372022,-53,7514319120561293,-53,4731238969387144,-58,7668210011584211,-59,-6337260906175129,-53,6399563016346596,-53 +71,-8095465537495927,-53,-7896896295511240,-54,6813364742986367,-60,5985352978256215,-54,8494480319597814,-53,8425168175358714,-59,8347435589634612,-59,-8946609924137079,-53,-8278348198864786,-56 +72,-8907761189757353,-53,-5304355845449297,-55,4844282922191554,-58,-8994042202139525,-53,-7370361202985777,-57,5023805572246149,-58,4909792853766124,-58,-8946104753487277,-53,-8288042757769581,-56 +73,5209280030172565,-53,7347912526340498,-53,8947143600723303,-61,-8884454597275458,-53,-5913965051592043,-55,6453697640922198,-59,4837063594716773,-59,6665863224653288,-54,8367526102677054,-53 +74,8946957143843435,-53,8319984445635350,-56,4669918745094932,-66,-5990105354291072,-55,8880771141029506,-53,8811854251652667,-59,5920973445809573,-60,8096111740200266,-53,7894421660583774,-54 +75,6922341892053884,-54,8315530957081203,-53,5548469743999534,-60,8842725309009733,-53,6851399064806017,-55,5714597850637126,-60,5709070633719992,-60,7794220319139279,-54,8120355994054906,-53 +76,8001918341466138,-53,8267155687792561,-54,7148395924629846,-59,6511633517752295,-53,6221891309278716,-53,8179961307551210,-59,7264050072885312,-59,7895377656500038,-53,8667125855378765,-54 +77,-7304994962192126,-55,-8819751728513594,-53,5133342469042251,-59,6567777493194521,-53,6163521913659308,-53,8991792132590433,-60,4658469302956257,-59,6632924099818282,-53,6093328985280701,-53 +78,4723524316474427,-53,-7669276175966750,-53,6463845998783612,-62,-6513107828458196,-53,-6221308342245393,-53,8478336021385068,-60,7518784682166197,-60,-6114545085788538,-54,-8472264901790143,-53 +79,-6857529111760787,-55,8842433277647976,-53,5587533603691283,-60,8365821292283395,-53,6674724257574795,-54,8762664906609869,-60,6028288314788812,-60,-4820706652802036,-55,8926083664835095,-53 +80,-8553010906861605,-53,5644605035375929,-54,6479929794078678,-59,5213853433194723,-55,-8911935807652216,-53,5732145437385270,-59,6261672990327447,-59,6622712102526340,-56,-8968542166032174,-53 +81,-4751107609714672,-53,7651236530393555,-53,7889009652663333,-59,-8045364738339720,-53,-8097314271333173,-54,6438298051697994,-59,6561541969378192,-59,-8112659352093073,-53,-7824036160908851,-54 +82,8560686049787219,-53,5597923374255595,-54,6433908755871024,-59,8648964526118847,-53,5026491330944719,-54,5951403162806783,-59,6225595909788579,-59,8600117716138883,-53,5350944316648942,-54 +83,5142498624996478,-54,8631875253169522,-53,6086733552538698,-59,-8445733671600714,-53,6253103115759503,-54,4946719657540179,-58,7559008736648242,-59,4509402761414419,-55,8935590018507831,-53 +84,-5600113766747984,-53,7053835849659139,-53,6941572338747887,-59,-7737478556963185,-57,-8994063550636911,-53,6523785789794678,-60,8169906167833129,-60,-8074329625756491,-56,-8950245569019261,-53 +85,7811477461227682,-53,-8968798026564007,-54,6393735245487840,-61,4787237308910350,-53,-7628650578009791,-53,8016040600769051,-59,7685593939818194,-59,4999642568123023,-53,-7491247667825760,-53 +86,6022687940075693,-53,6697064323037221,-53,5038488835237120,-59,-7804756123871605,-53,8991196252355698,-54,4533508745868005,-59,5031372529924830,-59,-7399407507543064,-53,5135389625229471,-53 +87,8948455428761609,-53,-8201092914131071,-56,7971817149904061,-60,-8292357930978970,-53,7032993428094064,-54,6720764217852006,-61,5265526274930293,-60,-8112875838304429,-53,7825392463089586,-54 +88,5782564442790957,-54,8530458811332443,-53,4785897593591870,-60,5929341145282696,-53,6779237137670419,-53,7705511944893200,-59,5975890445222911,-60,6485549805940412,-54,8403086914969738,-53 +89,-7261562325331577,-53,-5326900111016299,-53,4904068351762196,-58,-5148305534069251,-53,7389235293051826,-53,4935716182015066,-58,4930578654711480,-58,-7289456568895077,-53,-5288640706150672,-53 +90,7644829625200391,-53,4762973645591137,-53,8875009674754436,-62,-8450643837799684,-53,6231523606687310,-54,5825591374531142,-59,4602774555687988,-59,4789288777977999,-53,7628052121369128,-53 +91,-6185404689624092,-56,-8972869927089630,-53,8924087810802743,-59,6536441196506080,-53,6195471853213701,-53,4604378050065115,-58,8949949406780587,-59,-6107739784716587,-56,-8973694878857009,-53 +92,7403654831726760,-54,-8210653445589951,-53,6797931112144875,-59,-8039387548320536,-53,8123200278055996,-54,4598342169211331,-60,7641394870572238,-60,-8444868228028060,-53,6264271746083061,-54 +93,8942156383259892,-53,-8574564560863662,-56,8746788339200142,-59,6371204054531448,-53,6366883104122431,-53,7183499932736176,-62,8620975097649955,-62,6480095095145888,-53,6256014902439898,-53 +94,7173432727547938,-54,8262047433434845,-53,7761010728178697,-60,6215819033320099,-55,-8871802482872044,-53,4950692097641433,-59,7968160150284816,-60,7288394625025433,-54,8236843396824359,-53 +95,8152085343625108,-53,-7660965418436697,-54,5977362845289395,-61,-8184750992859140,-53,-7520062597038576,-54,5206050727446437,-60,5042184731163374,-60,7061218334151165,-53,-5591715507517893,-53 +96,-6777669703548942,-55,8845959779914348,-53,5573357824618522,-59,9006369470848276,-53,6102516555379996,-59,4897317209223457,-59,5160410250984272,-59,8994941583047753,-53,7404438156393425,-57 +97,-4620800125316640,-55,8932656028199724,-53,6782290007492114,-60,-6777800900173649,-54,-8344814191258499,-53,6094556848812763,-59,8510655181198616,-60,-6302161446165970,-55,8868082663744220,-53 +98,-6413050379186343,-53,-6323021443439572,-53,4727243193893847,-58,8805998070606811,-55,8733765027325001,-53,8477752080714093,-60,7740932346816440,-59,7610911025804153,-56,8955998871319742,-53 +99,-8807047782919219,-53,-7529609877119647,-55,4758079088922589,-58,-7363303152870186,-53,5187505361809276,-53,8842992769044588,-61,8109968846433505,-59,-9006292163885786,-53,-8603353745339531,-62 +100,-7085014689512724,-54,8280515258269805,-53,7393921001258723,-59,5007563026318305,-54,8651013890968498,-53,4601240290722946,-58,8283182157001636,-59,-7795403808273599,-55,8792887331404302,-53 +101,5375694982337102,-65,-6481958845993197,-61,9007163570295682,-53,-6438546729547440,-66,-5568548850432336,-65,9007199117851013,-53,9007197225267388,-53,-6168870396522237,-67,-6179564384935653,-63 +102,-5660586521155651,-59,7676331989532235,-59,9005966319830368,-53,5462432849769068,-59,-5824166842316969,-60,9006679928579330,-53,9006351414547577,-53,-8917943608683560,-60,6532548534949204,-59 +103,8411990612276466,-60,-7620675012007754,-59,9006172387846469,-53,5734433733722613,-59,8360130263978877,-60,9006516768682513,-53,9006237078279749,-53,8501647995857027,-60,-7274805500247634,-59 +104,7275077652345339,-59,5117845237726663,-59,9006126927229656,-53,6029201241496458,-60,7000726101782987,-59,9006411845833016,-53,9006180432751617,-53,6860892535346576,-59,5300914633936573,-59 +105,7753316299569421,-60,-8034134624329478,-59,9006120737012831,-53,-5648708241879017,-59,-6598072048049056,-59,9006176760722601,-53,9006150934390342,-53,7305269020592439,-60,-8000384277840487,-59 +106,5659827413985561,-63,-4885374811816115,-60,9007116694356773,-53,5458821427510693,-68,-5709469698801713,-68,9007199251515154,-53,9007144817676969,-53,4622628744115507,-63,-7932982030825839,-61 +107,7476095359024225,-61,7043487239832656,-61,9007109890022628,-53,-6775348332535672,-59,-8108220346938870,-62,9006563177174396,-53,9006860634327580,-53,-4979383287293778,-59,-6993461316184382,-63 +108,-7729223868393255,-61,-7067119046956725,-60,9006979432111067,-53,-4700914645935947,-59,4727624910865242,-60,9006824028841826,-53,9006947492612091,-53,-8290608599226333,-60,4720947157859257,-61 +109,-5779710294223299,-59,-4924127703237967,-62,9006741385732338,-53,-8111842617001153,-62,5816701638308093,-62,9007178155970778,-53,9007072200107767,-53,-6116221894145106,-60,4839149580210679,-64 +110,-8535172880744023,-60,7876974508280715,-60,9006742198045782,-53,-8485089290857435,-60,8230442551675140,-59,9006037194643912,-53,9006062243072736,-53,-8486418584538077,-60,8117051966327144,-59 +111,7201112827941142,-59,-7710006133858627,-60,9006295024387887,-53,-5133753967165408,-61,-7612382191647260,-62,9007164659714083,-53,9006404269287188,-53,6712463085262782,-59,-7375587701029570,-60 +112,-4675306365535924,-59,6030929151174511,-59,9006410048431940,-53,6028210242558294,-59,5305412099769207,-60,9006611378860134,-53,9006586064172713,-53,-7179479567399340,-60,5688383435053941,-59 +113,7929932346527551,-60,-7201180116108300,-60,9006810489324313,-53,-7168755780437862,-62,7767912018103028,-63,9007185177852655,-53,9006818654843770,-53,7843300931833150,-60,-7128360228266743,-60 +114,6867499751194133,-60,-8633730186008738,-61,9006976320103166,-53,-6402234121438141,-59,-5934423502427925,-59,9006166411855353,-53,9006521349995085,-53,-4702929105286054,-59,-5282150784520287,-59 +115,-4995027232500379,-62,7794951612473149,-65,9007193770266230,-53,4967094937027289,-58,8877650897436352,-62,9005844990155571,-53,9005881773019180,-53,4899044151356701,-58,8775779480656449,-62 +116,7691497038606196,-60,6280464264737661,-59,9006464216814683,-53,5676811624250957,-62,-7599646877211590,-59,9006409674744600,-53,9006449288751163,-53,6037042741432828,-62,-7400402149093593,-59 +117,4812328981743053,-58,-8386124517378351,-62,9005928847455158,-53,-8131783742024024,-66,6211577692740582,-60,9007068472441650,-53,9006607967469504,-53,6600461578992785,-59,7953450673878321,-64 +118,-4819201547394939,-59,-4505838464431575,-60,9006815705331522,-53,-4852108756247018,-59,-5132018582260482,-60,9006790943580210,-53,9006809798716131,-53,-4827384792018633,-59,-4661552221829548,-60 +119,-6802715363049122,-59,-8522942540221323,-60,9006325927211452,-53,5767184868572632,-59,5310624123290767,-59,9006366235212665,-53,9006341700445979,-53,-6744539769453317,-59,-8434336853626021,-60 +120,7763548894852344,-60,6144783119376896,-61,9006963057247577,-53,-4595665419794685,-59,-7091690205553602,-59,9006231386816171,-53,9006694745679226,-53,-5802512712777418,-60,-5367296077465602,-59 +121,-5190400487869711,-63,8180763688476697,-65,9007197607099598,-53,8518885204850962,-63,6075039551727939,-61,9007164152087033,-53,9007174230828249,-53,6357303904617316,-63,5197790261743522,-61 +122,6479856668511110,-59,7319479446398773,-60,9006448654191161,-53,-6192876929006989,-61,-8945695726588114,-63,9007162532978552,-53,9006671974698538,-53,5413020163636779,-59,6198217135291797,-60 +123,5510164546734743,-59,-8947490186822301,-60,9006516502575567,-53,-5248781629910470,-59,-6555652360523709,-61,9006789474668235,-53,9006526491503612,-53,5454827039002300,-59,-8918334627289095,-60 +124,-6957613054901674,-61,7444471423446652,-59,9006407133265295,-53,8926372909596017,-59,5640825552955954,-60,9006011502587361,-53,9006152229390130,-53,8213254266424461,-59,6259328657198902,-60 +125,8947796589107545,-61,-6393305639294704,-61,9007096816178505,-53,-8655385026753620,-59,-7231813321360789,-60,9006006682544563,-53,9006366769717661,-53,-7096054985562713,-59,-6654237046556266,-60 +126,-5915275169579052,-60,7860691520611610,-62,9007067616779109,-53,5408421828146783,-60,7397818653264823,-62,9007088558681787,-53,9007075918269255,-53,-5705153448868594,-60,7852109177766380,-62 +127,7184518615049765,-60,-6144367328451770,-63,9007022368122521,-53,4587617988630397,-60,8986246522014642,-63,9007123671984666,-53,9007104979888204,-53,5234455691333476,-60,5217529285479565,-63 +128,-7421921719541002,-59,4813194723086114,-59,9006138683125058,-53,6381371716005741,-59,-4984940323191286,-61,9006626303339745,-53,9006540269413616,-53,-5677946621068567,-59,8095351298682186,-60 +129,8358235394413965,-59,-6366590456131747,-64,9006251887575837,-53,6261107284851803,-61,-5073907352607340,-60,9007078823123701,-53,9006568447080990,-53,6781877731550182,-59,-5932389285140882,-62 +130,-7005658174252793,-61,-6383331252839332,-59,9006605439716941,-53,-8970195671275483,-61,-6060330615973928,-62,9007123321040907,-53,9006787542482891,-53,-7422483687806702,-61,-5189892033771388,-59 +131,-8926997154563717,-59,-8721212245759812,-62,9006103063590247,-53,5887083976736358,-60,6944430386448179,-59,9006428224530683,-53,9006166021258627,-53,-8682180391430515,-59,-7395329762536729,-62 +132,-7720200412609109,-66,7853762495452410,-64,9007198389090765,-53,8566076766660211,-59,6030466521190675,-60,9006081517923357,-53,9006245897092512,-53,7904771032280450,-59,5605800681214163,-60 +133,6890382313873592,-59,7447478430796766,-60,9006367855535980,-53,4816438806124057,-61,4857714691885892,-61,9007159617344838,-53,9006633642880371,-53,5630228706696933,-59,6335298977282368,-60 +134,-6046780971992439,-60,-5210527836809567,-61,9007052374777045,-53,5776045865115069,-59,6345049520292608,-60,9006610681522515,-53,9006941285016632,-53,7578794976149063,-60,8648699353176843,-61 +135,7475554263002767,-60,-7319971712076989,-59,9006283694833073,-53,-5206010472597178,-60,-8761478491647491,-60,9006847335998952,-53,9006793596818841,-53,-8483656509021017,-64,-5464550473788398,-59 +136,5418306770898283,-60,4700224583879507,-59,9006800373036070,-53,-4755614793345645,-60,-4656682946257389,-58,9005947012409278,-53,9006212429868915,-53,-8030152934596714,-61,-8293408521421514,-59 +137,-6028269928602277,-60,-4628824329333647,-62,9007071591862106,-53,4619646302381706,-60,8308028787117589,-59,9006191451071178,-53,9006931363283020,-53,-8265891482416919,-70,8891926768968859,-60 +138,4962792858624953,-62,-6003237588722771,-59,9006705608103975,-53,-5784933452227605,-63,5233460951900748,-60,9007104684671848,-53,9006949794034668,-53,6765547778269842,-63,-8538828248662443,-60 +139,-4505458955898898,-58,-4706244065080447,-60,9006023715513064,-53,-7599430528418973,-60,4717441305210331,-60,9006928181156083,-53,9006080569690083,-53,-8819645691069221,-59,-8720652500398833,-61 +140,-6264151395149909,-60,-8464501698580320,-60,9006823545842234,-53,8697968106787914,-61,7285586689627643,-60,9006955328189944,-53,9006834130690477,-53,-6180586731263791,-60,-8340489317485662,-60 +141,7687453236689402,-60,6259707364962604,-64,9006998505697018,-53,-7543219259700429,-61,5823548342994405,-59,9006691427570185,-53,9006860330174001,-53,-7636179416299075,-62,4908832388626717,-59 +142,-6532312610771512,-62,5038643989246679,-59,9006846141205938,-53,-5552873488713754,-60,7203786027174117,-61,9007050826114357,-53,9006995507484878,-53,-7128668432819734,-61,6886974464431927,-60 +143,7446335300660685,-61,-7543786300769889,-59,9006380994505953,-53,-5179783466317497,-59,-7016916119892872,-62,9006825204182164,-53,9006650981509015,-53,4665721087702607,-62,-6333565567993904,-59 +144,-7995907278933493,-64,5241792084192042,-65,9007198317662712,-53,-7901288846956532,-59,-6903231318051515,-61,9006312756874427,-53,9006546389207335,-53,-6784813709693659,-59,-5848118642287467,-61 +145,4592765895153910,-60,8979925237931763,-61,9007059482165116,-53,6424853551637831,-62,-5887826585479626,-62,9007183172695447,-53,9007154846946312,-53,6391519318537047,-61,6804651895314007,-62 +146,-8823649296849835,-59,-8004954008370219,-60,9005926899580760,-53,7328979072598747,-60,-5799339343506281,-59,9006561439112344,-53,9006246918117121,-53,-7236361140566970,-59,-8462009805398012,-60 +147,-8894039392889116,-62,-5361887990026311,-67,9007182497869911,-53,-6789981779468948,-60,-6431981775058139,-59,9006482346954289,-53,9006797207314592,-53,-5547385076198795,-60,-4687412030011033,-59 +148,-5575895989392784,-59,-5551931993636620,-63,9006776256583474,-53,-7824776397706523,-61,6301263784969378,-63,9007145291172446,-53,9006783061999218,-53,-5531251769491280,-59,-5405735371848965,-63 +149,8720991882246274,-59,7884185482452667,-64,9006167626045399,-53,5297855581266432,-65,-6174612996782939,-66,9007199130337212,-53,9006394127440955,-53,7704572097646595,-59,6774834893448805,-64 +150,-8927408951411619,-62,-4870358174005830,-61,9007162285872849,-53,-4640539105191184,-60,8903489363123337,-63,9007122095705644,-53,9007129753985557,-53,-8903128841734642,-61,6676513779412930,-63 +151,-5117494527018330,-60,-8786262641473046,-60,9006848958601116,-53,-8382472536233694,-61,7124091267848341,-61,9007096747482153,-53,9006968561850889,-53,-4953485167642686,-60,-6599290568241850,-60 +152,5110437728258474,-65,-4684136335172193,-65,9007199095731377,-53,-6680852834165353,-59,-5178390380290196,-60,9006503471937140,-53,9006905761940080,-53,-8653513606865970,-60,-6852467101648829,-61 +153,-7574297525542693,-61,-5390546548775059,-61,9007126047101014,-53,-7382805468646492,-61,4726384045863351,-59,9006850333838560,-53,9006925863476805,-53,-7402833536863982,-61,8184689868268462,-60 +154,-6366666449680582,-60,4544436493510779,-58,9005942289097748,-53,-4603059263867430,-60,-4512277777779976,-61,9007110219898831,-53,9006085159919785,-53,-6266414408359641,-60,8507856180807128,-59 +155,-8315462892464006,-61,5323463686543271,-59,9006756605475222,-53,6490737557567408,-63,-8745715303305593,-64,9007196012109923,-53,9007160900512520,-53,-6989344518443292,-63,6498261367223945,-61 +156,-7973902878605453,-64,5552789713683767,-60,9007093944737597,-53,8606181471032152,-65,-8146674810860307,-71,9007199009622350,-53,9007135942763560,-53,-5228981361705714,-64,8620794902868581,-61 +157,-8999419439622098,-65,4618315707476009,-60,9007126721556263,-53,-4832547899485344,-60,5687690147887049,-60,9007010522512040,-53,9007075579101943,-53,-5968059317753585,-61,5253372457557127,-60 +158,-8990704857832094,-60,6224070642221082,-59,9006400335119345,-53,-5359219536704428,-60,-6498682060214756,-68,9007101940705833,-53,9006616415947807,-53,-8318733710571542,-60,5069929798731620,-59 +159,-6414108758732472,-60,6336798424388515,-59,9006515636454997,-53,-8229495029802859,-61,5277117701927083,-59,9006764468903011,-53,9006546027353101,-53,-6163059261904596,-60,6221101279609277,-59 +160,-7388382729501383,-66,-8843048243601638,-61,9007132971787952,-53,7589390320864933,-59,-8862223923569642,-60,9006152483472549,-53,9006849277914675,-53,7660263228764492,-60,-6679317902344886,-60 +161,8160789164630476,-61,8523520901804055,-61,9007081305572417,-53,-4690259926733629,-59,7721594911883860,-60,9006699094391751,-53,9007000215223470,-53,-8446367990395229,-61,6396105417106366,-60 +162,8747427571126986,-65,-4953611146748571,-60,9007115862309567,-53,-7754552956921649,-59,5700438183127324,-59,9005943822584274,-53,9006044776618213,-53,-7476080969583626,-59,5411857720320665,-59 +163,6492221696153478,-59,6277121132933159,-59,9006093962176750,-53,7266319797227956,-60,-4695934538896716,-60,9006945645424928,-53,9006395374995129,-53,6008540994690608,-59,4817742618598352,-59 +164,4609208611145705,-58,-5405664948661829,-61,9006022743439323,-53,4563686075947026,-61,8483534944161341,-60,9006937764512752,-53,9006349423684022,-53,7906171855576214,-59,-7082642880723324,-63 +165,6208922480844416,-60,4620203329427660,-59,9006779334014922,-53,-5945468894124777,-62,-8029372916201887,-59,9006317983604999,-53,9006563068418619,-53,-6105294071786963,-63,-6840687343675955,-59 +166,8575810275914064,-60,-6557618938088597,-61,9006913647259410,-53,-6533955133969470,-60,-8651912640070172,-59,9006040049647139,-53,9006741286213762,-53,-6541189855306114,-64,-5809429356372946,-59 +167,-4698624813098136,-59,5219094690074418,-59,9006530873180766,-53,-4821897936011026,-59,5127079757704772,-59,9006527869418804,-53,9006528770472514,-53,-4787563613492514,-59,5152708037301813,-59 +168,5155333800591746,-62,-7697865293356140,-64,9007192842479658,-53,-6504789471968069,-59,-7612159529789327,-60,9006429458235720,-53,9006952290598148,-53,-7138021638155764,-60,-4683865628600057,-60 +169,-6164954155820854,-62,4757731895523800,-61,9007172033027975,-53,6139905764853507,-59,7964970845479474,-61,9006634591457353,-53,9006885212026956,-53,8931432218545501,-60,7188032240803105,-61 +170,-5162378933717069,-60,-4957132723938446,-62,9007103756423083,-53,-4597041127479358,-58,5434833954344559,-61,9006028548806312,-53,9006566339953821,-53,-6803322473912666,-59,5147618840866498,-62 +171,8590520187773888,-65,-7409887096836197,-60,9007012978430883,-53,-5975980962660514,-61,4695033365630232,-58,9005973950098437,-53,9006508705002440,-53,-4805331551329146,-61,7036222229814210,-59 +172,8310893488776074,-60,6627124496928198,-60,9006816422510289,-53,8402484552129729,-62,4653722557763948,-62,9007179718141983,-53,9006871990473364,-53,7708248774250784,-60,6096926201552450,-60 +173,-4878276846857441,-58,6813650013894514,-61,9005869762750104,-53,7023578845738105,-62,6104755620257763,-59,9006683717673596,-53,9006657329879037,-53,-5208208376751296,-59,7172268685219563,-60 +174,-5770895330741121,-66,7701190259205803,-62,9007186668156722,-53,6849284972348615,-61,6555478310492472,-60,9007013913649663,-53,9007114561677865,-53,7778956658187616,-62,4605958339585340,-60 +175,6678434395569710,-59,-5037151577906936,-60,9006508798711717,-53,-8664320661649140,-65,7284494127150363,-62,9007187769638071,-53,9007084507124563,-53,5698106064600644,-60,-4730987485120705,-62 +176,-7292238725344689,-62,-5074925114024863,-60,9007100732711741,-53,6836167896697542,-61,5986957608493432,-60,9007038225699983,-53,9007094811441097,-53,5682908573817237,-61,4769924526254457,-60 +177,5443201362097866,-61,-6128291803431457,-59,9006665164466968,-53,-6366055187770099,-60,-5324175773371980,-59,9006677758561009,-53,9006676237745146,-53,4749536552392075,-61,-6097609460704831,-59 +178,6330092338869096,-59,-6331379421158845,-59,9006112866970752,-53,-8294173803017467,-62,-8531109241664796,-61,9007123039957013,-53,9006827067100938,-53,5781059974777635,-60,-8742281069166161,-60 +179,-5536336973342860,-59,-6408349093207122,-59,9006227242622919,-53,-6324914591348895,-60,8754403491585249,-60,9006804040402272,-53,9006505203551929,-53,-5206511556293157,-59,-4909399451755012,-59 +180,-7555769127569884,-61,7492429965869211,-60,9006960696886621,-53,7877259712007820,-61,4539768996156349,-61,9007129238112941,-53,9007059451156708,-53,-7251166010470074,-62,6162483410124068,-60 +181,7462735464620299,-62,5615348908466090,-60,9007080625573690,-53,-4933879030954202,-64,4590940297257315,-64,9007198653613665,-53,9007094664182257,-53,6915724235377232,-62,5280180155859283,-60 +182,-5239033902347196,-60,7983628617416009,-62,9007092761329195,-53,-4598273432882495,-58,-5570318816696662,-60,9005947814357080,-53,9006789665641592,-53,-5460494311495868,-59,-5089005460286755,-62 +183,-7800893318114347,-60,4788049800596429,-59,9006682361363427,-53,-6957290560164989,-60,-6205083917063781,-60,9006904797699606,-53,9006840220295641,-53,-7656882149958749,-60,6880284019042596,-60 +184,-6542577256368048,-61,-5037103988119827,-59,9006819128841109,-53,5585735356853252,-60,-5634209902492130,-60,9006985987071486,-53,9006844933356385,-53,-5593471648771855,-61,-4918165975456289,-59 +185,-6633280945402677,-60,-8866510238670515,-59,9005984661666361,-53,4971685774929807,-59,-7067590479946738,-59,9006187251087900,-53,9006134605841705,-53,-8469714454156002,-61,-8606339316901698,-59 +186,6758267896074501,-60,-8383662947730820,-59,9006091886238522,-53,4538840348482208,-58,-6861083352598396,-60,9005922883480801,-53,9005942273565106,-53,8958280328298282,-59,-7068686469469442,-60 +187,-8416006427317507,-59,7991879056631839,-64,9006238444010272,-53,-6949159907650955,-60,5285759203979262,-60,9006940973501619,-53,9006545738467532,-53,-6871980335697819,-59,7980483032077355,-62 +188,7338935160514948,-59,-5406086785492487,-59,9006073162551484,-53,7683623126529155,-61,8579344089510201,-60,9006899858647804,-53,9006200095451636,-53,7053684174160568,-59,-4895548309830255,-59 +189,6716026118662904,-64,5354487114484099,-59,9006810091207109,-53,8123553346532967,-59,4631137106128926,-59,9006014149902928,-53,9006312454843158,-53,6527668924918876,-59,4777091195948211,-59 +190,-6690475768114512,-60,6836519778559701,-60,9006889233758982,-53,-5071067869680602,-65,-5476456147925541,-65,9007199070420992,-53,9007141139761266,-53,-5979192146891167,-61,5732285585625189,-61 +191,8295796559123905,-60,7456994479874134,-61,9006918977572099,-53,-5433252251943322,-59,-6814202458563058,-59,9006169832889560,-53,9006419787944535,-53,-4581857502056882,-59,-6043061414805360,-59 +192,-4674935466885815,-59,-8126709076912704,-61,9006847116147421,-53,5412011638276893,-59,8114139257943184,-59,9005909922398793,-53,9006601484847911,-53,6309851808221992,-60,5844023601180586,-59 +193,8232111666614510,-60,4543957350009401,-58,9005850242106777,-53,4841846582929156,-63,6206997479169191,-64,9007197503755226,-53,9007134919210769,-53,8322279179382123,-62,7657586692511269,-61 +194,8990157465385659,-60,-4525994897943886,-59,9006647780614210,-53,8740015695230166,-64,-4942577223098402,-59,9006867162004576,-53,9006665544608401,-53,8659589505130494,-60,-4542308616901593,-59 +195,-6624765306170238,-60,-5541830106535675,-60,9006946498608608,-53,8950628648173209,-62,7703342212092813,-59,9006378025520242,-53,9006750211842236,-53,4669416979759448,-63,5748697423291352,-59 +196,-5836235649558595,-61,-8144528286744993,-60,9006945653732646,-53,-6396830207104673,-66,-7252622838189219,-63,9007196436243489,-53,9007005443443206,-53,-5046174532833258,-61,-7129965901632837,-60 +197,4698236761670729,-58,5550588788434406,-63,9006000941046816,-53,-8173304519908017,-62,-7825819218744415,-60,9006977605080288,-53,9006633303666847,-53,6402222139495312,-59,-7020236883518904,-62 +198,-4569294824456030,-63,-8587798887831914,-60,9006948270246841,-53,-7426844061585687,-63,4639458051274363,-65,9007196263480784,-53,9007141527686624,-53,-6051658475860880,-63,-8115599031388558,-61 +199,7363958090631146,-63,8041623648196751,-59,9006319930111034,-53,-7015201662523632,-59,6994773051430826,-59,9005869112806789,-53,9006075572976938,-53,-5579023156096196,-59,7195999475589308,-59 +200,-5225951668749630,-60,-7453126560787327,-60,9006918510896090,-53,-7603878805817749,-61,8291846096725021,-59,9006218426665620,-53,9006497128540788,-53,-7929355510828203,-61,6919254483928827,-59 diff --git a/test/grid/geometry/data/accusphgeom/gca_gca_pairs_seed20251104_N100_with_baseline.csv b/test/grid/geometry/data/accusphgeom/gca_gca_pairs_seed20251104_N100_with_baseline.csv new file mode 100644 index 000000000..8941026a4 --- /dev/null +++ b/test/grid/geometry/data/accusphgeom/gca_gca_pairs_seed20251104_N100_with_baseline.csv @@ -0,0 +1,32 @@ +pairs_id,ref_angle_deg,a0_x_sig,a0_x_exp,a0_y_sig,a0_y_exp,a0_z_sig,a0_z_exp,a1_x_sig,a1_x_exp,a1_y_sig,a1_y_exp,a1_z_sig,a1_z_exp,b0_x_sig,b0_x_exp,b0_y_sig,b0_y_exp,b0_z_sig,b0_z_exp,b1_x_sig,b1_x_exp,b1_y_sig,b1_y_exp,b1_z_sig,b1_z_exp,baseline_x_sig,baseline_x_exp,baseline_y_sig,baseline_y_exp,baseline_z_sig,baseline_z_exp +8,5.85151e-08,-7230928401775691,-53,6413551006731969,-54,-8616241346204776,-54,8773798999433595,-53,-6000682488046241,-57,-8009426984445850,-55,-7230928402586527,-53,6413550996183423,-54,-8616241351334767,-54,8773799000690885,-53,-6000682357192868,-57,-8009426968536627,-55,-8679261493192915,-55,8181549207548938,-54,-7725742864056452,-53 +14,8.70857e-08,-4607551742761558,-53,-7277539273274117,-54,-6830773080982109,-53,6068630138184276,-53,5502869841781176,-53,7488602049766939,-54,-4607551734638633,-53,-7277539289680749,-54,-6830773082091328,-53,6068630128670322,-53,5502869851389287,-53,7488602052365278,-54,5399770290349084,-55,4985477006629652,-54,-8549476488032171,-53 +21,2.51618e-07,-7848987668729198,-53,7391109489544601,-54,4843926555650740,-54,-4699941368622357,-56,-8962842026919421,-53,5378199410351573,-56,-7848987668539763,-53,7391109489604955,-54,4843926556786475,-54,-4699941463474771,-56,-8962842028808165,-53,5378199126013745,-56,-7899292735369896,-53,7141700366785971,-54,4890742955135284,-54 +24,4.02326e-07,-4568328231438871,-56,-5420625993444073,-54,8570749909878658,-53,-6116375928311063,-53,6025041063082346,-53,-5447371515274217,-54,-4568327936189426,-56,-5420625903420298,-54,8570749926571625,-53,-6116375934808137,-53,6025041055158319,-53,-5447371521151577,-54,-6650054535212328,-53,5934449427713896,-53,-5198035141420018,-55 +27,5.38487e-07,-5234064471250982,-55,5058579375404553,-53,-7336770317160175,-53,7902251572577888,-53,-6669376776424640,-54,5500513376701818,-54,-5234064502419141,-55,5058579335254894,-53,-7336770343452975,-53,7902251586143424,-53,-6669376636627692,-54,5500513468250617,-54,7959464020805770,-54,7744128462582881,-54,-7092142844636924,-53 +29,6.11719e-07,6050986485628052,-54,-5624975274774618,-54,-8004120335292124,-53,-8271075057631019,-53,-8456986651598985,-55,5744185531146434,-54,6050986468807466,-54,-5624975236079039,-54,-8004120345269573,-53,-8271075024522282,-53,-8456986956264200,-55,5744185609702737,-54,8263033716329409,-56,-7484661030342715,-54,-8127592586860685,-53 +33,8.91218e-07,-6217993416889531,-53,-5375584469257136,-54,-5936494688416347,-53,8468705800669433,-53,6650362224383374,-55,-5156143178421444,-54,-6217993403770544,-53,-5375584572815451,-54,-5936494678714016,-53,8468705767162147,-53,6650362753381212,-55,-5156143227983011,-54,-6704166870089359,-54,-4705053759241793,-54,-8022393180327608,-53 +34,9.48276e-07,-4962019623479789,-53,-5584302200557263,-54,6979390510927226,-53,-7402494061160030,-55,7122803430921892,-56,-8769954796524428,-53,-4962019616170019,-53,-5584302245707275,-54,6979390507092851,-53,-7402494213353407,-55,7122804370969596,-56,-8769954776565984,-53,-6025832852997469,-53,-5940435175330124,-54,5999732021037956,-53 +35,9.95453e-07,-6035617260517938,-53,6276240380711044,-53,-4608587368240598,-54,7967925325362540,-54,-5818023032045022,-53,5604309762634636,-53,-6035617359078084,-53,6276240269121175,-53,-4608587459803159,-54,7967925495800519,-54,-5818022935560021,-53,5604309802218937,-53,-5796477142969564,-54,-7348950478093914,-56,8478560672644647,-53 +37,1.58465e-06,7109237221237507,-53,-5391249802838377,-53,4936087698042736,-55,-7325286922064621,-54,6319459771478421,-54,-7598077937289827,-53,7109237104464673,-53,-5391249958781681,-53,4936087663801729,-55,-7325286912718776,-54,6319459783959252,-54,-7598077936947270,-53,-6932332681741704,-54,6041380115989781,-54,-7745370287181335,-53 +42,6.23097e-06,-7588425694061722,-54,5300337922171494,-54,7727237007089296,-53,5408937037834848,-53,7294951169381972,-54,-6210391323775678,-53,-7588426574445265,-54,5300338281723904,-54,7727236729290293,-53,5408937802886953,-53,7294950544480860,-54,-6210390840961762,-53,-7744496150009437,-56,6740256028889544,-53,5895883564613215,-53 +44,8.33965e-06,4682519649110699,-54,-7138809952082747,-53,-4968453815481641,-53,-5285473185278374,-56,4695891434211099,-53,7657789400935486,-53,4682518032116495,-54,-7138810329968187,-53,-4968453653509528,-53,-5285464430986344,-56,4695891945672347,-53,7657789181709099,-53,8009098880075897,-54,-7929304319014417,-53,5959000202527590,-55 +45,9.41793e-06,8237100306028881,-53,7022501152606908,-54,-7801277242848944,-56,-7777560375867427,-54,-7705965765779950,-53,-5147861179623380,-54,8237100704763774,-53,7022499948427653,-54,-7801267641690998,-56,-7777560736256870,-54,-7705965493685948,-53,-5147862264352867,-54,-7024749570822305,-58,-8077493784355160,-53,-7958779314590154,-54 +46,1.51358e-05,-8035181008987649,-55,7196128282874947,-53,-5030916197430068,-53,-7477428736703008,-57,-8558733024802324,-54,7911905879718196,-53,-8035186374763230,-55,7196127591634174,-53,-5030916650542349,-53,-7477422757881264,-57,-8558732639694402,-54,7911906005938418,-53,-5306357025023269,-55,-5877246113821485,-54,8410368237430861,-53 +47,2.93867e-05,-8530315408545074,-55,-8548698092705777,-53,-7484925409163612,-55,-8505865213578574,-56,8942315233371556,-53,-5915409831234063,-58,-8530307527993322,-55,-8548697921080605,-53,-7484937526617150,-55,-8505884745109989,-56,8942315020689593,-53,-5915289701358967,-58,-7001500447725934,-53,-6239928543688813,-54,-4730164144880664,-53 +51,6.56984e-05,8275974108886025,-53,5167277393821351,-54,-4883728922649398,-54,-7782304739439814,-53,-8433407601596845,-54,6675064540130152,-55,8275976163428860,-53,5167273828800138,-54,-4883718768116487,-54,-7782306117123403,-53,-8433405211054708,-54,6675050921788293,-55,-6020764263020120,-53,-6698722592085563,-53,5412348747689106,-59 +55,0.000159124,4893062498398695,-53,-8738676079592702,-54,-6172233217836116,-53,-5954009323050466,-54,7366386763692950,-54,7661668758944494,-53,4893069065022559,-53,-8738652712226967,-54,-6172236283009425,-53,-5954016406967466,-54,7366374159614903,-54,7661670412262375,-53,-5356451187211240,-56,5259089675872662,-54,8588761644687132,-53 +56,0.000232053,8489692534713884,-53,-8426129931190372,-56,5637513389912148,-54,-7989824472859641,-54,-8510125299767764,-55,-7787420952999829,-53,8489684044474492,-53,-8426392313840217,-56,5637540021519676,-54,-7989814084509932,-54,-8510045039288690,-55,-7787429099367662,-53,6377195679142295,-55,-7371982033656451,-54,-8062339585556965,-53 +59,0.00057029,5781018189649026,-53,-7080086847675367,-55,-6676562756755131,-53,-4910678895854829,-53,5937140050511939,-54,6942799377885296,-53,5781038845985583,-53,-7080167487032430,-55,-6676539526433992,-53,-4910706842040514,-53,5937194599449982,-54,6942767949327463,-53,7224228704922008,-53,5697803458625058,-55,-5187590747906738,-53 +61,0.000767608,8251478203846825,-54,7967596070464569,-53,6326215350792608,-56,-6052848968698497,-53,5425533121280526,-56,6635715023714418,-53,8251479346515285,-54,7967595719333065,-53,6326219807145638,-56,-6052926909053857,-53,5425916326945822,-56,6635639032788121,-53,8173101628292853,-54,7982689126669748,-53,6723880617402731,-56 +70,0.00723011,4808962270786330,-56,6097550991378913,-57,-8979034354739469,-53,5423237434572450,-53,8165150607383119,-56,7118737270293392,-53,4807007705752224,-56,6114984972798435,-57,-8979004394930222,-53,5423388220981957,-53,8159770776237459,-56,7118718780114812,-53,8506946025451435,-53,7375198673046507,-55,4631159905674187,-54 +71,0.00890077,5769443369962036,-54,-5983614074820145,-53,6083122702890450,-53,-5790623327371145,-54,6593948792725823,-53,-5409865811981166,-53,5767188889012695,-54,-5983906168917601,-53,6083369827308603,-53,-5788239974269031,-54,6594257583702071,-53,-5410127062691128,-53,5211102010518357,-55,4530259865131377,-54,8619837631472525,-53 +74,0.0206073,-8035659862333313,-54,-7363807833259551,-53,-6560797284980930,-54,-7917940611768033,-56,8425555425553090,-53,6053119466146400,-54,-8035005826375428,-54,-7362865851567044,-53,-6565825221238476,-54,-7921074297809246,-56,8424427095051799,-53,6059142060728708,-54,-8880049322165750,-53,-4836184681008166,-56,-5526538576809171,-55 +76,0.0462026,8127173703048925,-57,7441123885394209,-56,8944633338416150,-53,8639627923593602,-57,8409091324007294,-55,-8741768269249357,-53,8062154403355115,-57,7452323460812561,-56,8944717496730583,-53,8731602008241724,-57,8401170055227091,-55,-8741887316731750,-53,8484083315109138,-55,5760563192395373,-53,6591420238834665,-53 +77,0.052799,5221222816770053,-53,-4950274785268426,-54,6909570579987614,-53,-4968381977795672,-57,5287825588365527,-54,-8604820653289657,-53,5221435213196362,-53,-4944297667067441,-54,6910479930790452,-53,-4975824390786103,-57,5274735625260514,-54,-8606812142843384,-53,7594585201434241,-53,-7582744394822946,-55,8912527867656076,-54 +81,0.099423,5136209688356046,-53,-7570367369932403,-54,6357780502888556,-53,-6641817027751932,-53,4674075341242368,-54,-5617308956303604,-53,5139981731784467,-53,-7586164797614735,-54,6350019852519077,-53,-6646291071864332,-53,4692812768224429,-54,-5608104000752893,-53,-6739521340379483,-55,-6597567917735689,-53,5896008331225631,-53 +82,0.173987,6554513973516502,-54,-7987770620112503,-54,-7378218293563006,-53,6401872806519504,-56,6043074187390434,-55,8843463246735159,-53,6528975833162550,-54,-8024627229178541,-54,-7373879339620572,-53,6524199114179981,-56,6131344766951964,-55,8838267408201787,-53,7404977277993908,-53,-5127941779437070,-53,6537454199100429,-62 +84,0.356562,-7451089160378882,-53,6599354785688124,-56,4993036832193222,-53,8396296442464927,-53,5695352472672690,-54,6353058498534731,-55,-7441483107366564,-53,6280708453029850,-56,5013741470761668,-53,8395252491023875,-53,5704009808254162,-54,6344058074638568,-55,8125490787250087,-53,6291989587332959,-54,4564539249100196,-54 +93,3.5037,8495418835663911,-53,5453966742131135,-54,4933248690565520,-55,-8570297643620917,-53,-8511367560966975,-55,-7101465318633758,-55,8572736635333399,-53,5115703395540285,-54,8372219540489569,-56,-8634536639855365,-53,-7949279712094154,-55,-6480709755306040,-55,-8040774305772088,-53,5865117273986101,-56,-7984508305031541,-54 +94,4.49269,-7368955444388499,-53,4662451496900326,-53,4512063879102644,-54,-7641843075768684,-55,-7528642108591652,-53,-4560627814860603,-53,-7383904688886592,-53,4545090976327561,-53,4878416946900325,-54,-7448872215335746,-55,-7149907908367128,-53,-5151756907740391,-53,-8600420037733444,-53,5033249332855021,-54,7284015936230180,-56 +99,9.5294,-7742073097949454,-55,7446055199367508,-56,-8747405273392888,-53,8198771485748285,-54,-7182485598965534,-53,7137667844998008,-54,-5584749836427226,-56,6228565086417920,-55,-8844072287110214,-53,7613143603838169,-54,-7330708031740297,-53,7183416815087230,-54,8200993966311024,-54,-7840059179852882,-53,6751633831473711,-55 diff --git a/test/grid/geometry/test_accusphgeom_baseline.py b/test/grid/geometry/test_accusphgeom_baseline.py new file mode 100644 index 000000000..94565c2c6 --- /dev/null +++ b/test/grid/geometry/test_accusphgeom_baseline.py @@ -0,0 +1,135 @@ +"""Baseline regression tests ported from the AccuSphGeom C++ test suite. + +Test data and expected results are taken directly from: + https://github.com/hongyuchen1030/AccuSphGeom + +Specific C++ tests mirrored here: + tests/test_gca_gca_intersection_baseline.cpp — 31 near-tangent GCA pairs + tests/test_gca_constlat_intersection_baseline.cpp — 200 arc/latitude cases + +The C++ library uses ultra-tight tolerances (3–100 ULP) backed by Shewchuk +adaptive precision and a geogram fallback. This Python port implements only +the EFT tier, so the tolerances here reflect what double-precision EFT can +achieve: + + GCA-GCA intersection: 3e-8 (C++ reference: 1e-8) + GCA-const-lat intersection: 1e-13 (C++ reference: 3–100 ULP ≈ 7e-16–2e-14) +""" + +import math +import os + +import numpy as np +import pytest + +from uxarray.grid.intersections import gca_const_lat_intersection, gca_gca_intersection + +_DATA_DIR = os.path.join(os.path.dirname(__file__), "data", "accusphgeom") +_GCA_GCA_CSV = os.path.join( + _DATA_DIR, "gca_gca_pairs_seed20251104_N100_with_baseline.csv" +) +_GCA_CONSTLAT_CSV = os.path.join( + _DATA_DIR, "gca_constlat_cases_with_baseline.csv" +) + + +def _sigexp(sig, exp): + return math.ldexp(int(sig), int(exp)) + + +def _parse_vec3(fields, start): + return np.array( + [ + _sigexp(fields[start], fields[start + 1]), + _sigexp(fields[start + 2], fields[start + 3]), + _sigexp(fields[start + 4], fields[start + 5]), + ] + ) + + +def _parse_scalar(fields, start): + return _sigexp(fields[start], fields[start + 1]) + + +# ── GCA-GCA intersection ────────────────────────────────────────────────────── + + +def _load_gca_gca(): + rows = [] + with open(_GCA_GCA_CSV) as f: + next(f) + for line in f: + line = line.strip() + if not line: + continue + fields = line.split(",") + assert len(fields) == 32 + pair_id = int(fields[0]) + a0 = _parse_vec3(fields, 2) + a1 = _parse_vec3(fields, 8) + b0 = _parse_vec3(fields, 14) + b1 = _parse_vec3(fields, 20) + baseline = _parse_vec3(fields, 26) + rows.append((pair_id, a0, a1, b0, b1, baseline)) + return rows + + +@pytest.fixture(scope="module") +def gca_gca_rows(): + return _load_gca_gca() + + +def test_gca_gca_row_count(gca_gca_rows): + assert len(gca_gca_rows) == 31 + + +@pytest.mark.parametrize("idx", range(31)) +def test_gca_gca_intersection_baseline(gca_gca_rows, idx): + pair_id, a0, a1, b0, b1, baseline = gca_gca_rows[idx] + result = gca_gca_intersection(np.stack([a0, a1]), np.stack([b0, b1])) + assert result.shape[0] >= 1, f"pair_id={pair_id}: expected intersection, got none" + err = float(np.linalg.norm(result[0] - baseline)) + assert err < 1e-15, f"pair_id={pair_id}: err={err:.3e} ≥ 1e-15" + + +# ── GCA-const-lat intersection ──────────────────────────────────────────────── + + +def _load_gca_constlat(): + rows = [] + with open(_GCA_CONSTLAT_CSV) as f: + next(f) + for line in f: + line = line.strip() + if not line: + continue + fields = line.split(",") + assert len(fields) == 19 + case_id = int(fields[0]) + a0 = _parse_vec3(fields, 1) + a1 = _parse_vec3(fields, 7) + z0 = _parse_scalar(fields, 13) + bx = _parse_scalar(fields, 15) + by = _parse_scalar(fields, 17) + rows.append((case_id, a0, a1, z0, bx, by)) + return rows + + +@pytest.fixture(scope="module") +def gca_constlat_rows(): + return _load_gca_constlat() + + +def test_gca_constlat_row_count(gca_constlat_rows): + assert len(gca_constlat_rows) == 200 + + +@pytest.mark.parametrize("idx", range(200)) +def test_gca_constlat_intersection_baseline(gca_constlat_rows, idx): + case_id, a0, a1, z0, bx, by = gca_constlat_rows[idx] + result = gca_const_lat_intersection(np.stack([a0, a1]), z0) + assert not np.all(np.isnan(result[0])), f"case_id={case_id}: no intersection returned" + dx = result[0, 0] - bx + dy = result[0, 1] - by + err = math.sqrt(dx * dx + dy * dy) + assert err < 5e-15, f"case_id={case_id}: err_xy={err:.3e} ≥ 5e-15" diff --git a/test/grid/geometry/test_eft_accuracy.py b/test/grid/geometry/test_eft_accuracy.py new file mode 100644 index 000000000..58f831b4d --- /dev/null +++ b/test/grid/geometry/test_eft_accuracy.py @@ -0,0 +1,170 @@ +"""Accuracy regression tests for the compensated EFT primitives in +``uxarray.utils.computing`` — specifically the compensated sum-of-squares used by +the GCA / constant-latitude intersection kernel. + +An earlier revision computed the naive ``Σ h² + Σ l²`` instead of the +correct compensated squared norm ``Σ (h + l)² = Σ h² + 2·Σ h·l`` (matching +AccuSphGeom's ``numeric::sum_of_squares_c``. The correct compensated result +reaches ~1e-30 relative accuracy on well-conditioned inputs; the naive form +only reaches ~1e-15. +""" + +from fractions import Fraction + +import numpy as np +import pytest + +from uxarray.utils.computing import ( + _HAS_FMA, + _sum_of_squares_c, + _two_prod_veltkamp, + accucross, + two_prod, +) + +# Correct compensated result: ~1e-30 rel err. Naive Σh²+Σl²: ~1e-15. +_SUM_SQ_REL_TOL = 1e-24 + + +def _unit(v): + return v / np.linalg.norm(v) + + +def _make_normals(n, seed): + """Compensated cross products (hi, lo) of well-separated arc pairs. + + ``|dot(a, b)| < 0.999`` keeps ``|n|`` away from the noise floor so the + compensated squared norm is meaningful. (Extreme near-parallel arcs lose + accuracy for *any* algorithm and are out of scope for this primitive test.) + """ + rng = np.random.default_rng(seed) + hi = np.empty((n, 3)) + lo = np.empty((n, 3)) + i = 0 + while i < n: + a = _unit(rng.standard_normal(3)) + b = _unit(rng.standard_normal(3)) + if abs(np.dot(a, b)) > 0.999: + continue + xh, yh, zh, xl, yl, zl = accucross(a[0], a[1], a[2], b[0], b[1], b[2]) + hi[i] = (xh, yh, zh) + lo[i] = (xl, yl, zl) + i += 1 + return hi, lo + + +def _rel_err_vs_exact(res, hi_tup, lo_tup): + """|(res_hi + res_lo) − Σ (h + l)²| / Σ (h + l)², computed exactly.""" + true = sum((Fraction(h) + Fraction(l)) ** 2 for h, l in zip(hi_tup, lo_tup)) + if true == 0: + return 0.0 + return abs(float((Fraction(res[0]) + Fraction(res[1]) - true) / true)) + + +@pytest.fixture(scope="module") +def normals(): + return _make_normals(2000, seed=20260716) + + +@pytest.mark.parametrize("ncomp", [2, 3]) +def test_sum_of_squares_compensated_accuracy(normals, ncomp): + """The generic keeps the 2·Σh·l cross term for each tuple length used by the + const-lat kernel (N=2 -> denominator, N=3 -> |n|²). A naive Σh²+Σl² would + regress to ~1e-15 and fail this bound.""" + hi, lo = normals + worst = 0.0 + for i in range(hi.shape[0]): + h = tuple(float(x) for x in hi[i, :ncomp]) + lo_t = tuple(float(x) for x in lo[i, :ncomp]) + worst = max(worst, _rel_err_vs_exact(_sum_of_squares_c(h, lo_t), h, lo_t)) + assert worst < _SUM_SQ_REL_TOL, ( + f"N={ncomp}: _sum_of_squares_c max rel err {worst:.2e} ≥ " + f"{_SUM_SQ_REL_TOL:.0e} (naive Σh²+Σl² regressed the cross term?)" + ) + + +def test_sum_of_squares_keeps_cross_term(normals): + """Directly assert the result tracks the compensated ``Σh² + 2·Σh·l`` and + NOT the naive ``Σh² + Σl²``, on cases where the two formulas are exactly + distinguishable.""" + hi, lo = normals + checked = 0 + for i in range(hi.shape[0]): + h = (float(hi[i, 0]), float(hi[i, 1])) + lo_t = (float(lo[i, 0]), float(lo[i, 1])) + big_h = [Fraction(x) for x in h] + big_l = [Fraction(x) for x in lo_t] + correct = sum(big_h[k] * big_h[k] for k in range(2)) + 2 * sum( + big_h[k] * big_l[k] for k in range(2) + ) + naive = sum(big_h[k] * big_h[k] + big_l[k] * big_l[k] for k in range(2)) + # The two formulas differ by ~2·Σh·l ≈ 1e-16 relative — that IS the bug + # we test. Skip only cases where the gap falls below the compensated + # accuracy floor (so the two are genuinely indistinguishable there). + if correct == 0 or abs(float((naive - correct) / correct)) < 1e-20: + continue + res = _sum_of_squares_c(h, lo_t) + val = Fraction(res[0]) + Fraction(res[1]) + assert abs(val - correct) < abs(val - naive) + checked += 1 + assert checked > 100, "test did not exercise enough distinguishable cases" + + +def test_sum_of_squares_generic_any_n(): + """The single generic works for any tuple length — e.g. N=4, which no + hand-written ``_sum_sq_cN`` exists for — with the same compensated accuracy. + This is the point of one N-generic primitive instead of per-N helpers. + """ + rng = np.random.default_rng(7) + worst = 0.0 + for _ in range(2000): + h = tuple(float(x) for x in rng.standard_normal(4)) + # residual-scale low parts (~1 ulp of each high part) + lo_t = tuple(x * 1e-16 * float(rng.standard_normal()) for x in h) + worst = max(worst, _rel_err_vs_exact(_sum_of_squares_c(h, lo_t), h, lo_t)) + assert worst < _SUM_SQ_REL_TOL, ( + f"N=4: _sum_of_squares_c max rel err {worst:.2e} ≥ {_SUM_SQ_REL_TOL:.0e}" + ) + + +def _random_operands(rng, n): + for _ in range(n): + yield ( + float(rng.standard_normal() * rng.integers(1, 1 << 20)), + float(rng.standard_normal() * rng.integers(1, 1 << 20)), + ) + + +def test_two_prod_error_term_is_the_exact_residual(): + """``two_prod`` must return the EXACT rounding error of ``a * b``. + + This is the property ``computing._validate_fma`` exists to guarantee at + import time, asserted here directly against exact rational arithmetic so it + cannot silently lapse. A non-fused FMA lowering (or a non-compliant FMA) + yields ``e = 0.0``, which fails this outright. + + Note this is deliberately checked as ``p + e == a*b`` **exactly, over the + rationals** — NOT as the float expression ``p + e``, which rounds straight + back to ``p`` (since ``|e| <= ulp(p)/2``) and would make the assertion + vacuously true for any ``e`` whatsoever. + """ + rng = np.random.default_rng(20260716) + for a, b in _random_operands(rng, 20000): + p, e = two_prod(a, b) + assert Fraction(p) + Fraction(e) == Fraction(a) * Fraction(b), ( + f"two_prod({a!r}, {b!r}) = ({p!r}, {e!r}) is not an exact " + f"error-free transform" + ) + + +@pytest.mark.skipif(not _HAS_FMA, reason="no FMA path on this toolchain") +def test_two_prod_fma_matches_veltkamp_bit_for_bit(): + """The FMA and Veltkamp paths must agree bit-for-bit, so ``_HAS_FMA`` can + never change results. The exact residual is unique and representable, so + two correct implementations have no freedom to differ.""" + rng = np.random.default_rng(20260717) + for a, b in _random_operands(rng, 20000): + p, e = two_prod(a, b) + pv, ev = _two_prod_veltkamp(a, b) + assert np.float64(p).view(np.int64) == np.float64(pv).view(np.int64) + assert np.float64(e).view(np.int64) == np.float64(ev).view(np.int64) diff --git a/test/grid/geometry/test_intersections.py b/test/grid/geometry/test_intersections.py index 94b5be7a9..8407d1c8b 100644 --- a/test/grid/geometry/test_intersections.py +++ b/test/grid/geometry/test_intersections.py @@ -2,9 +2,9 @@ import pytest import uxarray as ux from uxarray.constants import ERROR_TOLERANCE -from uxarray.grid.arcs import extreme_gca_z +from uxarray.grid.arcs import extreme_gca_z, on_minor_arc from uxarray.grid.coordinates import _lonlat_rad_to_xyz, _xyz_to_lonlat_rad,_xyz_to_lonlat_rad_scalar -from uxarray.grid.intersections import gca_gca_intersection, gca_const_lat_intersection, _gca_gca_intersection_cartesian, get_number_of_intersections +from uxarray.grid.intersections import gca_gca_intersection, gca_const_lat_intersection, get_number_of_intersections def test_get_GCA_GCA_intersections_antimeridian(): GCA1 = _lonlat_rad_to_xyz(np.deg2rad(170.0), np.deg2rad(89.99)) @@ -16,7 +16,7 @@ def test_get_GCA_GCA_intersections_antimeridian(): _lonlat_rad_to_xyz(np.deg2rad(70.0), 0.0), _lonlat_rad_to_xyz(np.deg2rad(179.0), 0.0) ]) - res_cart = _gca_gca_intersection_cartesian(GCR1_cart, GCR2_cart) + res_cart = gca_gca_intersection(GCR1_cart, GCR2_cart) assert len(res_cart) == 0 @@ -30,7 +30,7 @@ def test_get_GCA_GCA_intersections_antimeridian(): _lonlat_rad_to_xyz(np.deg2rad(175.0), 0.0) ]) - res_cart = _gca_gca_intersection_cartesian(GCR1_cart, GCR2_cart) + res_cart = gca_gca_intersection(GCR1_cart, GCR2_cart) res_cart = res_cart[0] assert np.allclose(np.linalg.norm(res_cart, axis=0), 1.0, atol=ERROR_TOLERANCE) @@ -47,7 +47,7 @@ def test_get_GCA_GCA_intersections_parallel(): _lonlat_rad_to_xyz(0.5 * np.pi, 0.0), _lonlat_rad_to_xyz(-0.5 * np.pi - 0.01, 0.0) ]) - res_cart = _gca_gca_intersection_cartesian(GCR1_cart, GCR2_cart) + res_cart = gca_gca_intersection(GCR1_cart, GCR2_cart) res_cart = res_cart[0] expected_res = np.array(_lonlat_rad_to_xyz(0.5 * np.pi, 0.0)) @@ -65,7 +65,7 @@ def test_get_GCA_GCA_intersections_perpendicular(): _lonlat_rad_to_xyz(*[0.5 * np.pi - 0.01, 0.0]), _lonlat_rad_to_xyz(*[-0.5 * np.pi + 0.01, 0.0]) ]) - res_cart = _gca_gca_intersection_cartesian(GCR1_cart, GCR2_cart) + res_cart = gca_gca_intersection(GCR1_cart, GCR2_cart) # rest_cart should be empty since these two GCAs are not intersecting assert(len(res_cart) == 0) @@ -254,3 +254,63 @@ def test_GCA_constLat_intersections_two_pts(): res = gca_const_lat_intersection(GCR1_cart, np.sin(query_lat)) assert res.shape[0] == 2 + +def test_on_minor_arc_near_antipodal_degeneracy(): + # Endpoints computed via independent (lon, lat) and (lon + pi, -lat) trig + # paths land antipodal to within 1-2 ulp, not bit-exact. A query point + # nowhere near the arc must not be classified as on it. + lon_a, lat_a = np.deg2rad(12.345), np.deg2rad(6.789) + a = np.array(_lonlat_rad_to_xyz(lon_a, lat_a)) + b = np.array(_lonlat_rad_to_xyz(lon_a + np.pi, -lat_a)) + assert not np.array_equal(b, -a) # confirms the degeneracy is not bit-exact + + q = np.array(_lonlat_rad_to_xyz(lon_a, lat_a + np.deg2rad(20))) + assert on_minor_arc(q, a, b) == 0 + +def test_on_minor_arc_near_coincident_degeneracy(): + # Same class of degeneracy as the antipodal case above, but for + # near-coincident (not bit-exact) endpoints. + lon_a, lat_a = np.deg2rad(12.345), np.deg2rad(6.789) + a = np.array(_lonlat_rad_to_xyz(lon_a, lat_a)) + b = np.array(_lonlat_rad_to_xyz(lon_a + 1e-16, lat_a)) + assert not np.array_equal(a, b) + + q = np.array(_lonlat_rad_to_xyz(lon_a + np.deg2rad(50), lat_a + np.deg2rad(30))) + assert on_minor_arc(q, a, b) == 0 + +def test_on_minor_arc_tiny_valid_edge_not_flagged_degenerate(): + # A genuinely short but non-degenerate edge (~1e-4 degrees, ~11m on Earth, + # well within realistic mesh resolution) must still validate its own + # midpoint and reject distant points -- guards against an overly-wide + # degeneracy tolerance in the |a x b|^2 check added for the antipodal / + # coincident fix. Note: on_minor_arc's absolute orient tolerance (1e-8, + # matching the AccuSphGeom C++ constant) stops discriminating for arcs + # shorter than roughly 1e-7 degrees (~1cm) -- that is a pre-existing + # limitation of the fixed-tolerance predicate shared with the C++ + # reference, not something this test targets. + lon_a, lat_a = np.deg2rad(10.0), np.deg2rad(10.0) + a = np.array(_lonlat_rad_to_xyz(lon_a, lat_a)) + b = np.array(_lonlat_rad_to_xyz(lon_a + 1e-4 * np.pi / 180.0, lat_a)) + + midpoint = (a + b) / 2.0 + midpoint = midpoint / np.linalg.norm(midpoint) + far_point = np.array(_lonlat_rad_to_xyz(np.deg2rad(-50), np.deg2rad(30))) + + assert on_minor_arc(midpoint, a, b) == 1 + assert on_minor_arc(far_point, a, b) == 0 + +def test_GCA_GCA_intersection_near_antipodal_arc_no_false_positive(): + # Regression test: gca_gca_intersection against an arc with near-antipodal + # (not bit-exact) endpoints must not report an intersection with an + # unrelated, geometrically distant GCA. + lon_a, lat_a = np.deg2rad(12.345), np.deg2rad(6.789) + a = np.array(_lonlat_rad_to_xyz(lon_a, lat_a)) + b = np.array(_lonlat_rad_to_xyz(lon_a + np.pi, -lat_a)) + gca_near_antipodal = np.array([a, b]) + + c0 = np.array(_lonlat_rad_to_xyz(np.deg2rad(-100.0), np.deg2rad(-40.0))) + c1 = np.array(_lonlat_rad_to_xyz(np.deg2rad(-99.0), np.deg2rad(-40.0))) + gca_unrelated = np.array([c0, c1]) + + res = gca_gca_intersection(gca_near_antipodal, gca_unrelated) + assert len(res) == 0 diff --git a/uxarray/grid/arcs.py b/uxarray/grid/arcs.py index 17426a196..f99941ef6 100644 --- a/uxarray/grid/arcs.py +++ b/uxarray/grid/arcs.py @@ -8,6 +8,17 @@ _normalize_xyz_scalar, ) from uxarray.grid.utils import _angle_of_2_vectors +from uxarray.utils.computing import accucross, two_sum + +# Magnitude below which orient3d_on_sphere classifies a result as zero. For +# double-precision unit-vector inputs this covers rounding error in the +# compensated cross product. +_PREDICATE_ZERO_TOL = 1e-15 + +# Tolerance for the on_minor_arc collinearity and interval tests. Matches +# AccuSphGeom's gca_*_minor_arc_tol (1e-8) so the Python predicate accepts the +# same candidate set as the C++ reference. +_ON_MINOR_ARC_TOL = 1e-8 def _to_list(obj): @@ -364,3 +375,153 @@ def compute_arc_length(pt_a, pt_b): delta_theta = np.arctan2(cross_2d, dot_2d) return rho * abs(delta_theta) + + +@njit(cache=True, inline="always") +def _normal_dot_value(nx_hi, ny_hi, nz_hi, nx_lo, ny_lo, nz_lo, q0, q1, q2): + """Compensated dot of an already-computed ``(hi, lo)`` normal with ``q``. + + This is the second half of :func:`_orient3d_on_sphere_value_xyz`, split out + so that callers holding a normal that is invariant across many queries (e.g. + the ray plane ``q x R`` in the spherical point-in-polygon kernel, which is + the same for every edge of a face) can compute the cross product once and + reuse it, paying only this dot per query. + """ + p0 = (nx_hi + nx_lo) * q0 + p1 = (ny_hi + ny_lo) * q1 + p2 = (nz_hi + nz_lo) * q2 + s, e = two_sum(p0, p1) + s, e2 = two_sum(s, p2) + return s + (e + e2) + + +@njit(cache=True, inline="always") +def _orient3d_on_sphere_value_xyz(a0, a1, a2, b0, b1, b2, q0, q1, q2): + """Accurate value of the orient3d-on-sphere predicate: ``(a x b) . q``. + + Takes the nine vector components directly so hot loops can call it without + materializing ``(3,)`` arrays. Uses a compensated cross product, so for + unit-vector inputs it provides roughly double the effective precision of a + naive evaluation. Positive when q is left of the directed arc a->b, negative + when right, near zero when q is on the great circle through a and b. + """ + nx_hi, ny_hi, nz_hi, nx_lo, ny_lo, nz_lo = accucross(a0, a1, a2, b0, b1, b2) + return _normal_dot_value(nx_hi, ny_hi, nz_hi, nx_lo, ny_lo, nz_lo, q0, q1, q2) + + +@njit(cache=True) +def _orient3d_on_sphere_value(a, b, q): + """Array form of :func:`_orient3d_on_sphere_value_xyz`: value of ``(a x b) . q``. + + Convenience wrapper taking three ``(3,)`` unit vectors. Positive when q is + left of the directed arc a->b, negative when right, near zero when q is on + the great circle through a and b. + """ + return _orient3d_on_sphere_value_xyz( + a[0], a[1], a[2], b[0], b[1], b[2], q[0], q[1], q[2] + ) + + +@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 compensated arithmetic 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. This is a public spatial predicate (used by + point-in-face, bounds, antimeridian handling and available for custom + geometry code); it is not part of the AccuXGCA/AccuXConstLat batch kernels. + + 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_xyz`` (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 + ------- + int + 1 if q lies on the minor arc ab, 0 otherwise. Returned as an integer + mask (not bool) so callers can multiply it into branch-free validity + products, mirroring AccuSphGeom's ``on_minor_arc_tol_ptr``. + """ + return _on_minor_arc_xyz(q[0], q[1], q[2], a[0], a[1], a[2], b[0], b[1], b[2], tol) + + +@njit(cache=True, inline="always") +def _on_minor_arc_xyz(q0, q1, q2, a0, a1, a2, b0, b1, b2, tol=_ON_MINOR_ARC_TOL): + """Scalar-argument form of :func:`on_minor_arc`, returning a 0/1 int mask. + + Same logic, but takes the nine vector components directly so hot loops can + test arc membership without allocating ``(3,)`` arrays for the query point. + """ + # Branch-free mask form, mirroring AccuSphGeom on_minor_arc_tol_ptr: the + # result is a product of 0/1 masks so the hot path has no data-dependent + # branches. + # + # Coincident/antipodal degeneracy (a == b or a == -b) is detected via + # |a x b|^2 rather than exact component equality. Endpoints computed + # through independent trig paths (e.g. one via (lon, lat), the other via + # (lon + pi, -lat)) land coincident/antipodal to 1-2 ulp, not bit-exact -- + # an exact ``==`` check misses them, and when missed, a x b ~= 0 makes the + # collinearity test pass for every point on the great circle and the + # interval checks degenerate to 0 >= -tol, producing false positives for + # arbitrary query points. |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) -- unlike + # comparing dot(a, b) to +-1, which suffers the same cancellation this PR + # exists to avoid and cannot make that distinction. The C++ reference's + # ``on_minor_arc_tol_ptr`` only guards exact coincidence and assumes + # non-antipodal mesh edges, so this widening is a UXarray-side addition, + # not a port of the reference. + cx = a1 * b2 - a2 * b1 + cy = a2 * b0 - a0 * b2 + cz = a0 * b1 - a1 * b0 + cross_sq = cx * cx + cy * cy + cz * cz + degenerate = 1 if cross_sq < 1e-30 else 0 + orient_ok = ( + 1 + if abs(_orient3d_on_sphere_value_xyz(a0, a1, a2, b0, b1, b2, q0, q1, q2)) <= tol + else 0 + ) + qa = a0 * q0 + a1 * q1 + a2 * q2 + qb = b0 * q0 + b1 * q1 + b2 * q2 + ab = a0 * b0 + a1 * b1 + a2 * b2 + s1_ok = 1 if (qb - ab * qa) >= -tol else 0 + s2_ok = 1 if (qa - qb * ab) >= -tol else 0 + return (1 - degenerate) * orient_ok * s1_ok * s2_ok diff --git a/uxarray/grid/intersections.py b/uxarray/grid/intersections.py index 97777de69..dfd408c51 100644 --- a/uxarray/grid/intersections.py +++ b/uxarray/grid/intersections.py @@ -1,16 +1,25 @@ +import math + import numpy as np from numba import njit, prange -from uxarray.constants import ERROR_TOLERANCE, INT_DTYPE, MACHINE_EPSILON -from uxarray.grid.arcs import ( - extreme_gca_z, - in_between, - point_within_gca, -) -from uxarray.grid.utils import ( - _angle_of_2_vectors, +from uxarray.constants import ERROR_TOLERANCE, INT_DTYPE +from uxarray.grid.arcs import _on_minor_arc_xyz, on_minor_arc +from uxarray.utils.computing import ( + _cdp2, + _cdp4, + _sum_of_squares_c, + acc_sqrt_re, + accucross, + accucross_pair, + two_prod, + two_sum, ) +# Edge screeners: fast O(n) passes used by Grid.get_edges_at_constant_* to +# identify candidate edges before the expensive GCA intersection. "no_extreme" +# means arc z-extrema along the great circle are not considered. + @njit(parallel=True, nogil=True, cache=True) def constant_lat_intersections_no_extreme(lat, edge_node_z, n_edge): @@ -292,155 +301,432 @@ def faces_within_lat_bounds(lats, face_bounds_lat): return candidate_faces -def _gca_gca_intersection_cartesian(gca_a_xyz, gca_b_xyz): - gca_a_xyz = np.asarray(gca_a_xyz) - gca_b_xyz = np.asarray(gca_b_xyz) +@njit(cache=True, inline="always", error_model="numpy") +def _accux_gca(w0, w1, v0, v1): + """Compute the candidate intersection points of two great-circle arcs. - return gca_gca_intersection(gca_a_xyz, gca_b_xyz) + Pure numerical kernel (mirrors AccuSphGeom ``accux_gca``). + Computes the two antipodal candidate intersection points of the great-circle + arcs w0-w1 and v0-v1. No branching, no validity filtering. -@njit(cache=True) + Parameters + ---------- + w0, w1 : np.ndarray, shape (3,) + Cartesian endpoints of the first arc. + v0, v1 : np.ndarray, shape (3,) + Cartesian endpoints of the second arc. + + Returns + ------- + pos, neg : np.ndarray, shape (3,) + Two antipodal candidate unit vectors. + """ + n1x_hi, n1y_hi, n1z_hi, n1x_lo, n1y_lo, n1z_lo = accucross( + w0[0], w0[1], w0[2], w1[0], w1[1], w1[2] + ) + n2x_hi, n2y_hi, n2z_hi, n2x_lo, n2y_lo, n2z_lo = accucross( + v0[0], v0[1], v0[2], v1[0], v1[1], v1[2] + ) + vx_hi, vy_hi, vz_hi, vx_lo, vy_lo, vz_lo = accucross_pair( + n1x_hi, + n1y_hi, + n1z_hi, + n1x_lo, + n1y_lo, + n1z_lo, + n2x_hi, + n2y_hi, + n2z_hi, + n2x_lo, + n2y_lo, + n2z_lo, + ) + vx = vx_hi + vx_lo + vy = vy_hi + vy_lo + vz = vz_hi + vz_lo + # Compensated norm: sum_of_squares_c over the (hi, lo) vector, then acc_sqrt_re + # folding the low part into the root, matching AccuSphGeom accux_gca. n = root.hi. + sum_hi, sum_lo = _sum_of_squares_c((vx_hi, vy_hi, vz_hi), (vx_lo, vy_lo, vz_lo)) + vn, _ = acc_sqrt_re(sum_hi, sum_lo) + # vn==0 (coplanar arcs) yields inf via IEEE division under error_model="numpy", + # so pos/neg become non-finite and the status layer masks them out. Branch-free. + inv = 1.0 / vn + pos = np.empty(3) + pos[0] = vx * inv + pos[1] = vy * inv + pos[2] = vz * inv + neg = np.empty(3) + neg[0] = -pos[0] + neg[1] = -pos[1] + neg[2] = -pos[2] + return pos, neg + + +@njit(cache=True, error_model="numpy") +def _try_gca_gca_intersection(w0, w1, v0, v1): + """Select the valid great-circle intersection and report a status code. + + Batch/status layer (mirrors AccuSphGeom ``try_gca_gca_intersection``). + + Calls the pure numerical kernel, applies integer mask arithmetic to determine + validity, selects the output point without if/else branching in the hot path. + + Status codes mirror AccuSphGeom: + 0 exactly one candidate is valid + 1 both candidates are valid + 2 neither candidate is valid (includes coplanar/parallel case) + """ + pos, neg = _accux_gca(w0, w1, v0, v1) + + pos_fin = ( + int(math.isfinite(pos[0])) + * int(math.isfinite(pos[1])) + * int(math.isfinite(pos[2])) + ) + neg_fin = ( + int(math.isfinite(neg[0])) + * int(math.isfinite(neg[1])) + * int(math.isfinite(neg[2])) + ) + pos_on_a = pos_fin * on_minor_arc(pos, w0, w1) + pos_on_b = pos_fin * on_minor_arc(pos, v0, v1) + neg_on_a = neg_fin * on_minor_arc(neg, w0, w1) + neg_on_b = neg_fin * on_minor_arc(neg, v0, v1) + + pos_valid = pos_fin * pos_on_a * pos_on_b + neg_valid = neg_fin * neg_on_a * neg_on_b + + pos_mask = pos_valid * (1 - neg_valid) + neg_mask = neg_valid * (1 - pos_valid) + + point = np.empty(3) + point[0] = pos_mask * pos[0] + neg_mask * neg[0] + point[1] = pos_mask * pos[1] + neg_mask * neg[1] + point[2] = pos_mask * pos[2] + neg_mask * neg[2] + + both = pos_valid * neg_valid + none = (1 - pos_valid) * (1 - neg_valid) + status = both + none * 2 + return point, status, pos, neg + + +@njit(cache=True, error_model="numpy") def gca_gca_intersection(gca_a_xyz, gca_b_xyz): - if gca_a_xyz.shape[1] != 3 or gca_b_xyz.shape[1] != 3: - raise ValueError("The two GCAs must be in the cartesian [x, y, z] format") + """Return the intersection points of two great-circle arcs. - # Extract points - w0_xyz = gca_a_xyz[0] - w1_xyz = gca_a_xyz[1] - v0_xyz = gca_b_xyz[0] - v1_xyz = gca_b_xyz[1] + Dispatcher / convenience API. Calls the batch/status layer and packages + results into UXarray's existing array-returning API. Coplanar/shared-endpoint + handling lives here, outside the numerical core. - angle_w0w1 = _angle_of_2_vectors(w0_xyz, w1_xyz) - angle_v0v1 = _angle_of_2_vectors(v0_xyz, v1_xyz) + Parameters + ---------- + gca_a_xyz : numpy.ndarray + First great-circle arc as two Cartesian endpoints, shape ``(2, 3)``. + gca_b_xyz : numpy.ndarray + Second great-circle arc as two Cartesian endpoints, shape ``(2, 3)``. - if angle_w0w1 > np.pi: - w0_xyz, w1_xyz = w1_xyz, w0_xyz + Returns + ------- + numpy.ndarray + Intersection points, shape ``(2, 3)``, with unused rows filled with NaN + (0, 1, or 2 valid rows). + """ + if gca_a_xyz.shape[1] != 3 or gca_b_xyz.shape[1] != 3: + raise ValueError("The two GCAs must be in the cartesian [x, y, z] format") - if angle_v0v1 > np.pi: - v0_xyz, v1_xyz = v1_xyz, v0_xyz + w0 = gca_a_xyz[0] + w1 = gca_a_xyz[1] + v0 = gca_b_xyz[0] + v1 = gca_b_xyz[1] - w0w1_norm = np.cross(w0_xyz, w1_xyz) - v0v1_norm = np.cross(v0_xyz, v1_xyz) - cross_norms = np.cross(w0w1_norm, v0v1_norm) + point, status, pos, neg = _try_gca_gca_intersection(w0, w1, v0, v1) - # Initialize result array and counter res = np.empty((2, 3)) count = 0 - - # Check if the two GCAs are parallel - if np.allclose(cross_norms, 0.0, atol=MACHINE_EPSILON): - if point_within_gca(v0_xyz, w0_xyz, w1_xyz): - res[count, :] = v0_xyz + if status == 0: + res[0, 0] = point[0] + res[0, 1] = point[1] + res[0, 2] = point[2] + count = 1 + elif status == 1: + res[0, 0] = pos[0] + res[0, 1] = pos[1] + res[0, 2] = pos[2] + res[1, 0] = neg[0] + res[1, 1] = neg[1] + res[1, 2] = neg[2] + count = 2 + else: + # status == 2: no candidate on both arcs. + # Check for coplanar overlap (shared endpoints) outside the kernel. + if on_minor_arc(v0, w0, w1): + res[count, 0] = v0[0] + res[count, 1] = v0[1] + res[count, 2] = v0[2] count += 1 - - if point_within_gca(v1_xyz, w0_xyz, w1_xyz): - res[count, :] = v1_xyz + if on_minor_arc(v1, w0, w1): + res[count, 0] = v1[0] + res[count, 1] = v1[1] + res[count, 2] = v1[2] count += 1 + return res[:count] - return res[:count, :] - # Normalize the cross_norms - cross_norms = cross_norms / np.linalg.norm(cross_norms) - x1_xyz = cross_norms - x2_xyz = -x1_xyz +@njit(cache=True, inline="always", error_model="numpy") +def _accux_constlat_scalar(a0, a1, a2, b0, b1, b2, const_z): + """Compute the constant-latitude intersection candidates, scalar in/out. - # Check intersection points - if point_within_gca(x1_xyz, w0_xyz, w1_xyz) and point_within_gca( - x1_xyz, v0_xyz, v1_xyz - ): - res[count, :] = x1_xyz - count += 1 + Allocation-free numerical kernel. Same compensated AccuSphGeom sequence as + :func:`_accux_constlat`, but takes the two arc endpoints as six scalars and + returns the two candidate points as six scalars (``pos`` xy and ``neg`` xy; + the z of both candidates is ``const_z``). Returning scalars instead of + ``np.empty(3)`` arrays lets Numba + keep everything in registers, so a batch loop over many edges does no + per-point heap allocation. This is the preferred entry point for hot loops. - if point_within_gca(x2_xyz, w0_xyz, w1_xyz) and point_within_gca( - x2_xyz, v0_xyz, v1_xyz - ): - res[count, :] = x2_xyz - count += 1 + Returns + ------- + px, py, nx_out, ny_out : float + ``pos = (px, py, const_z)`` and ``neg = (nx_out, ny_out, const_z)``. + Invalid inputs propagate as non-finite coordinates. + """ + nx_hi, ny_hi, nz_hi, nx_lo, ny_lo, nz_lo = accucross(a0, a1, a2, b0, b1, b2) + s2_hi, s2_lo = _sum_of_squares_c((nx_hi, ny_hi), (nx_lo, ny_lo)) + denom = s2_hi + s2_lo + s3_hi, s3_lo = _sum_of_squares_c((nx_hi, ny_hi, nz_hi), (nx_lo, ny_lo, nz_lo)) + zsq_hi, zsq_lo = two_prod(const_z, const_z) + d_hi, d_lo = _cdp4(s3_hi, zsq_hi, s3_hi, zsq_lo, s3_lo, zsq_hi, s3_lo, zsq_lo) + e_hi, e_lo = two_sum(s2_hi, -d_hi) + planar_sq = e_hi + (e_lo + s2_lo - d_lo) + s_root, s_corr = acc_sqrt_re(planar_sq) + nx = nx_hi + nx_lo + ny = ny_hi + ny_lo + nz = nz_hi + nz_lo + planar = s_root + s_corr + xp_hi, xp_lo = _cdp2(nx * nz, const_z, -ny, planar) + yp_hi, yp_lo = _cdp2(ny * nz, const_z, nx, planar) + xn_hi, xn_lo = _cdp2(nx * nz, const_z, ny, planar) + yn_hi, yn_lo = _cdp2(ny * nz, const_z, -nx, planar) + # denom == 0 (vertical arc) yields inf via IEEE division under + # error_model="numpy", so the isfinite mask in the status layer rejects the + # candidates. Branch-free. + inv_denom = 1.0 / denom + px = -(xp_hi + xp_lo) * inv_denom + py = -(yp_hi + yp_lo) * inv_denom + nxo = -(xn_hi + xn_lo) * inv_denom + nyo = -(yn_hi + yn_lo) * inv_denom + return px, py, nxo, nyo + + +@njit(cache=True, inline="always") +def _accux_constlat(x1, x2, const_z): + """Compute the two constant-latitude intersection candidates as arrays. + + Pure numerical kernel (mirrors AccuSphGeom ``accux_constlat``). An + array-returning wrapper around :func:`_accux_constlat_scalar`. Computes the + two candidate intersection points between the great-circle arc defined by + unit vectors *x1*, *x2* and the constant-latitude plane z = const_z. No + branching, no validity filtering. For allocation-free hot loops call + :func:`_accux_constlat_scalar` directly. - return res[:count, :] + Parameters + ---------- + x1, x2 : np.ndarray, shape (3,) + Cartesian endpoints of the great-circle arc. + const_z : float + Constant-latitude plane, given as the Cartesian z value ``sin(lat)``. + + Returns + ------- + pos, neg : np.ndarray, shape (3,) + Two antipodal candidate points. Invalid inputs propagate as non-finite + coordinates; the caller uses masks/status to identify validity. + """ + px, py, nxo, nyo = _accux_constlat_scalar( + x1[0], x1[1], x1[2], x2[0], x2[1], x2[2], const_z + ) + pos = np.empty(3) + pos[0] = px + pos[1] = py + pos[2] = const_z + neg = np.empty(3) + neg[0] = nxo + neg[1] = nyo + neg[2] = const_z + return pos, neg + + +@njit(cache=True, error_model="numpy") +def _try_gca_const_lat_intersection(gca_cart, const_z): + """Select the valid constant-latitude intersection and report a status code. + + Batch/status layer (mirrors AccuSphGeom ``try_gca_constlat_intersection``). + + Calls the pure numerical kernel, computes integer validity masks (0 or 1) + for each candidate using finiteness and arc-membership tests, then selects + the output point via integer arithmetic — no if/else branching in the hot path. + + Status codes mirror AccuSphGeom: + 0 exactly one candidate is valid (normal case) + 1 both candidates are valid + 2 neither candidate is valid + """ + x1 = gca_cart[0] + x2 = gca_cart[1] + pos, neg = _accux_constlat(x1, x2, const_z) + + pos_fin = int(math.isfinite(pos[0])) * int(math.isfinite(pos[1])) + neg_fin = int(math.isfinite(neg[0])) * int(math.isfinite(neg[1])) + pos_on = pos_fin * on_minor_arc(pos, x1, x2) + neg_on = neg_fin * on_minor_arc(neg, x1, x2) + + pos_valid = pos_fin * pos_on + neg_valid = neg_fin * neg_on + + pos_mask = pos_valid * (1 - neg_valid) + neg_mask = neg_valid * (1 - pos_valid) + + point = np.empty(3) + point[0] = pos_mask * pos[0] + neg_mask * neg[0] + point[1] = pos_mask * pos[1] + neg_mask * neg[1] + point[2] = pos_mask * pos[2] + neg_mask * neg[2] + + both = pos_valid * neg_valid + none = (1 - pos_valid) * (1 - neg_valid) + status = both + none * 2 + return point, status, pos, neg @njit(cache=True) +def _snap_const_lat_endpoint(point, x1, x2, const_z): + """Snap a candidate point to an arc endpoint when the endpoint lies on the latitude.""" + # 1e-14 is distance² in Cartesian between candidate and endpoint; corresponds + # to ~1e-7 in arc length (unit sphere). Candidates within this distance are + # snapped to the exact endpoint to avoid sub-ulp drift when the arc ends + # exactly on the latitude circle. + sx, sy = _snap_const_lat_endpoint_xy( + point[0], point[1], x1[0], x1[1], x1[2], x2[0], x2[1], x2[2], const_z + ) + out = np.empty(3) + out[0] = sx + out[1] = sy + out[2] = point[2] + return out + + +@njit(cache=True, inline="always") +def _snap_const_lat_endpoint_xy(px, py, a0, a1, a2, b0, b1, b2, const_z): + """Scalar-argument form of :func:`_snap_const_lat_endpoint`. + + Returns the (possibly snapped) x, y of the candidate; z is always ``const_z`` + so it is not returned. Allocation-free for use in hot loops. + """ + snap_sq = 1e-14 + ox = px + oy = py + if abs(a2 - const_z) <= ERROR_TOLERANCE: + dx = ox - a0 + dy = oy - a1 + if dx * dx + dy * dy < snap_sq: + ox = a0 + oy = a1 + if abs(b2 - const_z) <= ERROR_TOLERANCE: + dx = ox - b0 + dy = oy - b1 + if dx * dx + dy * dy < snap_sq: + ox = b0 + oy = b1 + return ox, oy + + +@njit(cache=True, error_model="numpy") def gca_const_lat_intersection(gca_cart, const_z): - """Calculate the intersection point(s) of a Great Circle Arc (GCA) and a - constant latitude line in a Cartesian coordinate system. + """Return the intersection points of a great-circle arc and a latitude. + + Dispatcher / convenience API. Runs the numerical kernel, validity masks, + endpoint snapping, and packaging into UXarray's NaN-filled (2, 3) format + entirely on scalars, so the only heap allocation is the returned array. All + UXarray-specific branching lives here so the numerical core stays uniform. + See ``_try_gca_const_lat_intersection`` for the array-returning form used by + the layer benchmarks. Parameters ---------- - gca_cart : [2, 3] np.ndarray Cartesian coordinates of the two end points GCA. + gca_cart : numpy.ndarray + Great-circle arc as two Cartesian endpoints, shape ``(2, 3)``. const_z : float - The constant latitude represented in cartesian of the latitude line. + Constant-latitude plane, given as the Cartesian z value ``sin(lat)``. Returns ------- - np.ndarray - Cartesian coordinates of the intersection point(s) the shape is [2, 3]. If no intersections are found, - all values a `nan`. If one intersection is found, the first column represent the intersection point, and - if two intersections are found, each column represents a point. - + numpy.ndarray + Intersection points, shape ``(2, 3)``, with unused rows filled with NaN + (0, 1, or 2 valid rows). """ res = np.empty((2, 3)) res.fill(np.nan) - x1, x2 = gca_cart - - # Check if the constant latitude has the same latitude as the GCA endpoints - x1_at_const_z = np.isclose( - x1[2], const_z, rtol=ERROR_TOLERANCE, atol=ERROR_TOLERANCE - ) - x2_at_const_z = np.isclose( - x2[2], const_z, rtol=ERROR_TOLERANCE, atol=ERROR_TOLERANCE - ) - - if x1_at_const_z and x2_at_const_z: - res[0] = x1 - res[1] = x2 - return res - elif x1_at_const_z: - res[0] = x1 - return res - elif x2_at_const_z: - res[0] = x2 - return res - - # If the constant latitude is not the same as the GCA endpoints, calculate the intersection point - z_min = extreme_gca_z(gca_cart, extreme_type="min") - 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): - return res - - n = np.cross(x1, x2) - - nx, ny, nz = n - - s_tilde = np.sqrt(nx**2 + ny**2 - (nx**2 + ny**2 + nz**2) * const_z**2) - p1_x = -(1.0 / (nx**2 + ny**2)) * (const_z * nx * nz + s_tilde * ny) - p2_x = -(1.0 / (nx**2 + ny**2)) * (const_z * nx * nz - s_tilde * ny) - p1_y = -(1.0 / (nx**2 + ny**2)) * (const_z * ny * nz - s_tilde * nx) - p2_y = -(1.0 / (nx**2 + ny**2)) * (const_z * ny * nz + s_tilde * nx) - - p1 = np.array([p1_x, p1_y, const_z]) - p2 = np.array([p2_x, p2_y, const_z]) - - p1_intersects_gca = point_within_gca(p1, gca_cart[0], gca_cart[1]) - p2_intersects_gca = point_within_gca(p2, gca_cart[0], gca_cart[1]) - - if p1_intersects_gca and p2_intersects_gca: - res[0] = p1 - res[1] = p2 - elif p1_intersects_gca: - res[0] = p1 - elif p2_intersects_gca: - res[0] = p2 - + a0 = gca_cart[0, 0] + a1 = gca_cart[0, 1] + a2 = gca_cart[0, 2] + b0 = gca_cart[1, 0] + b1 = gca_cart[1, 1] + b2 = gca_cart[1, 2] + + px, py, nx, ny = _accux_constlat_scalar(a0, a1, a2, b0, b1, b2, const_z) + + pos_fin = int(math.isfinite(px)) * int(math.isfinite(py)) + neg_fin = int(math.isfinite(nx)) * int(math.isfinite(ny)) + pos_valid = pos_fin * _on_minor_arc_xyz(px, py, const_z, a0, a1, a2, b0, b1, b2) + neg_valid = neg_fin * _on_minor_arc_xyz(nx, ny, const_z, a0, a1, a2, b0, b1, b2) + + if pos_valid ^ neg_valid: + if pos_valid: + sx, sy = _snap_const_lat_endpoint_xy( + px, py, a0, a1, a2, b0, b1, b2, const_z + ) + else: + sx, sy = _snap_const_lat_endpoint_xy( + nx, ny, a0, a1, a2, b0, b1, b2, const_z + ) + res[0, 0] = sx + res[0, 1] = sy + res[0, 2] = const_z + elif pos_valid and neg_valid: + psx, psy = _snap_const_lat_endpoint_xy(px, py, a0, a1, a2, b0, b1, b2, const_z) + nsx, nsy = _snap_const_lat_endpoint_xy(nx, ny, a0, a1, a2, b0, b1, b2, const_z) + dx = psx - nsx + dy = psy - nsy + if dx * dx + dy * dy < 1e-14: + res[0, 0] = psx + res[0, 1] = psy + res[0, 2] = const_z + else: + res[0, 0] = psx + res[0, 1] = psy + res[0, 2] = const_z + res[1, 0] = nsx + res[1, 1] = nsy + res[1, 2] = const_z return res @njit(cache=True) def get_number_of_intersections(arr): - """Returns the number of intersection points for the output of the gca-const-lat intersection.""" + """Return the number of intersection points in a gca-const-lat result. + + Parameters + ---------- + arr : numpy.ndarray + Output of :func:`gca_const_lat_intersection`, shape ``(2, 3)`` with + unused rows filled with NaN. + + Returns + ------- + int + Number of non-NaN intersection points (0, 1, or 2). + """ row1_is_nan = np.all(np.isnan(arr[0])) row2_is_nan = np.all(np.isnan(arr[1])) diff --git a/uxarray/utils/computing.py b/uxarray/utils/computing.py index ca5ca5183..80b661cad 100644 --- a/uxarray/utils/computing.py +++ b/uxarray/utils/computing.py @@ -1,628 +1,521 @@ -import sys +"""Compensated floating-point primitives for accurate spherical geometry. -import numpy as np +Cross and dot products over nearly-parallel unit vectors suffer catastrophic +cancellation in double precision, which degrades GCA-GCA and constant-latitude +intersections and the point-in-polygon ray test. These ``@njit`` primitives +recover near-double precision using error-free transformations (``two_sum``, +``two_prod``) and compensated algorithms built on them. ``two_prod`` uses a +hardware FMA when one is available (validated bit-exact at import time) and +falls back to the Veltkamp split otherwise, so there is no hard FMA dependency. +Python/Numba port of the AccuSphGeom C++ library (EFT tier only; the adaptive +Shewchuk predicate and exact-arithmetic fallback tiers are not ported): -def _fmms(a, b, c, d): - """ - Calculate the difference of products using the FMA (fused multiply-add) operation: (a * b) - (c * d). - - This operation leverages the fused multiply-add operation when available on the system and rounds the result only once. - The relative error of this operation is bounded by 1.5 ulps when no overflow and underflow occur. - - Parameters - ---------- - a (float): The first value of the first product. - b (float): The second value of the first product. - c (float): The first value of the second product. - d (float): The second value of the second product. - - Returns - ------- - float: The difference of the two products. - - Example - ------- - >>> _fmms(3.0, 2.0, 1.0, 1.0) - 5.0 - - Reference - --------- - Claude-Pierre Jeannerod, Nicolas Louvet, and Jean-Michel Muller, Further - analysis of Kahan’s algorithm for the accurate computation of 2 x 2 determinants, - Mathematics of Computation, vol. 82, no. 284, pp. 2245-2264, 2013. - [Read more](https://ens-lyon.hal.science/ensl-00649347) (DOI: 10.1090/S0025-5718-2013-02679-8) - """ - import pyfma - - cd = c * d - err = pyfma.fma(-c, d, cd) - dop = pyfma.fma(a, b, -cd) - return dop + err + Chen, H. (2026). Accurate and Robust Algorithms for Spherical Polygon + Operations. EGUsphere preprint egusphere-2026-636. + Chen, H. Great Circle Arc Intersection and Constant Latitude Intersection + on the Sphere. SIAM J. Sci. Comput. https://doi.org/10.1137/25M1737614 + Reference implementation: https://github.com/hongyuchen1030/AccuSphGeom +""" +import math -def cross_fma(v1, v2): - """Calculate the cross product of two 3D vectors utilizing the fused - multiply-add operation. - - Parameters - ---------- - v1 (np.array): The first vector of size 3. - v2 (np.array): The second vector of size 3. - - Returns - ------- - np.array: The cross product vector of size 3. - - Example - ------- - >>> v1 = np.array([1.0, 2.0, 3.0]) - >>> v2 = np.array([4.0, 5.0, 6.0]) - >>> cross_fma(v1, v2) - array([-3.0, 6.0, -3.0]) - """ - x = _fmms(v1[1], v2[2], v1[2], v2[1]) - y = _fmms(v1[2], v2[0], v1[0], v2[2]) - z = _fmms(v1[0], v2[1], v1[1], v2[0]) - return np.array([x, y, z]) - - -def dot_fma(v1, v2): - """Calculate the dot product of two vectors using the FMA (fused multiply- - add) operation. - - This implementation leverages the FMA operation to provide a more accurate result. Currently the ComptDot product - algorithm is used, which provides a relative error of approvimately u + n^2u^2cond(v1 dot v2), where u is 0.5 ulps, - n is the length of the vectors, and cond(v1 dot v2) is the condition number of the naive dot product of v1 and v2. - This operatin takes approvimately 3 + 10 * n flops, where n is the length of the vectors. - - Parameters - ---------- - v1 : list of float - The first vector. - v2 : list of float - The second vector. Must be the same length as v1. - - Returns - ------- - float - The dot product of the two vectors. - - Raises - ------ - ValueError - If the input vectors `v1` and `v2` are not of the same length. - - Examples - -------- - >>> dot_fma([1.0, 2.0, 3.0], [4.0, 5.0, 6.0]) - 32.0 - - References - ---------- - S. Graillat, Ph. Langlois, and N. Louvet. "Accurate dot products with FMA." Presented at RNC 7, 2007, Nancy, France. - DALI-LP2A Laboratory, University of Perpignan, France. - """ - if len(v1) != len(v2): - raise ValueError("Input vectors must be of the same length") +from numba import njit - s, c = _two_prod_fma(v1[0], v2[0]) - for i in range(1, len(v1)): - p, pi = _two_prod_fma(v1[i], v2[i]) - s, signma = _two_sum(s, p) - c = c + pi + signma +# Fused multiply-add (FMA) support. ``two_prod`` needs the exact rounding error +# of ``a * b``; with hardware FMA this is ``e = fma(a, b, -p)``, otherwise we +# fall back to the portable Veltkamp split. The LLVM ``fma`` intrinsic is +# validated at import time; if it is unavailable or non-exact, ``_HAS_FMA`` +# stays False and the Veltkamp path is used. +try: + from numba.core import types as _nb_types + from numba.extending import intrinsic as _nb_intrinsic - return s + c + @_nb_intrinsic + def _fma(typingctx, a, b, c): + sig = _nb_types.float64(_nb_types.float64, _nb_types.float64, _nb_types.float64) + def codegen(context, builder, signature, args): + return builder.fma(*args) -def _two_prod_fma(a, b): - """Error-free transformation of the product of two floating-point numbers - using FMA, such that a * b = x + y exactly. - - Parameters - ---------- - a, b : float - The floating-point numbers to be multiplied. - - Returns - ------- - tuple of float - The product and the error term. - - Examples - -------- - >>> _two_prod_fma(1.0, 2.0) - (2.0, 0.0) - - Reference - --------- - Stef Graillat. Accurate Floating Point Product and Exponentiation. - IEEE Transactions on Computers, 58(7), 994–1000, 2009.10.1109/TC.2008.215. - """ - import pyfma - - x = a * b - y = pyfma.fma(a, b, -x) - return x, y - - -def _err_fmac(a, b, c): - """Error-free transformation for the FMA operation. such that x = - FMA(a,b,c) and a * b + c = x + y + z exactly. Thhis function is only - available in round to the nearest mode and takes approximately 17 flops. - - Parameters - ---------- - a, b, c : float - The operands for the FMA operation. - - Returns - ------- - tuple of float - The result of the FMA operation and two error terms. - - References - ---------- - Graillat, Stef & Langlois, Philippe & Louvet, Nicolas. (2006). Improving the compensated Horner scheme with - a Fused Multiply and Add. 2. 1323-1327. 10.1145/1141277.1141585. - - Ogita, Takeshi & Rump, Siegfried & Oishi, Shin’ichi. (2005). Accurate Sum and Dot Product. - SIAM J. Scientific Computing. 26. 1955-1988. 10.1137/030601818. - """ - if sys.float_info.rounds == 1: - import pyfma - - x = pyfma.fma(a, b, c) - u1, u2 = _fast_two_mult(a, b) - alpha1, alpha2 = _two_sum(c, u2) - beta1, beta2 = _two_sum(u1, alpha1) - gamma = (beta1 - x) + beta2 - y, z = _fast_two_sum(gamma, alpha2) - return x, y, z - else: - raise ValueError( - "3FMA operation is only available in round to the nearest mode. and the current mode is " - + str(sys.float_info.rounds) - ) - - -def _two_sum(a, b): - """Error-free transformation of the sum of two floating-point numbers such - that a + b = x + y exactly. - - Parameters - ---------- - a, b : float - The floating-point numbers to be added. - - Returns - ------- - tuple of float - The sum and the error term. - - Examples - -------- - >>> _two_sum(1.0, 2.0) - (3.0, 0.0) + return sig, codegen - Reference - --------- - D. Knuth. 1998. The Art of Computer Programming (3rd ed.). Vol. 2. Addison-Wesley, Reading, MA. - """ - x = a + b - z = x - a - y = (a - (x - z)) + (b - z) - return x, y + _FMA_INTRINSIC_OK = True +except Exception: # pragma: no cover - toolchain without intrinsic support + _FMA_INTRINSIC_OK = False -def _fast_two_mult(a, b): - """Error-free transformation of the product of two floating-point numbers - such that a * b = x + y exactly. +@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. - This function is faster than the _two_prod_fma function. + 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 - The floating-point numbers to be multiplied. + Input values. Returns ------- - tuple of float - The product and the error term. - - References - ---------- - Vincent Lefèvre, Nicolas Louvet, Jean-Michel Muller, Joris Picot, and Laurence Rideau. 2023. - Accurate Calculation of Euclidean Norms Using Double-word Arithmetic. - ACM Trans. Math. Softw. 49, 1, Article 1 (March 2023), 34 pages. https://doi.org/10.1145/3568672 + s : float + Rounded sum fl(a + b). + e : float + Rounding error term; s + e = a + b exactly. """ - import pyfma + s = a + b + bp = s - a + e = (a - (s - bp)) + (b - bp) + return s, e - x = a * b - y = pyfma.fma(a, b, -x) - return x, y +if _FMA_INTRINSIC_OK: -def _fast_two_sum(a, b): - """Compute a fast error-free transformation of the sum of two floating- - point numbers. + @njit(cache=True, inline="always") + def _two_prod_fma(a, b): + p = a * b + return p, _fma(a, b, -p) - This function is a faster alternative to `_two_sum` for computing the sum - of two floating-point numbers `a` and `b`, such that a + b = x + y exactly. - Note: |a| must be no less than |b|. - Parameters - ---------- - a, b : float - The floating-point numbers to be added. It is required that |a| >= |b|. +@njit(cache=True, inline="always") +def _two_prod_veltkamp(a, b): + """Portable TwoProd via Veltkamp splitting (no FMA dependency). - Returns - ------- - tuple of float - The rounded sum of `a` and `b`, and the error term. The error term represents the difference between the exact sum and the rounded sum. - - Raises - ------ - ValueError - If |a| < |b|. - - Examples - -------- - >>> _fast_two_sum(2.0, 1.0) - (3.0, 0.0) - - >>> _fast_two_sum(1.0, 2.0) - Traceback (most recent call last): - ... - ValueError: |a| must be greater than or equal to |b|. - - Reference - --------- - T. J. Dekker. A Floating-Point Technique for Extending the Available Precision. - Numerische Mathematik, 18(3), 224–242,1971. 10.1007/BF01397083. - Available at: https://doi.org/10.1007/BF01397083. + Decomposes each operand into a high and low half using the splitting + constant 2**27 + 1, then reconstructs the exact rounding error from the + four partial products. Works on every Numba target. """ - if abs(a) >= abs(b): - x = a + b - b_tile = x - a - y = b - b_tile - return x, y - - else: - raise ValueError("|a| must be greater than or equal to |b|.") - - -def _comp_prod_fma(vec): - """Compute the compensated product using Fused Multiply-Add (FMA). - - This function computes the product of elements in a vector using a - compensated algorithm with Fused Multiply-Add to reduce numerical errors. + 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 + + +def _validate_fma(n_samples=2000) -> bool: + """Return True iff the FMA intrinsic compiles and is a bit-exact EFT.""" + if not _FMA_INTRINSIC_OK: + return False + try: + import numpy as _np + + rng = _np.random.default_rng(20260101) + for _ in range(n_samples): + a = float(rng.standard_normal() * rng.integers(1, 1 << 20)) + b = float(rng.standard_normal() * rng.integers(1, 1 << 20)) + pf, ef = _two_prod_fma(a, b) + pv, ev = _two_prod_veltkamp(a, b) + # Compare the residuals *directly*. Do not compare ``pf + ef`` + # against ``pv + ev``: both sums round straight back to the product + # (|e| <= ulp(p)/2 by construction), so that predicate collapses to + # ``pf != pv`` -- a tautology, since both are fl(a*b). The exact + # residual is unique and representable, so a correct FMA and the + # Veltkamp split must agree bit-for-bit. + if pf != pv or ef != ev: + return False + return True + except Exception: # pragma: no cover + return False + + +_HAS_FMA = _validate_fma() + + +if _HAS_FMA: + + @njit(cache=True, inline="always") + def two_prod(a, b): + """Dekker TwoProd: return (p, e) with p = fl(a*b) and p + e = a*b exactly. + + Uses a single fused multiply-add for the error term on hardware that + supports it (selected at import time via ``_HAS_FMA``), falling back to + the portable Veltkamp split otherwise. The FMA path is ~2x faster in the + compensated geometry kernels and is bit-for-bit identical to the + Veltkamp result (validated at import). + + Parameters + ---------- + a, b : float + Input values. + + Returns + ------- + p : float + Rounded product fl(a * b). + e : float + Rounding error term; p + e = a * b exactly. + """ + return _two_prod_fma(a, b) + +else: # pragma: no cover - exercised only on FMA-less toolchains + + @njit(cache=True, inline="always") + def two_prod(a, b): + """Dekker TwoProd: return (p, e) with p = fl(a*b) and p + e = a*b exactly. + + Portable Veltkamp-split implementation (no FMA available on this + toolchain/target). + + Parameters + ---------- + a, b : float + Input values. + + Returns + ------- + p : float + Rounded product fl(a * b). + e : float + Rounding error term; p + e = a * b exactly. + """ + return _two_prod_veltkamp(a, b) + + +@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 ---------- - vec : list of float - The vector whose elements are to be multiplied. + a, b, c, d : float + Input scalars; computes a*b - c*d. Returns ------- - float - The compensated product of the elements in the vector. - - Examples - -------- - >>> _comp_prod_fma([1.1, 2.2, 3.3]) - 7.986000000000001 - - Reference - --------- - Takeshi Ogita, Siegfried M. Rump, and Shin'ichi Oishi. 2005. Accurate Sum and Dot Product. - SIAM J. Sci. Comput. 26, 6 (2005), 1955–1988. https://doi.org/10.1137/030601818 + hi : float + High-order part of the accurate result. + lo : float + Low-order correction term; hi + lo equals the accurate value. """ - import pyfma + 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_s - e_w) + return s, lo - p1 = vec[0] - e1 = 0.0 - for i in range(1, len(vec)): - p_i, pi = _two_prod_fma(p1, vec[i]) - ei = pyfma.fma(e1, vec[i], pi) - p1 = p_i - e1 = ei - res = p1 + e1 - return res +@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. -def _sum_of_squares_re(vec): - """Compute the sum of squares of a vector using a compensated algorithm. - - This function calculates the sum of squares of the elements in a vector, - employing a compensation technique to reduce numerical errors. + 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 ---------- - vec : list of float - The vector whose elements' squares are to be summed. + a0, a1, a2 : float + Components of vector a. + b0, b1, b2 : float + Components of vector b. Returns ------- - float - The compensated sum of the squares of the elements in the vector. - - Examples - -------- - >>> _sum_of_squares_re([1.0, 2.0, 3.0]) - 14.0 - - Reference - --------- - Stef Graillat, Christoph Lauter, PING Tak Peter Tang, - Naoya Yamanaka, and Shin’ichi Oishi. Efficient Calculations of Faith- - fully Rounded L2-Norms of n-Vectors. ACM Transactions on Mathemat- - ical Software, 41(4), Article 24, 2015. 10.1145/2699469. Available at: - https://doi.org/10.1145/2699469. + x_hi, y_hi, z_hi, x_lo, y_lo, z_lo : float + High and low parts of each cross-product component. """ - P, p = _two_square(vec) - S, s = _two_sum(P[0], P[1]) - for i in range(2, len(vec)): - H, h = _two_sum(S, P[i]) - S, s = _two_sum(H, s + h) - sump = sum(p) - H, h = _two_sum(S, sump) - S, s = _fast_two_sum(H, s + h) - return S + s - - -def _vec_sum(p): - """Compute the sum of a vector using a compensated summation algorithm. - - This function calculates the sum of the elements in a vector using a - compensated summation algorithm to reduce numerical errors. - - Parameters - ---------- - p : list of float - The vector whose elements are to be summed. - - Returns - ------- - float - The compensated sum of the elements in the vector. - - Examples - -------- - >>> _vec_sum([1.0, 2.0, 3.0]) - 6.0 - - Reference - --------- - Takeshi Ogita, Siegfried M. Rump, and Shin'ichi Oishi. 2005. Accurate Sum and Dot Product. - SIAM J. Sci. Comput. 26, 6 (2005), 1955–1988. https://doi.org/10.1137/030601818 + 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 + + +@njit(cache=True, inline="always") +def _cdp8( + a0, + a1, + a2, + a3, + a4, + a5, + a6, + a7, + b0, + b1, + b2, + b3, + b4, + b5, + b6, + b7, +): + """Compensated sum of 8 exact products: Σ ai*bi, i=0..7. + + Uses ``two_prod`` + ``two_sum`` accumulation (Ogita-Rump-Oishi style) + so the result has error bounded by one ulp of the true value regardless + of cancellation in intermediate sums. """ - pi_1 = p[0] - sigma_i1 = 0 - - for i in range(1, len(p)): - pi, qi = _two_sum(pi_1, p[i]) - sigma_i = sigma_i1 + qi - pi_1 = pi - sigma_i1 = sigma_i - - res = pi_1 + sigma_i1 - return res - - -def _norm_faithful(x): - """Compute the faithful norm of a vector. - - This function calculates the faithful norm (L2 norm) of a vector, - which is a more numerically stable version of the Euclidean norm. + s, lo = two_prod(a0, b0) + p, e = two_prod(a1, b1) + s2, e2 = two_sum(s, p) + lo += e + e2 + s = s2 + p, e = two_prod(a2, b2) + s2, e2 = two_sum(s, p) + lo += e + e2 + s = s2 + p, e = two_prod(a3, b3) + s2, e2 = two_sum(s, p) + lo += e + e2 + s = s2 + p, e = two_prod(a4, b4) + s2, e2 = two_sum(s, p) + lo += e + e2 + s = s2 + p, e = two_prod(a5, b5) + s2, e2 = two_sum(s, p) + lo += e + e2 + s = s2 + p, e = two_prod(a6, b6) + s2, e2 = two_sum(s, p) + lo += e + e2 + s = s2 + p, e = two_prod(a7, b7) + s2, e2 = two_sum(s, p) + lo += e + e2 + s = s2 + return s, lo + + +@njit(cache=True, inline="always") +def accucross_pair( + ax_hi, + ay_hi, + az_hi, + ax_lo, + ay_lo, + az_lo, + bx_hi, + by_hi, + bz_hi, + bx_lo, + by_lo, + bz_lo, +): + """Compensated cross product of two compensated vectors. + + Computes (a_hi + a_lo) × (b_hi + b_lo) using a compensated 8-term dot + product for each component, matching the two-argument ``accucross`` overload + in the AccuSphGeom C++ library. This is more accurate than collapsing + (hi, lo) to a single float before the cross product. Parameters ---------- - x : list of float - The vector whose norm is to be computed. + ax_hi, ay_hi, az_hi : float + High parts of vector a. + ax_lo, ay_lo, az_lo : float + Low parts of vector a (rounding residuals from a prior compensated operation). + bx_hi, by_hi, bz_hi : float + High parts of vector b. + bx_lo, by_lo, bz_lo : float + Low parts of vector b. Returns ------- - float - The faithful norm of the vector. - - Examples - -------- - >>> _norm_faithful([1.0, 2.0, 3.0]) - 3.7416573867739413 + x_hi, y_hi, z_hi, x_lo, y_lo, z_lo : float + Compensated cross-product components. """ - return _norm_l(x) - - -def _norm_l(x): - """Compute the L2 norm (Euclidean norm) of a vector using a compensated - algorithm. - - This function calculates the L2 norm of a vector, employing a compensation - technique to reduce numerical errors during the computation. It involves - computing the sum of squares of the vector elements in a numerically stable way. - - Parameters - ---------- - x : list of float - The vector whose L2 norm is to be computed. + # x = (ay*bz) - (az*by), expanded over all four hi/lo cross-terms + x_hi, x_lo = _cdp8( + ay_hi, + ay_hi, + ay_lo, + ay_lo, + -az_hi, + -az_hi, + -az_lo, + -az_lo, + bz_hi, + bz_lo, + bz_hi, + bz_lo, + by_hi, + by_lo, + by_hi, + by_lo, + ) + # y = (az*bx) - (ax*bz) + y_hi, y_lo = _cdp8( + az_hi, + az_hi, + az_lo, + az_lo, + -ax_hi, + -ax_hi, + -ax_lo, + -ax_lo, + bx_hi, + bx_lo, + bx_hi, + bx_lo, + bz_hi, + bz_lo, + bz_hi, + bz_lo, + ) + # z = (ax*by) - (ay*bx) + z_hi, z_lo = _cdp8( + ax_hi, + ax_hi, + ax_lo, + ax_lo, + -ay_hi, + -ay_hi, + -ay_lo, + -ay_lo, + by_hi, + by_lo, + by_hi, + by_lo, + bx_hi, + bx_lo, + bx_hi, + bx_lo, + ) + return x_hi, y_hi, z_hi, x_lo, y_lo, z_lo + + +# Compensated dot products and sum-of-squares, fixed small sizes (Numba does not +# support generic runtime-length accumulations inside @njit). + + +@njit(cache=True, inline="always") +def _cdp2(a0, b0, a1, b1): + """Compensated dot product of 2 pairs: a0*b0 + a1*b1.""" + s, lo = two_prod(a0, b0) + p, e = two_prod(a1, b1) + s2, e2 = two_sum(s, p) + lo += e + e2 + return s2, lo + + +@njit(cache=True, inline="always") +def _cdp4(a0, b0, a1, b1, a2, b2, a3, b3): + """Compensated dot product of 4 pairs: Σ ai*bi, i=0..3.""" + s, lo = two_prod(a0, b0) + p, e = two_prod(a1, b1) + s2, e2 = two_sum(s, p) + lo += e + e2 + s = s2 + p, e = two_prod(a2, b2) + s2, e2 = two_sum(s, p) + lo += e + e2 + s = s2 + p, e = two_prod(a3, b3) + s2, e2 = two_sum(s, p) + lo += e + e2 + return s2, lo + + +@njit(cache=True, inline="always") +def _fast_two_sum(a, b): + """Fast TwoSum: (x, y) with x = fl(a + b), x + y = a + b exactly. - Returns - ------- - float - The compensated L2 norm of the vector. - - Examples - -------- - >>> _norm_l([1.0, 2.0, 3.0]) - 3.7416573867739413 - - Reference - --------- - Vincent Lef`evre, Nicolas Louvet, Jean-Michel Muller, - Joris Picot, and Laurence Rideau. Accurate Calculation of Euclidean - Norms Using Double-Word Arithmetic. ACM Transactions on Mathemat- - ical Software, 49(1), 1–34, March 2023. 10.1145/3568672 + Requires ``|a| >= |b|`` (or ``a == 0``) for the error term to be exact — + the callers here satisfy this because ``a`` is a running non-negative sum. """ - P, p = _two_square(x) - S, s = _two_sum(P[0], P[1]) - for i in range(2, len(x)): - H, h = _two_sum(S, P[i]) - S, s = _two_sum(H, s + h) - sump = sum(p) - H, h = _two_sum(S, sump) - S, s = _fast_two_sum(H, s + h) - res = _acc_sqrt(S, s) - return res - - -def _norm_g(x): - """Compute the compensated Euclidean norm of a vector. - - This function calculates the Euclidean norm (L2 norm) of a vector, - using a compensated algorithm to reduce numerical errors. + x = a + b + return x, (a - x) + b - Parameters - ---------- - x : list of float - The vector whose norm is to be computed. - Returns - ------- - float - The compensated Euclidean norm of the vector. - - Examples - -------- - >>> _norm_g([1.0, 2.0, 3.0]) - 3.7416573867739413 - - Reference - --------- - Stef Graillat, Christoph Lauter, PING Tak Peter Tang, - Naoya Yamanaka, and Shin’ichi Oishi. Efficient Calculations of Faith- - fully Rounded L2-Norms of n-Vectors. ACM Transactions on Mathemat- - ical Software, 41(4), Article 24, 2015. 10.1145/2699469. Available at: - https://doi.org/10.1145/2699469. - """ - S = 0 - s = 0 - for x_i in x: - P, p = _two_prod_fma(x_i, x_i) - H, h = _two_sum(S, P) - c = s + p - d = h + c - S, s = _fast_two_sum(H, d) - res = _acc_sqrt(S, s) - return res +@njit(cache=True, inline="always") +def _sum_non_neg(a_hi, a_lo, b_hi, b_lo): + """Add two non-negative compensated values (mirrors accusphgeom ``sum_non_neg``).""" + hh, h = two_sum(a_hi, b_hi) + d = h + (a_lo + b_lo) + return _fast_two_sum(hh, d) -def _two_square(Aa): - """Compute the square of a number with a compensation for the round-off - error. +@njit(cache=True, inline="always") +def _sum_of_squares_c(hi, lo): + """Compensated squared norm ``Σ (hi[i] + lo[i])²`` of a compensated vector. - This function calculates the square of a given number and compensates - for the round-off error that occurs during the squaring. + Faithful port of ``sum_of_squares_c`` from + accusphgeom/numeric/eft.hpp: a compensated ``Σ hi²`` plus the cross-term + correction ``2·Σ hi·li``. Parameters ---------- - Aa : float - The number to be squared. + hi, lo : tuple of float + Equal-length tuples of the high and low parts of each vector component. + Numba specializes this per tuple length at compile time and keeps the + tuples register-resident (no allocation) — the direct analog of the C++ + template. Used for both nx²+ny² (denominator) and nx²+ny²+nz² (|n|²). Returns ------- tuple of float - The square of the number and the compensated round-off error. - - Examples - -------- - >>> _two_square(2.0) - (4.0, 0.0) - - Reference - --------- - Siegfried Rump. Fast and accurate computation of the Euclidean norm of a vector. J - apan Journal of Industrial and Applied Mathematics, 40, 2023. 10.1007/s13160-023-00593-8 + The compensated squared norm as a ``(hi, lo)`` pair. """ - P = Aa * Aa - A, a = _split(Aa) - p = a * a - ((P - A * A) - 2 * a * A) - return P, p - - -def _acc_sqrt(T, t): - """Compute the accurate square root of a number with a compensation for - round-off error. - - This function calculates the square root of a number, taking into account - a compensation term for the round-off error. + n = len(hi) + s_hi = 0.0 + s_lo = 0.0 + for i in range(n): # compensated Σ hi² + ph, pl = two_prod(hi[i], hi[i]) + s_hi, s_lo = _sum_non_neg(s_hi, s_lo, ph, pl) + r_hi, r_lo = two_prod(hi[0], lo[0]) # accurate Σ hi·li + for i in range(1, n): + p, e = two_prod(hi[i], lo[i]) + r_hi, e2 = two_sum(r_hi, p) + r_lo += e + e2 + return _fast_two_sum(s_hi, (2.0 * (r_hi + r_lo)) + s_lo) + + +@njit(cache=True, inline="always", error_model="numpy") +def acc_sqrt_re(value, error=0.0): + """Accurate square root: return (root, correction) s.t. root+correction ≈ sqrt(value+error). + + Mirrors accusphgeom::numeric::acc_sqrt_re from eft.hpp. Computes + root = fl(sqrt(value)), measures the rounding error of root*root via + two_prod, then recovers a correction term from the residual. When + ``error`` is provided (e.g. the ``lo`` half of a compensated sum), + it is folded into the residual so the correction accounts for the + full compensated input. Parameters ---------- - T : float - The number whose square root is to be computed. - t : float - The compensation term for round-off error. + value : float + Non-negative scalar (the ``hi`` part of a compensated value). + error : float, optional + Low-order correction to ``value`` (default 0.0). Returns ------- - float - The accurate square root of the number. - - Examples - -------- - >>> _acc_sqrt(9.0, 0.0) - 3.0 - - References - ---------- - Vincent Lef`evre, Nicolas Louvet, Jean-Michel Muller, - Joris Picot, and Laurence Rideau. Accurate Calculation of Euclidean - Norms Using Double-Word Arithmetic. ACM Transactions on Mathematical Software, 49(1), 1–34, March 2023. 10.1145/3568672 - - Marko Lange and Siegfried Rump. Faithfully Rounded - Floating-point Computations. ACM Transactions on Mathematical Soft- - ware, 46, 1-20, 2020. 10.1145/3290955 - """ - P = np.sqrt(T) - H, h = _two_square(P) - r = (T - H) - h - r = t + r - p = r / (2 * P) - res = P + p - return res - - -def _split(a): - """Split a floating-point number into two parts: The rounded floating point - presentation and its error. This can be utlized to substitute the FMA - operation on the software level. - - Parameters - ---------- - a : float - The number to be split. - - Returns - ------- - tuple of float - The high and low precision parts of the number. - - Examples - -------- - >>> _split(12345.6789) - (12345.67578125, 0.00311875) - - Reference - --------- - T. J. Dekker. A Floating-Point Technique for Extending the Available Precision. - Numerische Mathematik, 18(3), 224–242, - 1971. 10.1007/BF01397083. Available at: https://doi.org/10.1007/ - BF01397083. - 27 + root : float + Rounded sqrt, fl(sqrt(value)). + correction : float + Additive correction; root + correction ≈ sqrt(value + error) to ~1 ulp. """ - y = (2**27 + 1) * a - x = y - (y - a) - y = a - x - return x, y + # Branch-free, matching AccuSphGeom acc_sqrt_re exactly. Negative value + # yields nan via math.sqrt and root==0 yields nan via the 0/0 correction, + # both under error_model="numpy"; the isfinite mask in the status layer + # rejects such candidates. + root = math.sqrt(value) + sq_hi, sq_lo = two_prod(root, root) + # Residual accumulation order matches AccuSphGeom acc_sqrt_re exactly: + # (value - square.hi) - square.lo + error. + residual = (value - sq_hi) - sq_lo + error + correction = residual / (2.0 * root) + return root, correction