diff --git a/src/utils/UnitSpherical/point.jl b/src/utils/UnitSpherical/point.jl
index 0bca0e4bb..721e744c0 100644
--- a/src/utils/UnitSpherical/point.jl
+++ b/src/utils/UnitSpherical/point.jl
@@ -114,8 +114,14 @@ GI.crs(::UnitSphericalPoint) = GFT.ProjString("+proj=cart +R=1 +type=crs") # TOD
 """
     spherical_distance(x::UnitSphericalPoint, y::UnitSphericalPoint)
 
-Compute the spherical distance between two points on the unit sphere.  
-Returns a `Number`, usually Float64 but that depends on the input type.
+Compute the great-circle distance (central angle, in radians) between two
+points on the unit sphere.  Returns a `Number` whose type follows from the
+inputs.
+
+Uses the `atan2(‖x × y‖, x · y)` form, which is numerically stable across
+the full `[0, π]` range — unlike `acos(x · y)`, which loses precision for
+nearly-identical points.  Adapted from Google's S2 geometry library
+(see [`Vector3::Angle`](https://github.com/google/s2geometry/blob/a4f0cf58a9cfc214585c39de6e3682384fac0917/src/s2/util/math/vector.h#L492)).
 
 # Extended help
 
@@ -135,7 +141,7 @@ julia> spherical_distance(UnitSphericalPoint(1, 0, 0), UnitSphericalPoint(1, 0,
 0.0
 ```
 """
-spherical_distance(x::UnitSphericalPoint, y::UnitSphericalPoint) = acos(clamp(x ⋅ y, -1.0, 1.0))
+spherical_distance(x::UnitSphericalPoint, y::UnitSphericalPoint) = atan(norm(cross(x, y)), x ⋅ y)
 
 # ## Random points
 Random.rand(rng::Random.AbstractRNG, ::Random.SamplerType{UnitSphericalPoint}) = rand(rng, UnitSphericalPoint{Float64})
diff --git a/src/utils/UnitSpherical/robustcrossproduct/RobustCrossProduct.jl b/src/utils/UnitSpherical/robustcrossproduct/RobustCrossProduct.jl
index 5d71e9e14..e96c76791 100644
--- a/src/utils/UnitSpherical/robustcrossproduct/RobustCrossProduct.jl
+++ b/src/utils/UnitSpherical/robustcrossproduct/RobustCrossProduct.jl
@@ -211,9 +211,13 @@ function exact_cross_product(a::AbstractVector{T}, b::AbstractVector{T}) where T
     big_b = BigFloat.(b; precision=512)
     result_xf = cross(big_a, big_b)
     
-    # Check if we got a non-zero result
-    # This is equivalent to s2's `s2pred::IsZero`.
-    if !all(<=(1e-300), abs.(result_xf))
+    # Check if we got a non-zero result. Matches S2's s2pred::IsZero at
+    # s2predicates_internal.h:77-79, which is a literal sign == 0 test on
+    # each component (i.e., all three components exactly zero). We must NOT
+    # use a magnitude threshold here; any non-zero BigFloat result — even
+    # deeply subnormal by Float64 standards — should go through the exact
+    # rescaling path in `normalizableFromExact`.
+    if !(iszero(result_xf[1]) && iszero(result_xf[2]) && iszero(result_xf[3]))
         return normalizableFromExact(result_xf)
     end
     
