@delusional/point-proximity

Given a list of candidate points and a target (point or great-circle arc), find the top N closest points or compute all distances on a spherical surface.

Installation

npm install @delusional/point-proximity

Usage

import init, {
  find_top_n,
  compute_all_distances,
  find_top_n_to_line,
  compute_all_distances_to_line,
} from "@delusional/point-proximity";

await init();

// Flat [lat, lng, lat, lng, ...] in degrees
const candidates = new Float64Array([
  51.5074, -0.1278,  // London
  48.8566,  2.3522,  // Paris
  40.7128, -74.0060, // New York
]);

// Top 2 closest to a target point
const top2 = find_top_n(candidates, 50.0, 0.0, 2);
// Float64Array: [idx, dist_km, idx, dist_km]  (sorted closest-first)

// Distance from every point to the target
const dists = compute_all_distances(candidates, 50.0, 0.0);
// Float64Array: [dist_km, dist_km, dist_km]  (same order as input)

// Top 2 closest to a great-circle arc (London → Paris)
const topLine = find_top_n_to_line(
  candidates, 51.5074, -0.1278, 48.8566, 2.3522, 2,
);
// Float64Array: [idx, dist_km, idx, dist_km]

// Distance from every point to the arc
const lineDists = compute_all_distances_to_line(
  candidates, 51.5074, -0.1278, 48.8566, 2.3522,
);
// Float64Array: [dist_km, dist_km, dist_km]

Quick overview

The algorithm converts all points to unit-sphere Cartesian coordinates using a custom SIMD polynomial sin/cos, then scores them against the target using negated dot products — a monotonic proxy for great-circle distance that avoids all asin/sqrt calls during ranking.

Four public functions:

• find_top_n — SIMD-parallel brute-force top-N closest points to a target point
• compute_all_distances — SIMD-parallel distance computation for every point
• find_top_n_to_line — SIMD-parallel top-N closest points to a great-circle arc
• compute_all_distances_to_line — parallel distance to arc for every point

Point proximity

Architecture

Candidates (flat Float64Array)     Target (lat, lng)
        |                                |
        v                                v
  f32x4 polynomial sin/cos         Cartesian f32
  (4 points per SIMD iteration)    (sin/cos once)
        |                                |
  Cartesian f32 conversion               |
  (cos_lat·cos_lng, cos_lat·sin_lng, sin_lat)
        |                                |
        +--- f32x4 negated dot product --+
        |    (3 FMAs, zero libm calls)
        v
  Rayon-parallel scoring
  with per-thread top-N buffers
        |
  Batched compaction
  (select_nth_unstable when full)
        |
  Merge-reduce
  (combine thread results)
        |
  Final sort top-N
        |
  haversine_distance_km (winners only, f64)
        |
  [idx, dist_km, ...]

The scoring trick

For two unit vectors a and b on the sphere, −dot(a, b) is monotonic with great-circle distance: the smaller the dot product, the farther apart they are. This means we can rank all candidates by their dot product with the target and get the exact same ordering as ranking by haversine distance — but without any transcendental functions.

The SIMD kernel computes tx·cx + ty·cy + tz·cz for 4 candidates simultaneously — three multiply-adds on f32x4, zero libm calls. The entire scoring loop is trig-free.

SIMD polynomial sin/cos

The biggest bottleneck in brute-force proximity is converting (lat, lng) pairs to Cartesian coordinates, which requires sin and cos. We replace the libm transcendentals with a Cephes-style minimax polynomial evaluated on f32x4:

• Range reduction: Cody-Waite two-step reduction mod π/2. The high part is the exact f32 representation of π/2; the low part (−4.37e-8) captures the remainder, keeping the reduced argument r accurate to full f32 precision.
• Polynomial: 7th-degree minimax for sin, 6th-degree for cos, both on [−π/4, π/4]. Coefficients are the standard Cephes values.
• Quadrant reconstruction: The integer quotient j = round(x / (π/2)) determines which quadrant we're in. Mask-based swap and negate on i32x4 reconstructs full-range sin/cos from the reduced-range polynomials — no branches.
• Accuracy: ~1.2e-7 max absolute error (full f32 precision). On the unit sphere this is ~0.7 m positional error — more than sufficient for ranking.
• Scalar fallback: fast_sin_cos_f32 provides the same polynomial for single values, used in remainder loops (1–3 leftover points that don't fill a SIMD lane).

This processes 4 trig evaluations per SIMD instruction — ~3–5× faster than 4 sequential scalar f64::sin_cos calls.

Haversine only for winners

The expensive haversine_distance_km conversion (f64, exact) is computed only for the final top-N results. The ranking kernel stays completely libm-free, and the N winners (typically 10–100) are a negligible fraction of the input.

Chord-to-km for all distances

compute_all_distances needs actual km values for every point, so it can't use the ranking-only dot product. Instead it:

1. SIMD Cartesian conversion (same polynomial sin/cos)
2. f64 dot product per point
3. Chord formula: h = (1 − dot)/2, d = 2R · asin(√h)

This is equivalent to haversine but saves ~2 trig calls per point by reusing the SIMD Cartesian coordinates instead of recomputing sin/cos from scratch.

Parallel scoring with batched compaction

The scoring loop is parallelised with Rayon over groups of 4 candidates. Each thread maintains a local Vec<(score, index)> buffer. When it exceeds top_n + buffer_size entries, select_nth_unstable (quickselect, O(n)) compacts it down to the best top_n. This amortises the selection cost across many insertions.

After all threads finish, a reduce step merges thread-local results and keeps the global top-N. Final sort_unstable on only the N winners produces the sorted output.

Buffer sizing: buffer_size = max(top_n, 1024).min(65536). This balances compaction frequency against memory use — for small N we still batch enough to amortise, for large N we don't allocate excessively.

Partial sort

When you need the top 10 out of 10 million, sorting all 10 million is wasteful. Rust's select_nth_unstable finds the k-th smallest element in O(n) average time (quickselect), then we sort only the top-k slice in O(k log k). Total: O(n + k log k) instead of O(n log n).

Line proximity (point-to-arc)

The line functions measure distance from each candidate to the nearest point on a great-circle arc — the shortest path between two points on the sphere. If the candidate's projection onto the great circle falls between the endpoints, the distance is the cross-track distance. Otherwise, it's the distance to the nearer endpoint.

Architecture

Candidates (flat Float64Array)     Arc endpoints (start, end)
        |                                |
        v                                v
  f32x4 polynomial sin/cos         Precompute LineArc:
  (4 points per SIMD iteration)    A, B as Cartesian f32
        |                          n = normalize(A × B)
        v                                |
  P = (cos_lat·cos_lng,                  |
       cos_lat·sin_lng,                  |
       sin_lat)                          |
        |                                |
        +---- SIMD scoring kernel -------+
        |
        v
  For each P (4 at a time):
    ct = dot(P, n)             cross-track to great-circle plane
    Q = P − ct·n               projection onto great-circle plane
    s1 = dot(n, A × Q)         is Q past endpoint A?
    s2 = dot(n, Q × B)         is Q past endpoint B?
        |
    [on-arc?] ──yes──> score = −√(1 − ct²)
        |
        no
        |
    score = −max(dot(P,A), dot(P,B))
        |
        v
  Same parallel fold + batched compaction + merge-reduce
        |
  line_score_to_km_f64 (winners only, f64)
        |
  [idx, dist_km, ...]

The great-circle normal

Given arc endpoints A and B as unit vectors, the great-circle plane has normal n = normalize(A × B). This is precomputed once and splatted into f32x4 for the SIMD kernel.

When |A × B|² < ε (start ≈ end), the arc is degenerate — effectively a single point. The kernel falls back to standard point-proximity scoring against endpoint A.

Cross-track distance

For a point P, the signed distance to the great-circle plane is ct = dot(P, n). Since P is a unit vector and n is a unit normal, |ct| is the sine of the angular distance from P to the great circle. The geodesic cross-track distance is R · arcsin(|ct|).

On-arc test via scalar triple products

The projection of P onto the great-circle plane is Q = P − ct · n. To check whether Q falls between A and B on the arc, we use two scalar triple products:

• s1 = dot(n, A × Q) ≥ 0 — Q is on B's side of A
• s2 = dot(n, Q × B) ≥ 0 — Q is on A's side of B

Both non-negative means Q lies on the arc. This is branchless in SIMD — the final score is selected via on_arc.select(on_arc_score, off_arc_score).

Scoring

• On-arc: −√(1 − ct²). This is −cos(angular_distance), monotonic with distance to the great circle. The sqrt is unavoidable here (unlike point proximity) because the on-arc and off-arc scores must be comparable on the same scale.
• Off-arc: −max(dot(P, A), dot(P, B)). Standard negated dot product to the nearer endpoint — same as point proximity.

The SIMD kernel processes ~97 ops per 4 points: sin/cos for Cartesian conversion, cross product, dot products, projection, triple products, conditional sqrt, and mask-select.

Final km conversion

Since the SIMD scoring kernel doesn't record which case (on-arc vs off-arc) applied for each point, the final km conversion recomputes the full geometry in f64 for the N winners:

• On-arc: R · arcsin(|ct|) (cross-track distance)
• Off-arc: 2R · arcsin(√((1 − max_dot) / 2)) (chord distance to nearest endpoint)

This is only computed for winners (typically N = 10–100), so the cost is negligible.

Data format

Points are passed as a flat Float64Array in [lat, lng, lat, lng, ...] format. No per-point objects, no JSON parsing. The WASM function reads directly from the caller's linear memory — zero-copy ingestion.

Complexity

| Function | Time | Space | |----------|------|-------| | find_top_n | O(n + k log k) | O(k) per thread | | compute_all_distances | O(n) | O(n) | | find_top_n_to_line | O(n + k log k) | O(k) per thread | | compute_all_distances_to_line | O(n) | O(n) |

Where n = number of candidates, k = top_n. Every candidate must be scored (no spatial index), and candidate selection uses select_nth_unstable (introselect, O(n)) with batched compaction. The final k winners are sorted in O(k log k). All functions are parallelised across threads via Rayon.