Loosen the distance contract to 64 eps and add the cases over 8

src/distances.jl

@@ -2,19 +2,39 @@
     dist_sqr_r_s_t(x::Simplex{N,T}, y::Simplex{M,T})::Tuple{T, Vec3{T}, SVector{N-1,T}, SVector{M-1,T}}
 
 Return a tuple `(d², r, s, t)` of:
-1. `d²`: the minimum squared distance between the two simplexes.
-2. `r`: the minimum separation vector, pointing from a closest point on `x` to a closest point on `y`.
-3. `s`: the barycentric coordinates of a closest point on `x`.
-4. `t`: the barycentric coordinates of a closest point on `y`.
+1. `d²`: the minimum squared distance between the two simplexes, computed as `r⋅r`.
+2. `r`: the separation vector from a closest point on `x` to a closest point on `y`.
+3. `s`: the barycentric coordinates of that point on `x`.
+4. `t`: the barycentric coordinates of that point on `y`.
+
+These are computed in floating point, so they are approximate, as described
+under Accuracy below.
 
 `s` holds the barycentric weights of vertices `x[2:N]`; the weight of `x[1]`
-is implied as `1 - sum(s)`, so a closest point on `x` is:
+is implied as `1 - sum(s)`, so the point on `x` is:
 
     (1 - sum(s))*x[1] + sum(s[i]*x[i+1] for i in 1:N-1)
 
-Each weight is in `[0, 1]` with `sum(s) ≤ 1`. For a `Point`, `s` is empty;
-for a `LineSeg`, a closest point is `(1 - s[1])*x[1] + s[1]*x[2]`.
-`t` works the same way for `y`.
+For a `Point`, `s` is empty; for a `LineSeg`, the point is
+`(1 - s[1])*x[1] + s[1]*x[2]`. `t` works the same way for `y`.
+
+# Accuracy
+
+With `ℓ` the largest absolute coordinate of `x` and `y`, for finite
+coordinates below `floatmax(T)^(1/4)/4` (1e9 for `Float32`, 3e76 for
+`Float64`):
+
+- `√d²` is within `64*eps(T)*ℓ + floatmin(T)^(1/4)` of the exact minimum distance.
+- Each weight in `s` and `t` is in `[0, 1]`, and the computed `sum(s)` and
+  `sum(t)` are at most one.
+
+The error scales with `ℓ`, not with the distance, so simplexes that touch far
+from the origin can get a small positive `d²`.
+
+With a coordinate above `floatmax(T)^(1/4)/4`, or an infinite one, the results
+may be meaningless.
+
+With a NaN coordinate, `d²` is NaN.
 """
 function dist_sqr_r_s_t end
 
@@ -31,7 +51,7 @@ For finite coordinates below `(floatmax(T)^(1/4))/4` (1e9 in `Float32`,
    `sum(t)` are at most one.
 3. `d²` is computed from `r` alone, the same way in every method:
    `fma(r[3], r[3], fma(r[2], r[2], r[1]*r[1]))`.
-4. `|√d² - d| ≤ 8*eps(T)*ℓ + floatmin(T)^(1/4)`.
+4. `|√d² - d| ≤ 64*eps(T)*ℓ + floatmin(T)^(1/4)`.
 5. With `X(s)` and `Y(t)` the points the weights give, computed exactly,
    `‖r - (Y(t) - X(s))‖` is within the same bound.
 
@@ -53,7 +73,7 @@ this at the edges of the range and at hard points.
 
 Return the minimum squared distance between two simplexes.
 
-This is the `d²` of [`dist_sqr_r_s_t`](@ref), bit for bit.
+This is the `d²` of [`dist_sqr_r_s_t`](@ref), bit for bit, with the same accuracy.
 """
 @inline function dist_sqr(x::Simplex{N, T}, y::Simplex{M, T}) where {T, N, M}
     d2, r, s, t = dist_sqr_r_s_t(x, y)

test/test-mindistance.jl

@@ -264,6 +264,47 @@ end
     end
 end
 
+# A point on the face of a triangle, within 1e-69 of it, with coordinates
+# near 1e7. The distance came out as 9.96, which is 8.07*eps*ℓ, just over 8.
+@testset "Float32 point on the face of a triangle with large coordinates" begin
+    A = Float32.(SA[-10351314, -10323355, -10169968])
+    B = Float32.(SA[ 10351314,   6399932,   9067019])
+    C = Float32.(SA[  9871146,  10323356,  10169969])
+    P = Float32.(SA[  9718146,  10124128,  10000623])
+    x, y = SA[P], SA[A, B, C]
+    tol = 64*eps(Float32)*maximum(v -> maximum(abs, v), (x..., y...))
+    ref = sqrt(dist_sqr(map(v -> big.(v), x), map(v -> big.(v), y)))
+    for (a, b) in ((x, y), (y, x))
+        d2, r, s, t = dist_sqr_r_s_t(a, b)
+        @test abs(sqrt(d2) - ref) ≤ tol
+        @test all(≥(0), s) && all(≥(0), t) && sum(s) ≤ 1 && sum(t) ≤ 1
+        # `r` is the separation of the points the weights give.
+        pa = (1 - sum(big.(s))) .* a[1] .+ sum(big.(s) .* a[2:end]; init = zero(SVector{3, BigFloat}))
+        pb = (1 - sum(big.(t))) .* b[1] .+ sum(big.(t) .* b[2:end]; init = zero(SVector{3, BigFloat}))
+        @test norm(r .- (pb .- pa)) ≤ tol
+    end
+end
+
+# The end point of a long line segment 0.0135 from a short one, with the long
+# one across the range. With the short one first, `t` on the long one from
+# the lines through them, which is 1 - 2e-8, rounded to 3 eps short of one,
+# 9.7 along it. The distance came out as 9.43, which is 8.12*eps*ℓ.
+@testset "Float32 end point of a long line segment near a short one" begin
+    x = SA[SA[1.9877764f6, 9.505106f6, 9.729709f6], SA[-9.696928f6, -5.704686f6, -9.494322f6]]
+    y = SA[SA[-9.676423f6, -5.6547575f6, -9.56445f6], SA[-9.731833f6, -5.7896785f6, -9.374942f6]]
+    tol = 64*eps(Float32)*maximum(v -> maximum(abs, v), (x..., y...))
+    ref = sqrt(dist_sqr(map(v -> big.(v), x), map(v -> big.(v), y)))
+    for (a, b) in ((x, y), (y, x))
+        d2, r, s, t = dist_sqr_r_s_t(a, b)
+        @test abs(sqrt(d2) - ref) ≤ tol
+        @test 0 ≤ s[1] ≤ 1 && 0 ≤ t[1] ≤ 1
+        # `r` is the separation of the points the weights give.
+        pa = (1 - big(s[1])) .* a[1] .+ big(s[1]) .* a[2]
+        pb = (1 - big(t[1])) .* b[1] .+ big(t[1]) .* b[2]
+        @test norm(r .- (pb .- pa)) ≤ tol
+    end
+end
+
 # Hard points: one coordinate carrying the separation, separations that
 # cancel, nearly parallel line segments, and points on the border between
 # two closest features, where rounding picks either. Each is tested as given,