From d351e94e210492d52a0acd599f8508e054526c9c Mon Sep 17 00:00:00 2001 From: nhz2 Date: Sun, 1 Nov 2026 20:09:00 -0500 Subject: [PATCH 1/2] Support exact Rational{BigInt} distances The distance kernels now run on `Rational{BigInt}` coordinates and give exact results, useful as a reference. With floats, degenerate simplexes give NaN or infinite weights that `clamp01nan` makes points of the simplex, but rationals throw on `0//0`. `_clamp01_div` and `_inv0` give zero there instead. Float32 and Float64 results are unchanged. A new test checks the exact results: `d2` is the minimum when every vertex of each simplex is on the far side of the plane through its closest point, perpendicular to `r`. Bump the version to 0.3.1-dev. Co-Authored-By: Claude Opus 5.5 --- Project.toml | 2 +- src/distances.jl | 36 ++++++++++++++++++++++++------------ test/test-mindistance.jl | 34 ++++++++++++++++++++++++++++++++++ 3 files changed, 59 insertions(+), 13 deletions(-) diff --git a/Project.toml b/Project.toml index 934fd34..85ae3bf 100644 --- a/Project.toml +++ b/Project.toml @@ -1,6 +1,6 @@ name = "SimplexCellLists" uuid = "fbe69cd6-f244-4ce4-a90b-8979d86c9ea4" -version = "0.3.0" +version = "0.3.1-dev" authors = ["Nathan Zimmerberg"] [workspace] diff --git a/src/distances.jl b/src/distances.jl index 3939c84..4279aee 100644 --- a/src/distances.jl +++ b/src/distances.jl @@ -7,8 +7,8 @@ Return a tuple `(d², r, s, t)` of: 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. +For floating point coordinates these are approximate, as described under +Accuracy below. For `Rational{BigInt}` coordinates they are exact, which is useful as a reference. `s` holds the barycentric weights of vertices `x[2:N]`; the weight of `x[1]` is implied as `1 - sum(s)`, so the point on `x` is: @@ -108,6 +108,18 @@ end # Ties keep the earlier minimum. A NaN `d2` is kept, so a NaN coordinate gives a NaN `d2`. @inline _min_d2_r_s_t(best, cand) = _select_d2_r_s_t((cand[1] < best[1]) | isnan(cand[1]), cand, best) +# A weight `num/den` clamped to `[0, 1]`. With floats, a degenerate simplex +# gives a NaN or infinite quotient, which `clamp01nan` makes a point of the +# simplex. Exact types such as `Rational` throw on `0//0` instead, so that is +# zero, which is a point of the simplex too. +@inline _clamp01_div(num, den) = clamp01nan(num/den) +@inline _clamp01_div(num::Rational, den::Rational) = iszero(den) ? zero(num) : clamp(num/den, zero(num), one(num)) + +# `inv(x)`, but zero for an exact zero `x` of a type such as `Rational`, +# where `0*inv(0)` throws instead of giving NaN. +@inline _inv0(x) = inv(x) +@inline _inv0(x::Rational) = iszero(x) ? x : inv(x) + # The methods @inline function dist_sqr_r_s_t(x::Point{T}, y::Point{T}) where T @@ -119,7 +131,7 @@ end r = y[2] - y[1] # Written as `y - x[1]` so it folds to `y[1]` when `x[1]` is positive zero. p = y[1] - x[1] - t = clamp01nan(-dot3(p, r)/dot3(r, r)) + t = _clamp01_div(-dot3(p, r), dot3(r, r)) d = fma3(t, r, p) dot3(d, d), d, SVector{0,T}(), SVector{1,T}(t) end @@ -172,12 +184,12 @@ the triangle, and then the edges are exact or within the length of the edge. AX = A - x[1] d2in, rin, u, v = _tri_face_candidate(AX, E0, E1, cross3(E0, E1)) # The closest points on the edges `AC` and `BC`. - u1 = clamp01nan(-dot3(E1, AX)/dot3(E1, E1)) + u1 = _clamp01_div(-dot3(E1, AX), dot3(E1, E1)) r1 = fma3(u1, E1, AX) best = (dot3(r1, r1), r1, zero(T), u1) E2 = E1 - E0 BX = AX + E0 - u2 = clamp01nan(-dot3(BX, E2)/dot3(E2, E2)) + u2 = _clamp01_div(-dot3(BX, E2), dot3(E2, E2)) r2 = fma3(u2, E2, BX) best = _min_d2_r_s_t(best, (dot3(r2, r2), r2, one(T) - u2, u2)) # Last, so a tie with an edge keeps the edge's exact zero weight. @@ -192,10 +204,10 @@ end # already and share `N` between calls. Return `(d2, r, u, v)`, with `r` from `X` to `X'`, `d2 = dot(r, r)`, and `u` # and `v` the weights of `B` and `C`. @inline function _tri_face_candidate(AX, E0, E1, N) - ia = inv(dot3(E0, E0)) + ia = _inv0(dot3(E0, E0)) c1 = dot3(E0, E1)*ia # The weights of `B` and `C` at `X'`, clamped into the triangle, `v` first. - v = clamp01nan(-dot3(N, cross3(E0, AX))/dot3(N, N)) + v = _clamp01_div(-dot3(N, cross3(E0, AX)), dot3(N, N)) u = clamp01nan(fma(-v, c1, -dot3(E0, AX)*ia)) # A select, as `Base.min` also checks for NaN, which can't happen here. w = one(u) - u @@ -260,14 +272,14 @@ Between Line Segments" (Geometric Tools, 2018, revised 2023). n = cross3(P, Q) # `t` of the closest points of the lines through the line segments, # clamped, and the closest point on `x` of the point at `t` on `y`. - t = clamp01nan(-dot3(n, cross3(P, Q0mP0))/dot3(n, n)) - s = clamp01nan(fma(t, b, dot3(P, Q0mP0))/dot3(P, P)) + t = _clamp01_div(-dot3(n, cross3(P, Q0mP0)), dot3(n, n)) + s = _clamp01_div(fma(t, b, dot3(P, Q0mP0)), dot3(P, P)) rin = fma3(-s, P, fma3(t, Q, Q0mP0)) # The closest points on `y` of the end points of `x`. - t0 = clamp01nan(-e/c) + t0 = _clamp01_div(-e, c) r0 = fma3(t0, Q, Q0mP0) best = (dot3(r0, r0), r0, zero(T), t0) - t1 = clamp01nan((b - e)/c) + t1 = _clamp01_div(b - e, c) r1 = fma3(t1, Q, Q0mP0 - P) best = _min_d2_r_s_t(best, (dot3(r1, r1), r1, one(T), t1)) # Last, so a tie with an end point keeps its exact weight. @@ -343,7 +355,7 @@ end @inline function _seg_tri_crossing_candidate(P0, P1, A, E0, E1, N) P = P1 - P0 AP0 = A - P0 - s = clamp01nan(dot3(N, AP0)/dot3(N, P)) + s = _clamp01_div(dot3(N, AP0), dot3(N, P)) d2, r, u, v = _tri_face_candidate(fma3(-s, P, AP0), E0, E1, N) d2, r, s, u, v end diff --git a/test/test-mindistance.jl b/test/test-mindistance.jl index 31cb23d..4d14a6f 100644 --- a/test/test-mindistance.jl +++ b/test/test-mindistance.jl @@ -79,6 +79,40 @@ end end end +# With `Rational{BigInt}` coordinates nothing is rounded, so the results +# should be exact. This checks that `d2` is exactly the minimum squared +# distance, without a reference solver. +# +# Let `X` and `Y` be the points on `a` and `b` that the weights `s` and `t` +# give, and `r = Y - X`. Take the plane through `X` and the plane through `Y`, +# both perpendicular to `r`; they are `‖r‖` apart. If every vertex `p` of `a` +# has `(p - X)⋅r ≤ 0`, all of `a` is on the far side of the first plane, as +# its points are weighted averages of its vertices. If every vertex `q` of `b` +# has `(q - Y)⋅r ≥ 0`, all of `b` is on the far side of the second plane. Then +# no pair of points is closer than `‖r‖`, so `d2 = r⋅r` is the minimum. A +# correct result always passes: if `(p - X)⋅r > 0`, moving `X` a little +# toward `p` would bring it closer to `Y`. +# +# The vertices come from a small pool on a coarse grid, so the simplexes are +# often degenerate, coplanar or parallel. +@testset "Rational{BigInt} exact distance $N to $M" for N in 1:3, M in 1:3 + for trial in 1:300 + pool = [SVector{3, Rational{BigInt}}(rand(-2:2, 3) .// rand(1:2, 3)) for _ in 1:5] + a = SVector{N}(ntuple(_ -> rand(pool), N)) + b = SVector{M}(ntuple(_ -> rand(pool), M)) + d2, r, s, t = dist_sqr_r_s_t(a, b) + X = sum(a .* [1 - sum(s); s;]) + Y = sum(b .* [1 - sum(t); t;]) + @test r == Y - X + @test d2 == r ⋅ r + @test d2 == dist_sqr(a, b) + @test all(p -> (p - X) ⋅ r ≤ 0, a) + @test all(q -> (q - Y) ⋅ r ≥ 0, b) + @test all(w -> 0 ≤ w ≤ 1, s) && sum(s) ≤ 1 + @test all(w -> 0 ≤ w ≤ 1, t) && sum(t) ≤ 1 + end +end + @testset "line to triangle edge case" begin a = SA[ SA[0.7962602087717146, 0.7678176789137129, 0.33365738129157096], From b145690dab0bef2d6abb01133e60539205a7678a Mon Sep 17 00:00:00 2001 From: nhz2 Date: Sun, 1 Nov 2026 20:20:31 -0500 Subject: [PATCH 2/2] Bump version to 0.3.1 --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index 85ae3bf..b659b9c 100644 --- a/Project.toml +++ b/Project.toml @@ -1,6 +1,6 @@ name = "SimplexCellLists" uuid = "fbe69cd6-f244-4ce4-a90b-8979d86c9ea4" -version = "0.3.1-dev" +version = "0.3.1" authors = ["Nathan Zimmerberg"] [workspace]