Harden spherical routines for near-antipodal inputs (S2 parity)#399
Merged
Conversation
Replace `acos(clamp(x ⋅ y, -1, 1))` with `atan(norm(cross(x, y)), x ⋅ y)`, which stays well-conditioned across the full `[0, π]` range — the acos form collapses to 0 once points are within ~1e-8 of each other. Adapted from Google's S2 `Vector3::Angle`. Add tests for basic correctness, near-identical points (where the naive form fails), near-antipodal points, symmetry, and the two `S2Point` cases from S2's `TEST(S1Angle, ConstructorsThatMeasureAngles)`. Co-Authored-By: Claude Opus 4.7 (1M context) <noreply@anthropic.com>
The `sin(Ω)` divisor in `slerp` collapses for near-antipodal inputs: cancellation-prone geometries drift off the unit sphere, and exactly antipodal inputs return the zero vector. Document basic correctness plus both failure modes, with `@test_broken` on the near-antipodal / antipodal cases so they fire "unexpected pass" once `slerp` is rewritten around S2's robust-cross-product tangent form. Also includes a regression test for ConservativeRegridding.jl#83, where the half-sphere equatorial band produced a truly antipodal diagonal that poisoned SphericalCap radii via a (0,0,0) midpoint. Co-Authored-By: Claude Opus 4.7 (1M context) <noreply@anthropic.com>
Add a new `@testset "s2geometry RobustCrossProd parity"` block that literally reproduces S2's three dedicated `RobustCrossProd` tests (`RobustCrossProdCoverage`, `RobustCrossProdMagnitude`, `RobustCrossProdError`) from `s2edge_crossings_test.cc` at commit `a4f0cf5`, along with the helpers they require (`ChoosePoint`, `PerturbLength`, `TestRobustCrossProdError`). The sign oracle uses a BigFloat 3x3 determinant (`bigsign`), so it is independent of `robust_cross_product`. Straddle-plane probes use that same oracle; tautological perpendicularity substitutes were avoided. The 5000-iteration randomized `RobustCrossProdError` loop aggregates per-iteration failures into a single `@test failures[] == 0` assertion so no failing iterations are silently swallowed. Five `@test_broken` entries register genuine Julia-side behavioral gaps vs S2, confirmed against `s2edge_crossings.cc`: - `RobustCrossProduct.jl:216` drops to the symbolic path on a magnitude threshold (1e-300) where S2 uses literal `IsZero`. - `utils.jl:104` converts to `Float64` before scaling, flushing subnormal BigFloats to zero and defeating the rescaling. - `RobustCrossProduct.jl:327` returns `(-a[2], a[1], 0)` in the `a`-fallback symbolic branch; S2 returns `(a[2], -a[1], 0)` (the inline comment misquotes S2's return value). These bugs compound: the ~80/5000 failures in the randomized loop are all attributable to subnormal inputs routing through the broken symbolic path. Co-Authored-By: Claude Opus 4.7 (1M context) <noreply@anthropic.com>
The parity tests added in the previous commit revealed three places where GeometryOps' port of `S2::RobustCrossProd` silently diverged from S2's reference. Fix them: 1. `RobustCrossProduct.jl:214-222` — the exact path fell through to symbolic when the BigFloat cross had every component below a `1e-300` magnitude threshold. S2's `s2pred::IsZero` is a literal per-component zero test, so subnormal-but-nonzero results like `(0, 0, 5e-324)` were wrongly routed to symbolic. Replace the threshold with per-component `iszero`. 2. `utils.jl:102-141` (`normalizableFromExact`) — the rescaling step converted the BigFloat cross to Float64 *before* applying the scaling multiply, flushing subnormal components to the Float64 subnormal range and destroying the axis. Compute the shift from BigFloat exponents, `ldexp` in BigFloat, then cast. 3. `RobustCrossProduct.jl:326-332` (`symbolic_cross_product_sorted` a-fallback) — returned `(-a[2], a[1], 0)`; S2 returns `(a[1], -a[0], 0)` in 0-indexing, i.e. `(a[2], -a[1], 0)` in 1-indexing. Both components were sign-flipped. The prior in-code comment even misquoted S2's return value as `(-a[1], a[0], 0)`; correct the comment too. Bug 3 also drove the ~80 per-iteration failures in the 5000-iter randomized `TestRobustCrossProdError` loop — those compound with subnormal inputs routed through the broken symbolic path. Co-Authored-By: Claude Opus 4.7 (1M context) <noreply@anthropic.com>
With the three bug fixes in place, the five `@test_broken` entries in the parity block now pass — convert them to `@test` and note the source-fix location in each block's comment. Also tighten two loose assertions in the pre-existing `S2 RobustCrossProdCoverage` testset that were silently masking the same bugs: both `(5e-324, 1, 0) × (0, 1, 0)` and the related subnormal case had been patched with a comment reading "this is what s2 has but we have it on the other axis, IDK why" and the expected axis swapped. Restore S2's expected `(0, 0, 1)`. Bump `bigsign`'s BigFloat precision from 256 to 1024 bits. At 256 bits, the randomized `TestRobustCrossProdError` loop reported ~28 spurious oracle disagreements on near-antipodal iterations — not `robust_cross_product` bugs (the identity / antisymmetry checks held), but float-rounding artifacts in the BigFloat-derived sign. 1024 bits clears the artifacts. Co-Authored-By: Claude Opus 4.7 (1M context) <noreply@anthropic.com>
Replace `sin((1-t)Ω)/sin(Ω) · a + sin(tΩ)/sin(Ω) · b` with `cos(r) · a + sin(r) · dir` where `r = t · spherical_distance(a, b)` and `dir = normalize(robust_cross_product(a, b) × a)`. The two forms are mathematically identical in the well-behaved regime, but the new form never divides by `sin(Ω)` and therefore stays well-conditioned for near-antipodal inputs. For exactly antipodal inputs the great circle is mathematically ambiguous, but `robust_cross_product`'s symbolic-perturbation branch yields a deterministic perpendicular, so `dir` is defined and the result is a unit vector on *some* great circle through both points. Adapted from S2's `Interpolate` / `GetPointOnLine` in `s2edge_distances.cc`. Flip the seven `@test_broken` tripwires in the `slerp` testset to `@test`: cancellation-prone near-antipodal (5 ε values), exactly antipodal, and the ConservativeRegridding#83 antipodal diagonal all now return unit-norm results. The orthogonal midpoint doctest changes by one ulp in the first component (`0.7071067811865475` → `0.7071067811865476`) because the operand order in the tangent-form sum differs from the old formula; both values round to the same real number. Update the jldoctest literal. Co-Authored-By: Claude Opus 4.7 (1M context) <noreply@anthropic.com>
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
1 participant
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
Summary
spherical_distance: switch toatan2(‖x×y‖, x·y)— numerically stable across the full[0, π]range, unlikeacos(x·y)which loses precision for near-identical points. Adapted from S2'sVector3::Angle.robust_cross_product: port S2's threeRobustCrossProdtests (Coverage,Magnitude,Error— the 5000-iter randomized identity/antisymmetry probe) as literal parity tripwires, which surfaced three genuine divergences from S2. Fix all three.slerp: rewrite around S2's tangent-vector formcos(r)·a + sin(r)·dirwithdir = normalize(robust_cross_product(a, b) × a). Mathematically identical to the oldsin((1−t)Ω)/sin(Ω)·a + sin(tΩ)/sin(Ω)·bin the well-behaved regime, but never divides bysin(Ω), so it's well-conditioned for near-antipodal inputs. For exactly antipodal inputs (ambiguous great circle),robust_cross_product's symbolic-perturbation branch yields a deterministic perpendicular; the result is a unit vector on some great circle through both points.Motivating user-visible bug
Regression test for ConservativeRegridding.jl#83: a half-sphere equatorial band produces a truly antipodal diagonal between
(0°, −15°)and(180°, +15°). Under the oldslerpthis returned(0, 0, 0), poisoning theSphericalCapradius incell_range_extentand silently breaking the dual DFS. Now returns a valid unit vector.robust_cross_productbugs fixedAll three were silent — the pre-existing tests had loose assertions or placeholder comments ("this is what s2 has but we have it on the other axis, IDK why") that masked them.
RobustCrossProduct.jl:214-222— exact path used a1e-300magnitude threshold; S2'ss2pred::IsZerois a literal per-component zero test. Non-zero subnormal cross results were wrongly routed to symbolic.utils.jl:102-141(normalizableFromExact) — converted BigFloat to Float64 before the scaling multiply, flushing subnormal axes to the Float64 subnormal range before rescale could preserve them. Fixed by scaling in BigFloat first.RobustCrossProduct.jl:326-332(symbolic_cross_product_sorteda-fallback) — returned(-a[2], a[1], 0); S2 returns(a[1], −a[0], 0)in 0-indexing =(a[2], −a[1], 0)in 1-indexing. Both signs flipped. The in-code comment even misquoted S2's return. This compounded with bugs 1–2 to produce the ~80/5000 per-iteration failures in the randomizedTestRobustCrossProdErrorloop.Tests added
test/utils/unitspherical.jlgrows 118 → ~430 lines of coverage forspherical_distanceandslerp, including near-identical, near-antipodal, exactly antipodal, and the ConservativeRegridding#83 regression case.test/utils/robustcrossproduct.jlgrows by ~430 lines with a new@testset "s2geometry RobustCrossProd parity"block — literal ports of S2's three dedicated tests (CoverageL191,MagnitudeL264,ErrorL321 froms2edge_crossings_test.ccat commita4f0cf5), plus the helpers they require. Uses a BigFloat-determinant sign oracle (bigsign) that is independent ofrobust_cross_productitself.Test plan
test/utils/unitspherical.jl: 889 pass / 1 broken (pre-existingSpherical capstripwire, unrelated) / 0 failtest/utils/robustcrossproduct.jl: all pass / 0 broken / 0 fail (includingfailures[] == 0in the 5000-iter randomized loop)jldoctests forspherical_distanceandslerpupdated where the rewritten formula differs by one ulp🤖 Generated with Claude Code