Repository navigation
Accurate Float32 distances - #24
Merged
Merged
Conversation
Rewrite the point-triangle and line-line kernels from Eberly's expanded quadratic to a split into normal and in-plane parts, with the normal from `cross3_accurate` (Kahan's FMA 2×2 determinant). Expanding the squared distance lost about half the digits in Float32 for large simplexes close together, and the Gram determinant `a*c - b*b` cancels for nearly parallel segments and slivers. Point-triangle computes each barycentric weight from its own edge. Testing `sbar + tbar ≤ 1` judged a Float32 point 0.2 past the short edge of a needle triangle as inside. Inside, the returned weights are capped so they stay a convex combination after rounding. error-analysis/ has a symbolic tracer that runs the kernels on traced numbers and gives first order distance error bounds for all inputs in terms of the input sizes, checked against BigFloat, and a write-up of the point-point, point-line, point-triangle and line-line bounds. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
The point-triangle kernel formed `Δ = ‖E0×E1‖²`, and returned Inf for a point over the face when `inv(Δ)` overflowed, or NaN in the plane, even for a unit sliver. It now takes the closest of three points of the triangle: the closest points on edges AC and BC, and the projection with its weights clamped into the triangle. Each candidate is a point of the triangle, so none is closer than the true distance, and there is no inside test. The weights solve `v` first, as `N⋅(E0×(X - A))/(N⋅N)` with `N = E0×E1`, and then `u` by back substitution against `E1 = c1*E0 + e`, `e ⟂ E0`, so their errors don't grow like one over the sine of an angle. Clamping `v` before solving `u` makes the projection exactly the closest point on edge AB when it is past AB, so AB needs no candidate of its own. A first version solved `e` by Gram-Schmidt with a second pass. It was bound by its chain of dependent operations: every SIMD width took the same ~75 cycles per iteration. The two cross products don't wait on a division, and with `fma3` for the `p + t*r` updates and a `clamp01nan` that compiles to one max and one min, 10k pairs take 4.8 μs in Float32 ×16 SIMD (9.0 for the Gram-Schmidt version) and 9.8 μs in Float64 ×8 (18.7). Also states the contract every `dist_sqr_r_s_t` method must meet, and adds `tiny_normal_cases` to test it where squared lengths and areas are out of range. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
For line segments on the same line to within rounding, `n = P×Q` and `Q0mP0×Q` are both rounding errors, so `sbar` and `tbar` were arbitrary and the inside test could pass. The distance between disjoint segments 0.56 apart then came out as 3e-8, the height along `n`, which is not a distance between points of the segments. For nearly parallel overlapping segments the distance was right, but the weights gave points 0.57 apart with `r` 1.6e-7 long. Integer test coordinates missed this, as an exactly collinear `n` is zero. The kernel now takes the closest of three pairs of points: the closest points on `y` of the end points of `x`, and the closest point on `x` of the point at `t` on `y`, with `t = -n⋅(P×Q0mP0)/(n⋅n)` of the closest points of the lines, clamped. Each is a pair of points of the line segments, so none is closer than the true distance and `r` is the separation of the points the weights give. With `t` clamped to zero or one, `s` is the closest point on `x` of an end point of `y`, so those two candidates are gone. An error in `t` moves the point on `y` by about a rounding across the line through `x`, so it doesn't grow like one over the angle, and plain `cross3` is enough for `n`. 10k pairs take 4.5 μs in Float32 ×16 SIMD (5.9 before) and 9.2 μs in Float64 ×8 (12.4), and a scalar call is about 28% faster. Adds a test of nearly collinear and nearly parallel line segments, which checks the distance and that `r` matches the weights. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
nhz2
changed the title
Nothing uses it since the line-line kernel takes the closest of three pairs of points. The distance is now `dot(r, r)` of a pair of points, accurate whichever pair `n` picks, so `n = P×Q` no longer needs Kahan's 2×2 determinants. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
Möller–Trumbore decided the distance was zero from the signs of triple products, which are rounding errors when a line segment is in the plane of a triangle, so a line segment beside a triangle could get a distance of zero. Both kernels now take the closest of pairs of points of the simplexes: the edges, the vertices with the face, and where an edge meets the plane of the face, branch free. The normal for the crossing needs cross3_accurate, so it is back. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
nhz2
changed the title
State c₁ = 8 and c₂ = 1 in the accuracy bound, say exactly how d² follows from r and how close r is to the separation the weights give, and require a NaN d² for any NaN coordinate. The running minimum now keeps a NaN d². Tests: tolerances down from 64 to 8 eps, the old worst-case table renamed to hard_points with stale comments removed, and new testsets for coordinates near the top of the range and for NaN coordinates. Co-Authored-By: Claude Opus 5.5 <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
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.
Rewrite the line and triangle distance kernels to have better accuracy. The following is the goal of these changes:
What a
dist_sqr_r_s_tmethod must do to be correct, for an element typeT. Heredis the exact distance between thesimplexes as given, and
ℓis the largest absolute coordinate.For finite coordinates below
(floatmax(T)^(1/4))/4(1e9 inFloat32,3e76 in
Float64):NaNorInf, even for degenerate simplexes(repeated vertices, zero length edges, collinear triangles).
sandtis in[0, 1], and the computedsum(s)andsum(t)are at most one.d²is computed fromralone, the same way in every method:fma(r[3], r[3], fma(r[2], r[2], r[1]*r[1])).|√d² - d| ≤ 8*eps(T)*ℓ + floatmin(T)^(1/4).X(s)andY(t)the points the weights give, computed exactly,‖r - (Y(t) - X(s))‖is within the same bound.The
floatmin(T)^(1/4)term (3e-10 inFloat32, 1e-77 inFloat64) is onlyfor when something underflows. Then accuracy may be lost, but only by that
much.
Above
(floatmax(T)^(1/4))/4anything goes, except that for anycoordinates, if one is
NaN,d²isNaN.The methods must also accept
SIMD.Vecelements and be bit identical to the scalar method.tiny_normal_casesandhard_pointsin test/test-mindistance.jl partially checkthis at the edges of the range and at hard points.