@@ -322,9 +326,11 @@ function symbolic_cross_product_sorted(a::AbstractVector{T}, b::AbstractVector{T
     @assert b[1] == 0 && b[2] == 0 "Expected both b[1] and b[2] to be zero"
     
     if a[1] != 0 || a[2] != 0
-        # Fix: This needs to match C++ code which returns (-a[1], a[0], 0) in 0-based indexing
-        # In Julia's 1-based indexing, this is (-a[2], a[1], 0)
-        return UnitSphericalPoint{T}(-a[2], a[1], 0)
+        # S2's SymbolicCrossProdSorted at s2edge_crossings.cc:251-253 returns
+        # `Vector3_d(a[1], -a[0], 0)` in C++ 0-based indexing.
+        # In Julia's 1-based indexing, that's (a[2], -a[1], 0). Bug fixed at
+        # src/utils/UnitSpherical/robustcrossproduct/RobustCrossProduct.jl:327.
+        return UnitSphericalPoint{T}(a[2], -a[1], 0)
     end
     
     # This is always non-zero in the S2 implementation
diff --git a/src/utils/UnitSpherical/robustcrossproduct/utils.jl b/src/utils/UnitSpherical/robustcrossproduct/utils.jl
index 2c4682bcc..2c7ba7107 100644
--- a/src/utils/UnitSpherical/robustcrossproduct/utils.jl
+++ b/src/utils/UnitSpherical/robustcrossproduct/utils.jl
@@ -100,27 +100,45 @@ loss of precision due to floating-point underflow.
 This matches S2's NormalizableFromExact function.
 """
 function normalizableFromExact(xf::AbstractVector{BigFloat})
-    # First try a simple conversion
+    # Mirrors S2's NormalizableFromExact at s2edge_crossings.cc:318-336.
+    #
+    # First try a simple conversion to Float64; if that already has a
+    # component >= 2^-242 it needs no rescaling. Otherwise we must rescale
+    # *in BigFloat* (or equivalently compute the shift from BigFloat
+    # exponents) — we cannot convert to Float64 first, because components
+    # like 5e-324 would flush to the subnormal range or zero, destroying
+    # the axis information before the scaling multiply can fix it.
+    # (Bug fixed at src/utils/UnitSpherical/robustcrossproduct/utils.jl:102.)
     x = Float64.(xf)
-    
+
     if isNormalizable(x)
         return x
     end
-    
-    # Find the largest exponent
-    max_exp = -1000000  # Very small initial value
+
+    # Find the largest BigFloat exponent among nonzero components.
+    found = false
+    max_exp = 0
     for i in 1:3
         if !iszero(xf[i])
-            max_exp = max(max_exp, exponent(xf[i]))
+            e = exponent(xf[i])
+            if !found || e > max_exp
+                max_exp = e
+                found = true
+            end
         end
     end
-    
-    if max_exp < -1000000  # No non-zero components
-        return zero(xf)
+
+    if !found
+        return zero(x)  # The exact result is (0, 0, 0).
     end
-    
-    # Scale to get components in a good range
-    return Float64.(ldexp.(Float64.(xf), -max_exp))
+
+    # Scale in BigFloat so the largest component is in [0.5, 1), then
+    # convert. This matches S2's `ldexp(xf[i], -exp)` inside ExactFloat.
+    return UnitSphericalPoint(
+        Float64(ldexp(xf[1], -max_exp)),
+        Float64(ldexp(xf[2], -max_exp)),
+        Float64(ldexp(xf[3], -max_exp)),
+    )
 end
 
 
diff --git a/src/utils/UnitSpherical/slerp.jl b/src/utils/UnitSpherical/slerp.jl
index 9ef259026..2a65ee128 100644
--- a/src/utils/UnitSpherical/slerp.jl
+++ b/src/utils/UnitSpherical/slerp.jl
@@ -20,9 +20,22 @@ interpolates along that path.
 """
     slerp(a::UnitSphericalPoint, b::UnitSphericalPoint, i01::Number)
 
-Interpolate between `a` and `b`, at a proportion `i01` 
+Interpolate between `a` and `b`, at a proportion `i01`
 between 0 and 1 along the path from `a` to `b`.
 
+Uses the tangent-vector form `cos(r)·a + sin(r)·dir` — where `r = i01 ·
+spherical_distance(a, b)` and `dir = normalize(robust_cross_product(a, b) × a)`
+is the unit tangent at `a` pointing toward `b`. This avoids the `1/sin(Ω)`
+divisor of the classic `sin((1-t)Ω)/sin(Ω) · a + sin(tΩ)/sin(Ω) · b`
+formulation, which collapses for near- and exactly-antipodal inputs.
+Adapted from Google's S2 geometry library (see [`S2::Interpolate`](https://github.com/google/s2geometry/blob/a4f0cf58a9cfc214585c39de6e3682384fac0917/src/s2/s2edge_distances.cc#L77)
+and [`S2::GetPointOnLine`](https://github.com/google/s2geometry/blob/a4f0cf58a9cfc214585c39de6e3682384fac0917/src/s2/s2edge_distances.cc#L47)).
+
+For exactly antipodal `a` and `b` the great circle is mathematically
+ambiguous; `robust_cross_product` returns a deterministic perpendicular via
+its symbolic-perturbation branch, so the result is still a well-defined unit
+vector on *some* great circle through both points.
+
 ## Examples
 
 ```jldoctest
@@ -30,27 +43,30 @@ julia> using GeometryOps.UnitSpherical
 
 julia> slerp(UnitSphericalPoint(1, 0, 0), UnitSphericalPoint(0, 1, 0), 0.5)
 3-element UnitSphericalPoint{Float64} with indices SOneTo(3):
- 0.7071067811865475
+ 0.7071067811865476
  0.7071067811865475
  0.0
 ```
 """
 function slerp(a::UnitSphericalPoint, b::UnitSphericalPoint, i01::Number)
-    if a == b
-        return a
-    end
+    i01 == 0 && return a
+    i01 == 1 && return b
+    a == b && return a
     Ω = spherical_distance(a, b)
-    sinΩ = sin(Ω)
-    return (sin((1-i01)*Ω) / sinΩ) * a + (sin(i01*Ω)/sinΩ) * b
+    dir = normalize(cross(robust_cross_product(a, b), a))
+    r = i01 * Ω
+    return normalize(cos(r) * a + sin(r) * dir)
 end
 
 function slerp(a::UnitSphericalPoint, b::UnitSphericalPoint, i01s::AbstractVector{<: Number})
-    if a == b
-        return fill(a, i01s)
-    end
+    a == b && return fill(a, size(i01s))
     Ω = spherical_distance(a, b)
-    sinΩ = sin(Ω)
-    return @. (sin((1 - i01s) * Ω) / sinΩ) * a + (sin(i01s * Ω) / sinΩ) * b
+    dir = normalize(cross(robust_cross_product(a, b), a))
+    return [begin
+        t == 0 ? a :
+        t == 1 ? b :
+        normalize(cos(t * Ω) * a + sin(t * Ω) * dir)
+    end for t in i01s]
 end
 
 #=
diff --git a/test/utils/robustcrossproduct.jl b/test/utils/robustcrossproduct.jl
index 3f311e365..eb7ad4dd4 100644
--- a/test/utils/robustcrossproduct.jl
+++ b/test/utils/robustcrossproduct.jl
@@ -233,21 +233,22 @@ end
         UnitSphericalPoint(0, 0, 1)
     )
     
-    # Very small perturbation - will use high-precision arithmetic
+    # Very small perturbation - will use high-precision arithmetic.
+    # Tightened to match S2 after fixing the subnormal-flush bug at
+    # src/utils/UnitSpherical/robustcrossproduct/utils.jl:102 and the
+    # IsZero/magnitude-threshold bug at RobustCrossProduct.jl:216.
     test_robust_cross_prod_result(
-        UnitSphericalPoint(5e-324, 1, 0), 
+        UnitSphericalPoint(5e-324, 1, 0),
         UnitSphericalPoint(0, 1, 0),
-        UnitSphericalPoint(1, 0, 0)
-        # UnitSphericalPoint(0, 0, 1) # this is what s2 has but we have it on the other axis, IDK why
+        UnitSphericalPoint(0, 0, 1)
     )
-    
-    # Extremely small differences that can't be represented in double precision
-    # In this case, our implementation may choose a different sign than S2's
+
+    # Extremely small differences that can't be represented in double precision.
+    # Tightened to match S2 (same bugs as above).
     test_robust_cross_prod_result(
-        UnitSphericalPoint(5e-324, 1, 0), 
+        UnitSphericalPoint(5e-324, 1, 0),
         UnitSphericalPoint(5e-324, 1 - DBL_ERR, 0),
-        # UnitSphericalPoint(0, 0, 1)  # We allow either sign in the test function
-        UnitSphericalPoint(1, 0, 0)
+        UnitSphericalPoint(0, 0, 1)
     )
     
     # Test requiring symbolic perturbation
@@ -452,3 +453,417 @@ end
     @test prec_counters[3] + prec_counters[4] > 0  # Some EXACT or SYMBOLIC cases
 end
 
+# =============================================================================
+# s2geometry RobustCrossProd parity
+#
+# Faithful port of three TESTs in S2's s2edge_crossings_test.cc at commit
+# a4f0cf58a9cfc214585c39de6e3682384fac0917:
+#
+#   - RobustCrossProdCoverage (L191-L240):
+#     https://github.com/google/s2geometry/blob/a4f0cf58a9cfc214585c39de6e3682384fac0917/src/s2/s2edge_crossings_test.cc#L191-L240
+#   - RobustCrossProdMagnitude (L264-L284):
+#     https://github.com/google/s2geometry/blob/a4f0cf58a9cfc214585c39de6e3682384fac0917/src/s2/s2edge_crossings_test.cc#L264-L284
+#   - RobustCrossProdError (L321-L347):
+#     https://github.com/google/s2geometry/blob/a4f0cf58a9cfc214585c39de6e3682384fac0917/src/s2/s2edge_crossings_test.cc#L321-L347
+#
+# Plus ported helpers `TestRobustCrossProdError` (L111-L180),
+# `TestRobustCrossProd` (L182-L189), `ChoosePoint` (L289-L304), and
+# `PerturbLength` (L308-L319).
+#
+# The S2 sign oracle is `s2pred::Sign(a, b, c)`, i.e. the sign of the 3x3
+# determinant of (a, b, c) with symbolic-perturbation fallback. We replace it
+# with `bigsign`: the sign of that determinant computed in BigFloat at 1024-bit
+# precision. BigFloat isn't literally exact, but at 1024 bits it is vastly
+# more precise than Float64 and, crucially, is *not* `robust_cross_product`
+# or any predicate that wraps it. When `bigsign` returns 0, the case is
+# genuinely a symbolic one that BigFloat can't resolve; we skip the oracle
+# assertion in that iteration per the brief.
+#
+# Precision note: at the default 256 bits, the RobustCrossProdError loop
+# hit ~28 residual bigsign failures (a/b near-antipodal cases where the
+# determinant is mathematically zero but rounds to a definite sign in 256-bit
+# cross-multiplications). Bumping to 1024 bits makes the oracle strong enough
+# that `det(a, -a, c)` for Float64-derived inputs resolves to exactly 0 via
+# cancellation, eliminating those oracle-precision-artifact failures.
+
+# Independent sign oracle. Computes det(a, b, c) in BigFloat at 1024-bit
+# precision and returns its sign (+1, -1, or 0). Never calls
+# `robust_cross_product`.
+function bigsign(a, b, c)
+    setprecision(BigFloat, 1024) do
+        ab = (BigFloat(a[1]), BigFloat(a[2]), BigFloat(a[3]))
+        bb = (BigFloat(b[1]), BigFloat(b[2]), BigFloat(b[3]))
+        cb = (BigFloat(c[1]), BigFloat(c[2]), BigFloat(c[3]))
+        d = ab[1]*(bb[2]*cb[3] - bb[3]*cb[2]) -
+            ab[2]*(bb[1]*cb[3] - bb[3]*cb[1]) +
+            ab[3]*(bb[1]*cb[2] - bb[2]*cb[1])
+        return d > 0 ? 1 : (d < 0 ? -1 : 0)
+    end
+end
+
+# BigFloat version of s2pred::IsZero(cross(ToExact(a), ToExact(b))). Used to
+# distinguish the "have_exact" case from the symbolic case in `test_rcpe!`.
+function exact_cross_is_zero(a, b)
+    ba = (BigFloat(a[1]; precision=512), BigFloat(a[2]; precision=512), BigFloat(a[3]; precision=512))
+    bb = (BigFloat(b[1]; precision=512), BigFloat(b[2]; precision=512), BigFloat(b[3]; precision=512))
+    c1 = ba[2]*bb[3] - ba[3]*bb[2]
+    c2 = ba[3]*bb[1] - ba[1]*bb[3]
+    c3 = ba[1]*bb[2] - ba[2]*bb[1]
+    return iszero(c1) && iszero(c2) && iszero(c3)
+end
+
+# Port of S2's ChoosePoint:
+# https://github.com/google/s2geometry/blob/a4f0cf58a9cfc214585c39de6e3682384fac0917/src/s2/s2edge_crossings_test.cc#L289-L304
+function s2_choose_point(rng::AbstractRNG)
+    while true
+        x = rand(rng, UnitSphericalPoint)
+        x = x ./ norm(x)
+        coords = ntuple(3) do i
+            v = x[i]
+            if rand(rng) < 0.25           # Denormalized
+                v *= 2.0^(-1022 - 53 * rand(rng))
+            elseif rand(rng) < 1/3        # Zero when squared
+                v *= 2.0^(-511 - 511 * rand(rng))
+            elseif rand(rng) < 0.5        # Simply small
+                v *= 2.0^(-100 * rand(rng))
+            end
+            v
+        end
+        p = UnitSphericalPoint(coords)
+        if sum(abs2, p) >= ldexp(1.0, -968)
+            return UnitSphericalPoint(p ./ norm(p))
+        end
+    end
+end
+
+# Port of S2's PerturbLength:
+# https://github.com/google/s2geometry/blob/a4f0cf58a9cfc214585c39de6e3682384fac0917/src/s2/s2edge_crossings_test.cc#L308-L319
+#
+# S2 uses ExactFloat arithmetic to check the length-squared; we use BigFloat at
+# 256 bits (default precision) — not literally exact, but vastly more precise
+# than Float64 and plenty of margin for the 4*DBL_EPSILON tolerance check.
+function s2_perturb_length(rng::AbstractRNG, p::AbstractVector)
+    scale = 1.0 - 2*eps(Float64) + rand(rng) * 4 * eps(Float64)
+    q = p .* scale
+    bq = (BigFloat(q[1]), BigFloat(q[2]), BigFloat(q[3]))
+    nq2 = bq[1]*bq[1] + bq[2]*bq[2] + bq[3]*bq[3]
+    if abs(nq2 - 1) <= 4 * eps(Float64)
+        return UnitSphericalPoint(q)
+    end
+    return UnitSphericalPoint(p)
+end
+
+# Port of S2's GetPointOnLine (s2edge_distances.cc L47-L53) + GetPointOnRay
+# (s2edge_distances.h L265-L281).
+#
+# S2's version is:
+#   dir_tangent = RobustCrossProd(a, b).CrossProd(a).Normalize()
+#   return (cos(r) * a + sin(r) * dir_tangent).Normalize()
+# We follow the same structure. Note S2 explicitly `.Normalize()`s the result
+# to keep it unit-length despite error from the inexact dot/cross sequence —
+# that final normalize is what keeps the output within `IsUnitLength`.
+function s2_get_point_on_line(a, dir, r)
+    # Tangent at `a` towards `dir`: (a × dir) × a, normalized. This is a
+    # simplification of S2's RobustCrossProd path and is fine here because
+    # `a` and `dir` come from s2_choose_point (normalized).
+    perp = cross(cross(a, dir), a)
+    npp = norm(perp)
+    npp == 0 && return UnitSphericalPoint(a)   # degenerate; caller retries
+    u = perp ./ npp
+    raw = cos(r) .* a .+ sin(r) .* u
+    return UnitSphericalPoint(raw ./ norm(raw))
+end
+
+# Port of S2's TestRobustCrossProdError:
+# https://github.com/google/s2geometry/blob/a4f0cf58a9cfc214585c39de6e3682384fac0917/src/s2/s2edge_crossings_test.cc#L111-L180
+#
+# Each assertion is either registered with `@test` (in `strict=true` mode,
+# used by single-point TestRobustCrossProd calls) or accumulated in the
+# `failures` counter (in `strict=false` mode, used by the 5000-iteration
+# random loop — there the aggregate count is asserted once at the end, so
+# failures are a visible, asserted quantity rather than silenced).
+#
+# Returns the Precision level reached. Julia has no 80-bit long double, so
+# the LONG_DOUBLE branch is absent; we report only DOUBLE, EXACT, or SYMBOLIC.
+function test_rcpe!(a::AbstractVector, b::AbstractVector;
+                    strict::Bool, failures::Ref{Int})
+    kRobustErr = RobustCrossProduct.ROBUST_CROSS_PROD_ERROR  # 6 * DBL_ERR
+
+    result = normalize(robust_cross_product(a, b))
+    check(cond::Bool) = strict ? (@test cond) : (cond || (failures[] += 1))
+
+    # S2 L132-L138. S2's `Sign(a, b, result)` on symbolic cases falls through
+    # to symbolic perturbations and still yields +1/-1. Our BigFloat oracle
+    # returns 0 on such cases (a,b exactly collinear under BigFloat), and
+    # per the brief we skip the probe in that iteration — BigFloat is not a
+    # symbolic-perturbation engine. Even if we skip bigsign, the magnitude
+    # tests later in this function (and the identity/equality checks at
+    # L142-L157) still fire.
+    offset = kRobustErr .* result
+    a90 = cross(result, a)
+    sab = bigsign(a, b, result)
+    if sab != 0
+        # S2 L134: s2pred::Sign(a, b, result) == 1.
+        check(sab == 1)
+        # S2 L135-L138: the probe points (a ± offset) and (a90 ± offset)
+        # straddle the plane through (a, b). Using bigsign as the oracle
+        # instead of S2's `result.DotProd(probe)` gives a check that is
+        # independent of the direction vector under test.
+        #
+        # Limitation: when a bigsign on a probe point returns 0, the probe
+        # lands exactly on the plane (a, b) under BigFloat precision (256
+        # bits). S2's `Sign` resolves this via symbolic perturbation; our
+        # oracle cannot. We treat bigsign==0 as "oracle can't resolve" and
+        # skip the probe (rather than failing): the sab==1 check above
+        # already asserts the main sign, and the negation-identity and
+        # antisymmetry checks below still pin the exact numerical result.
+        sp = bigsign(a, b, a .+ offset);   sp == 0 || check(sp ==  1)
+        sm = bigsign(a, b, a .- offset);   sm == 0 || check(sm == -1)
+        s90p = bigsign(a, b, a90 .+ offset); s90p == 0 || check(s90p ==  1)
+        s90m = bigsign(a, b, a90 .- offset); s90m == 0 || check(s90m == -1)
+    end
+    # NOTE on symbolic (sab == 0) case: we explicitly skip the straddle
+    # probes — BigFloat can't resolve the symbolic tie, and porting S2's
+    # symbolic perturbation would amount to reimplementing
+    # `symbolic_cross_product` (the function under test) as the oracle.
+
+    # S2 L141: "have_exact" is true iff the ExactFloat cross product of a, b
+    # is non-zero.
+    have_exact = !exact_cross_is_zero(a, b)
+
+    # S2 L142-L145: identities under negation (only when non-symbolic).
+    if have_exact
+        check(normalize(robust_cross_product(-a, b)) == -result)
+        check(normalize(robust_cross_product(a, -b)) == -result)
+    end
+
+    # S2 L146-L157: antisymmetry under arg swap, only when a != b.
+    if a != b
+        check(normalize(robust_cross_product(b, a)) == -result)
+    end
+
+    # S2 L163-L179: precision level. Julia has no LD path, so:
+    _, have_dbl = RobustCrossProduct.stable_cross_product(a, b)
+    if have_dbl
+        return DOUBLE
+    elseif have_exact
+        return EXACT
+    else
+        return SYMBOLIC
+    end
+end
+
+# Port of S2's TestRobustCrossProd:
+# https://github.com/google/s2geometry/blob/a4f0cf58a9cfc214585c39de6e3682384fac0917/src/s2/s2edge_crossings_test.cc#L182-L189
+function s2_test_robust_cross_prod(a, b, expected_result, expected_prec::Precision)
+    # S2 L185: Sign(a, b, expected_result) == 1. We use bigsign; if the
+    # expected case is symbolic (bigsign returns 0), the exact-equality
+    # check on the next line still nails the expected direction.
+    sab = bigsign(a, b, expected_result)
+    if sab != 0
+        @test sab == 1
+    end
+    # S2 L186: normalized result equals expected_result *exactly*.
+    @test normalize(robust_cross_product(a, b)) == expected_result
+    # S2 L187: the precision level reached by TestRobustCrossProdError must
+    # equal `expected_prec`. Julia has no LONG_DOUBLE path — treat cases
+    # that S2 expects LONG_DOUBLE for as satisfied by DOUBLE or EXACT.
+    failures = Ref(0)
+    prec = test_rcpe!(a, b; strict=true, failures=failures)
+    if expected_prec == LONG_DOUBLE
+        @test prec == EXACT || prec == DOUBLE
+    else
+        @test prec == expected_prec
+    end
+end
+
+@testset "s2geometry RobustCrossProd parity" begin
+    DBL_ERR_LOCAL = eps(Float64) / 2
+    # Long-double epsilon for platforms with 80-bit LD. Julia doesn't have
+    # 80-bit LD; S2 tests expecting LONG_DOUBLE get accepted as EXACT/DOUBLE.
+    LD_ERR_LOCAL = 2.0^-64 / 2
+
+    # Port of S2 RobustCrossProdCoverage:
+    # https://github.com/google/s2geometry/blob/a4f0cf58a9cfc214585c39de6e3682384fac0917/src/s2/s2edge_crossings_test.cc#L191-L240
+    @testset "RobustCrossProdCoverage" begin
+        # S2 L199-L200.
+        s2_test_robust_cross_prod(
+            UnitSphericalPoint(1.0, 0.0, 0.0),
+            UnitSphericalPoint(0.0, 1.0, 0.0),
+            UnitSphericalPoint(0.0, 0.0, 1.0), DOUBLE)
+        # S2 L201-L202.
+        s2_test_robust_cross_prod(
+            UnitSphericalPoint(20*DBL_ERR_LOCAL, 1.0, 0.0),
+            UnitSphericalPoint(0.0, 1.0, 0.0),
+            UnitSphericalPoint(0.0, 0.0, 1.0), DOUBLE)
+        # S2 L207-L208: 16*DBL_ERR — LONG_DOUBLE on S2 (or EXACT on platforms
+        # w/o LD); in Julia always DOUBLE/EXACT.
+        s2_test_robust_cross_prod(
+            UnitSphericalPoint(16*DBL_ERR_LOCAL, 1.0, 0.0),
+            UnitSphericalPoint(0.0, 1.0, 0.0),
+            UnitSphericalPoint(0.0, 0.0, 1.0), LONG_DOUBLE)
+
+        # S2 L211-L212: 5*LD_ERR — LONG_DOUBLE on S2.
+        s2_test_robust_cross_prod(
+            UnitSphericalPoint(5*LD_ERR_LOCAL, 1.0, 0.0),
+            UnitSphericalPoint(0.0, 1.0, 0.0),
+            UnitSphericalPoint(0.0, 0.0, 1.0), LONG_DOUBLE)
+        # S2 L213-L214: 4*LD_ERR — EXACT.
+        s2_test_robust_cross_prod(
+            UnitSphericalPoint(4*LD_ERR_LOCAL, 1.0, 0.0),
+            UnitSphericalPoint(0.0, 1.0, 0.0),
+            UnitSphericalPoint(0.0, 0.0, 1.0), EXACT)
+
+        # S2 L217-L218: 5e-324 subnormal — EXACT.
+        #
+        # Previously @test_broken. Fixed by:
+        #   - src/utils/UnitSpherical/robustcrossproduct/RobustCrossProduct.jl:216
+        #     (replace 1e-300 magnitude threshold with literal IsZero check)
+        #   - src/utils/UnitSpherical/robustcrossproduct/utils.jl:102
+        #     (scale BigFloat components via ldexp *before* casting to
+        #     Float64, so subnormals aren't flushed to zero)
+        #   - src/utils/UnitSpherical/robustcrossproduct/RobustCrossProduct.jl:327
+        #     (symbolic sign-flip fix, not directly on this path but related)
+        @test normalize(robust_cross_product(
+                        UnitSphericalPoint(5e-324, 1.0, 0.0),
+                        UnitSphericalPoint(0.0, 1.0, 0.0))) ==
+                     UnitSphericalPoint(0.0, 0.0, 1.0)
+
+        # S2 L221-L222: exact cross product underflows in double precision.
+        # Previously @test_broken. Same fixes as above — the BigFloat cross
+        # is (0, 0, ~-5.5e-340), which now goes through the scaled exact
+        # path correctly.
+        @test normalize(robust_cross_product(
+                        UnitSphericalPoint(5e-324, 1.0, 0.0),
+                        UnitSphericalPoint(5e-324, 1 - DBL_ERR_LOCAL, 0.0))) ==
+                     UnitSphericalPoint(0.0, 0.0, -1.0)
+
+        # S2 L225-L226: symbolic.
+        s2_test_robust_cross_prod(
+            UnitSphericalPoint(1.0, 0.0, 0.0),
+            UnitSphericalPoint(1.0 + eps(Float64), 0.0, 0.0),
+            UnitSphericalPoint(0.0, 1.0, 0.0), SYMBOLIC)
+        # S2 L227-L228: symbolic.
+        s2_test_robust_cross_prod(
+            UnitSphericalPoint(0.0, 1.0 + eps(Float64), 0.0),
+            UnitSphericalPoint(0.0, 1.0, 0.0),
+            UnitSphericalPoint(1.0, 0.0, 0.0), SYMBOLIC)
+        # S2 L229-L230: symbolic.
+        s2_test_robust_cross_prod(
+            UnitSphericalPoint(0.0, 0.0, 1.0),
+            UnitSphericalPoint(0.0, 0.0, -1.0),
+            UnitSphericalPoint(-1.0, 0.0, 0.0), SYMBOLIC)
+
+        # S2 L233-L239: symbolic perturbation cases that can't happen in
+        # practice but are implemented for completeness. S2 asserts exact
+        # equality on `SymbolicCrossProd(a,b)`.
+        #
+        # S2 L234-L235. Previously @test_broken due to sign-flip bug at
+        # src/utils/UnitSpherical/robustcrossproduct/RobustCrossProduct.jl:327.
+        # Fix: return `(a[2], -a[1], 0)` (matching S2's 0-based
+        # `(a[1], -a[0], 0)`) instead of the previous `(-a[2], a[1], 0)`.
+        @test RobustCrossProduct.symbolic_cross_product(
+                        UnitSphericalPoint(-1.0, 0.0, 0.0),
+                        UnitSphericalPoint(0.0, 0.0, 0.0)) ==
+                     UnitSphericalPoint(0.0, 1.0, 0.0)
+
+        # S2 L236-L237: same sign-flip path as above, routed through
+        # `b < a` so SymbolicCrossProd returns -sorted(b, a). Fixed by
+        # the same edit at RobustCrossProduct.jl:327.
+        @test RobustCrossProduct.symbolic_cross_product(
+                        UnitSphericalPoint(0.0, 0.0, 0.0),
+                        UnitSphericalPoint(0.0, -1.0, 0.0)) ==
+                     UnitSphericalPoint(1.0, 0.0, 0.0)
+
+        # S2 L238-L239: falls through to the final literal `(1, 0, 0)`
+        # return; this path is not affected by the sign-flip bug.
+        # `SymbolicCrossProd((0,0,0), (0,0,-1))`: not a<b (0>-1 at z),
+        # so `-sorted(b, a)` with inner args (0,0,-1), (0,0,0). All
+        # earlier branches zero; falls to final `(1,0,0)` in sorted,
+        # negated outer: `(-1, 0, 0)`.
+        @test RobustCrossProduct.symbolic_cross_product(
+                UnitSphericalPoint(0.0, 0.0, 0.0),
+                UnitSphericalPoint(0.0, 0.0, -1.0)) ==
+              UnitSphericalPoint(-1.0, 0.0, 0.0)
+    end
+
+    # Port of S2 RobustCrossProdMagnitude:
+    # https://github.com/google/s2geometry/blob/a4f0cf58a9cfc214585c39de6e3682384fac0917/src/s2/s2edge_crossings_test.cc#L264-L284
+    @testset "RobustCrossProdMagnitude" begin
+        # S2 uses EXPECT_DOUBLE_EQ (bitwise double equality) on the angle.
+        # Keep `==` (no tolerance) — this is the load-bearing assertion
+        # that cross-product magnitudes are preserved correctly for
+        # underflowing inputs.
+        c1 = robust_cross_product(
+                UnitSphericalPoint(1.0, 0.0, 0.0),
+                UnitSphericalPoint(1.0, 1e-100, 0.0))
+        c2 = robust_cross_product(
+                UnitSphericalPoint(1.0, 0.0, 0.0),
+                UnitSphericalPoint(1.0, 0.0, 1e-100))
+        @test atan(norm(cross(c1, c2)), dot(c1, c2)) == π/2
+
+        # Same with symbolic perturbations (S2 L279-L283).
+        c1s = robust_cross_product(
+                UnitSphericalPoint(-1e-100, 0.0, 1.0),
+                UnitSphericalPoint(1e-100, 0.0, -1.0))
+        c2s = robust_cross_product(
+                UnitSphericalPoint(0.0, -1e-100, 1.0),
+                UnitSphericalPoint(0.0, 1e-100, -1.0))
+        @test atan(norm(cross(c1s, c2s)), dot(c1s, c2s)) == π/2
+    end
+
+    # Port of S2 RobustCrossProdError:
+    # https://github.com/google/s2geometry/blob/a4f0cf58a9cfc214585c39de6e3682384fac0917/src/s2/s2edge_crossings_test.cc#L321-L347
+    @testset "RobustCrossProdError" begin
+        rng = MersenneTwister(0x5e4f5c4d)
+        iters = 5000
+        counts = Dict(DOUBLE => 0, LONG_DOUBLE => 0, EXACT => 0, SYMBOLIC => 0)
+        failures = Ref(0)
+
+        for _ in 1:iters
+            local a, b
+            while true
+                a = s2_perturb_length(rng, s2_choose_point(rng))
+                dir = s2_choose_point(rng)
+                r = (π/2) * 2.0^(-53 * rand(rng))
+                if rand(rng) < 1/3
+                    r *= 2.0^(-1022 * rand(rng))
+                end
+                b = s2_perturb_length(rng, s2_get_point_on_line(a, dir, r))
+                rand(rng, Bool) && (b = -b)
+                a != b && break
+            end
+            counts[test_rcpe!(a, b; strict=false, failures=failures)] += 1
+        end
+
+        # Tripwire: aggregated count of *any* check failing across all
+        # 5000 iterations must be zero. S2's TestRobustCrossProdError
+        # individually EXPECTs every assertion; we aggregate for
+        # readability and assert the aggregate here.
+        #
+        # Previously @test_broken under the three bugs flagged in
+        # RobustCrossProdCoverage above. Fixed by:
+        #   - src/utils/UnitSpherical/robustcrossproduct/RobustCrossProduct.jl:216
+        #     (IsZero rather than 1e-300 magnitude threshold)
+        #   - src/utils/UnitSpherical/robustcrossproduct/RobustCrossProduct.jl:327
+        #     (sign-flip in symbolic_cross_product_sorted)
+        #   - src/utils/UnitSpherical/robustcrossproduct/utils.jl:102
+        #     (scale BigFloats before Float64 cast in normalizableFromExact)
+        #
+        # Under seed 0x5e4f5c4d this now reaches zero after also bumping
+        # the `bigsign` BigFloat oracle to 1024 bits (see helper above);
+        # at the default 256 bits, ~28 near-antipodal cases still hit
+        # BigFloat-precision artifacts even with the cross-product bugs
+        # fixed.
+        @test failures[] == 0
+        # Always print the count for diagnostics, whether passing or
+        # failing — keeps the number a visible, tracked quantity.
+        @info "RobustCrossProdError aggregate" failures=failures[] iters=iters counts=counts
+
+        # Sanity: the random mix should exercise both DOUBLE and
+        # EXACT/SYMBOLIC paths.
+        @test counts[DOUBLE] > 0
+        @test counts[EXACT] + counts[SYMBOLIC] > 0
+    end
+end
+
diff --git a/test/utils/unitspherical.jl b/test/utils/unitspherical.jl
index 620fdfcc9..5b55e7615 100644
--- a/test/utils/unitspherical.jl
+++ b/test/utils/unitspherical.jl
@@ -5,6 +5,159 @@ using GeometryOps.UnitSpherical
 
 import GeoInterface as GI
 
+@testset "spherical_distance" begin
+    @testset "Basic correctness" begin
+        # Orthogonal axes → π/2
+        @test spherical_distance(UnitSphericalPoint(1.0, 0.0, 0.0),
+                                 UnitSphericalPoint(0.0, 1.0, 0.0)) ≈ π/2
+        # Identical points → 0
+        p = UnitSphericalPoint(1.0, 0.0, 0.0)
+        @test spherical_distance(p, p) == 0.0
+        # Antipodal → π
+        @test spherical_distance(UnitSphericalPoint(1.0, 0.0, 0.0),
+                                 UnitSphericalPoint(-1.0, 0.0, 0.0)) ≈ π
+        # Known π/4 separation
+        @test spherical_distance(UnitSphericalPoint(1.0, 0.0, 0.0),
+                                 UnitSphericalPoint(cos(π/4), sin(π/4), 0.0)) ≈ π/4
+        # Known π/3 separation
+        @test spherical_distance(UnitSphericalPoint(1.0, 0.0, 0.0),
+                                 UnitSphericalPoint(cos(π/3), sin(π/3), 0.0)) ≈ π/3
+    end
+
+    @testset "Near-identical points (atan2 stability vs acos)" begin
+        # For each ε, check that the atan2 form recovers ε with very small
+        # relative error. This is the case where `acos(clamp(x⋅y, -1, 1))`
+        # degrades — x⋅y = cos(ε) ≈ 1 - ε²/2, and double precision can't
+        # distinguish that from 1 once ε ≲ 1e-8.
+        a = UnitSphericalPoint(1.0, 0.0, 0.0)
+        for ε in (1e-4, 1e-6, 1e-8, 1e-10, 1e-12)
+            b = UnitSphericalPoint(cos(ε), sin(ε), 0.0)
+            d = spherical_distance(a, b)
+            @test d ≈ ε rtol=1e-12
+
+            # Contrast: the naive acos form collapses to 0 around ε ≈ 1e-8.
+            d_naive = acos(clamp(a ⋅ b, -1.0, 1.0))
+            if ε ≤ 1e-10
+                @test d_naive == 0.0 || abs(d_naive - ε) / ε > 1e-2
+            end
+        end
+    end
+
+    @testset "Near-antipodal points" begin
+        # Should return ≈ π - ε. The atan2 form stays well-conditioned here
+        # (cross norm ≈ sin(ε) ≈ ε, dot ≈ -1 + ε²/2).
+        a = UnitSphericalPoint(1.0, 0.0, 0.0)
+        for ε in (1e-2, 1e-4, 1e-6, 1e-8)
+            # Point at angle (π - ε) from a, rotated in the xy-plane
+            b = UnitSphericalPoint(cos(π - ε), sin(π - ε), 0.0)
+            @test spherical_distance(a, b) ≈ π - ε rtol=1e-10
+        end
+    end
+
+    # Ported from Google's S2 geometry library, TEST(S1Angle, ConstructorsThatMeasureAngles):
+    # https://github.com/google/s2geometry/blob/a4f0cf58a9cfc214585c39de6e3682384fac0917/src/s2/s1angle_test.cc#L144-L152
+    # S2's Vector3::Angle accepts non-unit vectors (the atan2 form is scale-invariant);
+    # UnitSphericalPoint is typed for unit vectors, so we normalize the (0,0,2) case.
+    @testset "s2geometry ConstructorsThatMeasureAngles" begin
+        @test spherical_distance(UnitSphericalPoint(1.0, 0.0, 0.0),
+                                 UnitSphericalPoint(0.0, 0.0, 1.0)) == π/2
+        @test spherical_distance(UnitSphericalPoint(1.0, 0.0, 0.0),
+                                 UnitSphericalPoint(1.0, 0.0, 0.0)) == 0.0
+    end
+
+    @testset "Symmetry and range" begin
+        using Random
+        Random.seed!(0xA7A7)
+        for _ in 1:50
+            a = UnitSphericalPoint(randn(), randn(), randn())
+            b = UnitSphericalPoint(randn(), randn(), randn())
+            a = a / norm(a); b = b / norm(b)
+            d_ab = spherical_distance(a, b)
+            d_ba = spherical_distance(b, a)
+            @test d_ab == d_ba
+            @test 0 ≤ d_ab ≤ π + 1e-12
+        end
+    end
+end
+
+@testset "slerp" begin
+    @testset "Basic correctness" begin
+        # Quarter-turn along the equator
+        @test slerp(UnitSphericalPoint(1.0, 0.0, 0.0),
+                    UnitSphericalPoint(0.0, 1.0, 0.0), 0.5) ≈
+              UnitSphericalPoint(sqrt(2)/2, sqrt(2)/2, 0.0)
+        # Endpoints
+        a = UnitSphericalPoint(1.0, 0.0, 0.0)
+        b = UnitSphericalPoint(0.0, 1.0, 0.0)
+        @test slerp(a, b, 0.0) ≈ a
+        @test slerp(a, b, 1.0) ≈ b
+        # Identical inputs — fast path
+        @test slerp(a, a, 0.5) == a
+    end
+
+    @testset "Random well-separated points stay on the unit sphere" begin
+        using Random
+        Random.seed!(0x5E5E)
+        for _ in 1:50
+            va = randn(3); va ./= norm(va)
+            vb = randn(3); vb ./= norm(vb)
+            # Skip anything within 1° of antipodal — that's the degenerate regime
+            # covered by the dedicated testsets below.
+            (va ⋅ vb) < -0.9998 && continue
+            a = UnitSphericalPoint(va[1], va[2], va[3])
+            b = UnitSphericalPoint(vb[1], vb[2], vb[3])
+            for t in (0.0, 0.25, 0.5, 0.75, 1.0)
+                @test isapprox(norm(slerp(a, b, t)), 1.0, atol=1e-12)
+            end
+        end
+    end
+
+    @testset "Well-conditioned near-antipodal geometry" begin
+        # When the perpendicular direction is captured cleanly (b lives
+        # exactly in the xz-plane, a on the x-axis), slerp stays accurate
+        # right up against π.
+        a = UnitSphericalPoint(1.0, 0.0, 0.0)
+        for θ_deg in (179.0, 179.99, 179.9999, 179.99999999)
+            θ = deg2rad(θ_deg)
+            b = UnitSphericalPoint(cos(θ), 0.0, sin(θ))
+            @test isapprox(norm(slerp(a, b, 0.5)), 1.0, atol=1e-10)
+        end
+    end
+
+    # Previously broken under the divide-by-sin(Ω) form; now pass via the
+    # tangent-vector form in src/utils/UnitSpherical/slerp.jl.
+    @testset "Cancellation-prone near-antipodal" begin
+        a = UnitSphericalPoint(1.0, 0.0, 0.0)
+        for ε in (1e-10, 1e-12, 1e-14, 1e-15, 1e-16)
+            v = [-1.0, 0.0, ε]; v ./= norm(v)
+            b = UnitSphericalPoint(v[1], v[2], v[3])
+            m = slerp(a, b, 0.5)
+            @test isapprox(norm(m), 1.0, atol=1e-10)
+        end
+    end
+
+    @testset "Exactly antipodal returns a unit vector" begin
+        a = UnitSphericalPoint(1.0, 0.0, 0.0)
+        b = UnitSphericalPoint(-1.0, 0.0, 0.0)
+        m = slerp(a, b, 0.5)
+        @test all(isfinite, m)
+        @test isapprox(norm(m), 1.0, atol=1e-10)
+    end
+
+    # Regression test for ConservativeRegridding.jl#83. A half-sphere
+    # equatorial band has a truly antipodal diagonal between
+    # (0°, −15°) and (180°, +15°); under the old slerp this produced
+    # a (0,0,0) midpoint that poisoned the SphericalCap radius.
+    @testset "ConservativeRegridding#83 antipodal diagonal" begin
+        to_sphere = UnitSphereFromGeographic()
+        p1 = to_sphere((0.0, -15.0))
+        p3 = to_sphere((180.0, +15.0))
+        @test p1 ⋅ p3 ≈ -1.0
+        m = slerp(p1, p3, 0.5)
+        @test isapprox(norm(m), 1.0, atol=1e-10)
+    end
+end
+
 @testset "Coordinate transforms" begin
     @testset "UnitSphereFromGeographic" begin
         # Test with GeoInterface Point