[Udonburo]

Refs #1604.

This fixes Line::closest_point returning NaN coordinates when a non-zero segment's squared length underflows or the projection's dot products overflow. The existing calculation stays in place unless both dot products round to zero (0/0) or either dot product is non-finite. A non-zero numerator over a zero squared length still produces an infinite projection that clamps to an endpoint.

The fallback scales the direction and query difference separately by powers of two and uses compensated dot products to avoid the cancellation from unit normalization. The tests cover exact and near cancellation, subnormal interior points, endpoint clamping in both directions, and the two suggested overflow regression cases. Both new tests fail with NaN coordinates when the overflow guard is removed and pass when it is restored. Polygon/RTree distance, coordinate subtraction errors and existing intersection-classification issues are outside this change.

Validation after merging main at a4a6087, on Windows x86_64 MSVC with Rust 1.92.0:

• All 45 tests matching closest_point pass in debug and release, including the 27 tests in the closest_point module.
• Formatting, default-features all-target clippy (-Dwarnings), and the all-target build check without default features pass.
• The full geo library suite reports 1,225 passed, 2 failed and 2 ignored. Both failures (test_polygon_densify and haversine_segment_length) reproduce with identical assertion values on unchanged main, using the same lockfile and a separate build directory.
• Doctests, all-features/PROJ validation and other-platform CI were not rerun in this local review.

In the original measurements, ordinary-path timings were similar. The fallback cost about 2–3x for nearby tiny queries and reached about 1.5 µs in the extreme f64 cases measured. These benchmarks were not rerun for this update.

Codex helped with reproduction, implementation, testing, numerical review, merge-conflict resolution and documentation updates.

• I agree to follow the project's code of conduct.
• I added an entry to CHANGES.md.

[urschrei]

This guard catches underflow but not overflow. Past about 1e154 the squared length is +inf, usually with the numerator, so the ordinary path is inf / inf and t is NaN. scaled_projection already pulls an overflowing pair into range, so we only need to add:

Suggested change

	let t = if squared_length == F::zero() && numerator == F::zero() {
	// Avoid 0/0 from squared-length underflow, but preserve non-zero
	// numerators so that infinite projections still clamp to an endpoint.
	let t = if (squared_length == F::zero() && numerator == F::zero())
	|| !squared_length.is_finite()
	|| !numerator.is_finite()
	{
	// Avoid 0/0 from underflow and inf/inf from overflow, but preserve a
	// non-zero numerator over a zero squared length so that infinite
	// projections still clamp to an endpoint.

[urschrei]

Oops, I accidentally committed my own suggested change. Sorry! That having been said, I think this is sound. I think we should hold off to see whether #1609 is reopened and then try this on top.

Suggested change was merged

[urschrei]

Two additional tests from Opus that cover the inf/inf path:

    // Squared length overflows to +inf, so the ordinary path is inf / inf.
    #[test]
    fn long_line_has_finite_closest_points() {
        let line = Line::<f64>::new((0.0, 0.0), (1e200, 1e200));
        let closest = line.closest_point(&Point::new(1e200, 0.0));
        assert_finite(closest);
        let Closest::SinglePoint(closest) = closest else {
            panic!("expected a single point, got {closest:?}");
        };
        assert_relative_eq!(closest.x(), 5e199, max_relative = 1e-15);
        assert_relative_eq!(closest.y(), 5e199, max_relative = 1e-15);
    }

    // 1e160 is finite; its square is not.
    #[test]
    fn line_with_overflowing_squared_length_projects_in_range() {
        let line = Line::<f64>::new((0.0, 0.0), (1e160, 0.0));
        let closest = line.closest_point(&Point::new(5e159, 1.0));
        assert_finite(closest);
        let Closest::SinglePoint(closest) = closest else {
            panic!("expected a single point, got {closest:?}");
        };
        assert_relative_eq!(closest.x(), 5e159, max_relative = 1e-15);
        assert_relative_eq!(closest.y(), 0.0);
    }

[michaelkirk]

Two additional tests from Opus that cover the inf/inf path:

Do you want those two tests added to this branch @urschrei?

[Udonburo]

Hi @urschrei, just wanted to follow up on this PR.

I saw the discussion in #1615 about handling Hegel adoption separately, and was wondering if this fix could move forward independently.

I'd be happy to add the two overflow regression tests you suggested, or any other focused coverage that would be helpful.

Thanks again for taking the time to review this!

[urschrei]

Thanks @Udonburo. To save another review round, I've squashed the branch and added a commit on top that fixes the remaining edge cases.

Changes

• The scaled fallback now runs whenever the squared length is subnormal, not only when it is zero. It uses the same test as line_segment_distance in geo-types. Before this, LINE(0 0,0 3e-162) projected POINT(0 2e-162) to y = 1.5e-162.
• NaN and infinite inputs keep the unscaled quotient. A NaN direction no longer panics in debug builds (attempt to shift right with overflow), and a NaN query gives NaN coordinates, as on main, instead of an endpoint.
• Where the difference of two finite coordinates overflows (LINE(-1e308 0,1e308 0)), the coordinates are halved first and the result is interpolated between the endpoints. closest_point and Euclidean.distance now agree there.
• The projection is rescaled by three exact power-of-two multiplications, replacing a loop of up to about 2,100 iterations.
• The tests use wkt! fixtures. Four new regression tests cover the cases above;
• geo-benches has a new closest_point criterion benchmark.

Performance

Criterion medians, Apple Silicon, cargo bench -p geo-benches --bench closest_point:

Benchmark	main	This PR
Random pairs, f64 (1,024 calls)	4.25 µs	5.25 µs (1.24×)
Random pairs, f32 (1,024 calls)	3.77 µs	5.33 µs (1.41×)
Norway LineString (64 queries)	9.33 ms	9.26 ms (0.99×)
Tiny segments, 1e-162	7.0 µs	19.4 µs
Long segments, 1e200	133 µs	98 µs

Note that while the perf impact looks large, on "real-world" data it's within noise. In the last two rows main returns NaN or wrong results, so those timings don't compare like with like.

Follow-up work we'll want to do

• closest_point and geo-types' line_segment_distance each have their own scaled projection. The trigger and overflow handling now match, but a shared helper would change the distance hot path and needs a geo-types release.
• Segments at extreme scales take 20–30 ns per call in the fallback, mostly because the clamp branches wait on its result.
• Add a hegel property for the full f64 domain asserting that closest_point does not panic, for f32 and f64, once / if hegel is adopted (Ensure Triangle::interior_point is always within triangle bounds (hegel) #1615).

[urschrei]

I've also squashed the commits, mine are currently in a separate commit.

[urschrei]

On the first follow-up item (a shared scaled-projection helper): I measured how far the two fallbacks disagree before writing one, and I don't think it's worth doing.

I compared Euclidean.distance(p, line) with the distance from p to line.closest_point(p) over 2 million random finite pairs. Each segment and query offset had an independent scale between 1e-320 and 1e307.

• Where both results are finite, they agree to within about 3 ulps of the largest coordinate.
• 1.5% of pairs give an infinite distance while closest_point is correct, e.g. LINE(4.6e113 1.98e113,1.22e114 -5.2e112) with POINT(-7.5e214 -6.7e214). These are the point-offset product overflows that Make the point-to-segment distance scale-invariant #1613 leaves unchecked on the raw path. They never reach the distance fallback, so a shared fallback would not change them.

A shared helper would only remove duplicate code. Only the projection parameter could be shared, it would change the distance fallback's last bits, and geo could not be released before a geo-types release. I'm dropping it, and the speed of the closest_point fallback is next.