diff --git a/crates/engine/src/cmd/addons/cleanup/tests.rs b/crates/engine/src/cmd/addons/cleanup/tests.rs
index b20eb0d59b..1184687366 100644
--- a/crates/engine/src/cmd/addons/cleanup/tests.rs
+++ b/crates/engine/src/cmd/addons/cleanup/tests.rs
@@ -341,6 +341,33 @@ fn fix_all_is_one_undo_step() {
     assert_eq!(s.doc().unwrap().doc.artboards.len(), doc.artboards.len());
 }
 
+#[test]
+fn a_long_straight_run_loses_its_inner_anchors_in_one_step() {
+    let mut s = session();
+    let n: u32 = 5_000;
+    let mut pts: Vec<(f64, f64)> = (0..n).map(|i| (10.0 + f64::from(i) * 0.05, 150.0)).collect();
+    pts.push((260.0, 190.0));
+    let id = line(&mut s, &pts);
+    let other = rect(&mut s, 300.0, 20.0, 40.0, 40.0);
+    select(&mut s, &[other]);
+    let check = s.execute("addon.cleanup.check", &json!({"kinds": ["collinearAnchors"], "scope": "document"})).unwrap();
+    assert_eq!(check["total"], n - 2, "{check}");
+    assert_eq!(check["findings"][0]["ids"], json!([id.0]));
+    let shown = s.execute("addon.cleanup.show", &json!({"kind": "collinearAnchors", "scope": "document"})).unwrap();
+    assert_eq!(shown["count"], n - 2);
+    assert_eq!(s.doc().unwrap().selection.objects, [id], "Show selects the path");
+    select(&mut s, &[other]);
+    let (steps, before) = (undo_steps(&s), node(&s, id).unwrap());
+    let r = s.execute("addon.cleanup.fix", &json!({"kinds": ["collinearAnchors"], "scope": "document"})).unwrap();
+    assert_eq!((r["total"].as_u64(), r["fixed"]["collinearAnchors"].as_u64()), (Some(u64::from(n) - 2), Some(u64::from(n) - 2)));
+    assert_eq!(undo_steps(&s), steps + 1);
+    assert_eq!(s.doc().unwrap().selection.objects, [other], "Fix keeps the selection");
+    let ends: Vec<(f64, f64)> = node(&s, id).unwrap().path_data().unwrap().subpaths[0].anchors.iter().map(|a| (a.p.x, a.p.y)).collect();
+    assert_eq!(ends, [(10.0, 150.0), (10.0 + f64::from(n - 1) * 0.05, 150.0), (260.0, 190.0)]);
+    s.execute("edit.undo", &json!({})).unwrap();
+    assert_eq!(node(&s, id).unwrap(), before);
+}
+
 #[test]
 fn the_scope_is_the_selection_or_the_document() {
     let mut s = session();
diff --git a/crates/pathops/src/addons/cleanup.rs b/crates/pathops/src/addons/cleanup.rs
index edb3489007..75cf45b583 100644
--- a/crates/pathops/src/addons/cleanup.rs
+++ b/crates/pathops/src/addons/cleanup.rs
@@ -3,7 +3,7 @@
 //! nearly closed paths; [`near_point_groups`] finds anchors to snap together, and [`rejoin`] welds
 //! open pieces whose ends meet back into paths.
 
-use std::collections::HashMap;
+use std::collections::{BTreeSet, HashMap};
 
 use kurbo::{Point, Vec2};
 use vectorcraft_geom::{Anchor, AnchorKind, EPS, SubPath};
@@ -23,24 +23,64 @@ fn moved(p: Point, a_p: Point, h: Point) -> Point {
 /// Remove the corner anchors (no handles, straight segments both sides) where a straight run goes
 /// on in the same direction: the sine of the angle it turns there is at most `tolerance`. Open ends
 /// stay, and a subpath keeps two anchors (three when closed) → how many went.
+///
+/// Each anchor is tested once, first to last, against the nearest anchors still kept on either
+/// side. On a closed subpath, an anchor with none kept before it has the last anchor as its
+/// previous neighbor, and the last anchor has the first one kept as its next neighbor. The kept
+/// anchors move into place once at the end, so the work grows linearly with the number of anchors.
 pub fn remove_collinear_anchors(sp: &mut SubPath, tolerance: f64) -> usize {
-    let mut removed = 0;
-    let mut i = 0;
-    while i < sp.anchors.len() && sp.anchors.len() > min_anchors(sp) {
-        let n = sp.anchors.len();
-        let inner = sp.closed || (i > 0 && i + 1 < n);
-        let (Some(prev), Some(cur), Some(next)) = (sp.anchors.get((i + n - 1) % n), sp.anchors.get(i), sp.anchors.get((i + 1) % n)) else { break };
+    let n = sp.anchors.len();
+    let min = min_anchors(sp);
+    if n <= min {
+        return 0;
+    }
+    let mut keep = vec![true; n];
+    let (mut left, mut first_kept, mut last_kept) = (n, None, None);
+    #[cfg(test)]
+    let mut work = 0;
+    for i in 0..n {
+        if left <= min {
+            break;
+        }
+        #[cfg(test)]
+        {
+            work += 1;
+        }
+        // Every anchor after `i` is still there, the last one included.
+        let p = last_kept.unwrap_or(n - 1);
+        let q = if i + 1 < n { i + 1 } else { first_kept.unwrap_or(i) };
+        let inner = sp.closed || (last_kept.is_some() && i + 1 < n);
+        let (Some(prev), Some(cur), Some(next)) = (sp.anchors.get(p), sp.anchors.get(i), sp.anchors.get(q)) else { break };
         let straight = !prev.has_out() && !cur.has_in() && !cur.has_out() && !next.has_in();
         let (a, b) = (cur.p - prev.p, next.p - cur.p);
         let len = a.hypot() * b.hypot();
         if inner && straight && len > EPS && a.dot(b) > 0.0 && a.cross(b).abs() / len <= tolerance {
-            sp.anchors.remove(i);
-            removed += 1;
+            if let Some(k) = keep.get_mut(i) {
+                *k = false;
+            }
+            left -= 1;
         } else {
-            i += 1;
+            first_kept.get_or_insert(i);
+            last_kept = Some(i);
         }
     }
-    removed
+    if left < n {
+        let mut kept = keep.iter();
+        sp.anchors.retain(|_| kept.next().copied().unwrap_or(true));
+        #[cfg(test)]
+        {
+            work += n;
+        }
+    }
+    #[cfg(test)]
+    COLLINEAR_WORK.with(|w| w.set(work));
+    n - left
+}
+
+#[cfg(test)]
+thread_local! {
+    /// Anchors tested plus anchors moved by the last [`remove_collinear_anchors`] on this thread.
+    static COLLINEAR_WORK: std::cell::Cell<usize> = const { std::cell::Cell::new(0) };
 }
 
 /// Merge each anchor lying within `distance` of the one before it (and a closed subpath's last
@@ -101,58 +141,203 @@ pub fn close_nearly_closed(sp: &mut SubPath, distance: f64) -> bool {
     true
 }
 
+/// The cell of a grid of `size`-wide squares that `p` lies in.
+fn grid_cell(p: Point, size: f64) -> (i64, i64) {
+    ((p.x / size).floor() as i64, (p.y / size).floor() as i64)
+}
+
+/// The cell `(x, y)` and the cells around it, each once.
+fn around((x, y): (i64, i64)) -> Vec<(i64, i64)> {
+    let mut out: Vec<(i64, i64)> = [-1, 0, 1].into_iter().flat_map(|dx| [-1, 0, 1].map(|dy| (x.saturating_add(dx), y.saturating_add(dy)))).collect();
+    out.sort_unstable();
+    out.dedup();
+    out
+}
+
+/// The root of `i`'s set in the union-find forest `parent`.
+fn root(parent: &mut [usize], mut i: usize) -> usize {
+    while let Some(&p) = parent.get(i).filter(|&&p| p != i) {
+        let up = parent.get(p).copied().unwrap_or(p);
+        if let Some(slot) = parent.get_mut(i) {
+            *slot = up;
+        }
+        i = p;
+    }
+    i
+}
+
 /// Groups of `points` to snap together: points within `distance` of one another (and so on, one
 /// to the next) that aren't all at the same place, skipping the pairs `linked` names (such as
 /// neighbouring anchors of one path). Each group lists its points in order; snapping them all to
 /// the first one makes them coincide. Non-finite points and a distance that isn't a positive number
 /// find nothing.
+///
+/// Two points are compared when their cells of a `distance`-wide grid touch. Each group grows from
+/// one point at a time, which looks only at points no group has taken yet, so a cluster of points
+/// all within `distance` of one another is grouped with work that grows linearly with its size.
 pub fn near_point_groups(points: &[Point], distance: f64, linked: impl Fn(usize, usize) -> bool) -> Vec<Vec<usize>> {
     if !(distance > 0.0 && distance.is_finite()) {
         return vec![];
     }
-    let cell = |p: Point| ((p.x / distance).floor() as i64, (p.y / distance).floor() as i64);
-    let mut grid: HashMap<(i64, i64), Vec<usize>> = HashMap::new();
-    for (i, p) in points.iter().enumerate().filter(|(_, p)| p.is_finite()) {
-        grid.entry(cell(*p)).or_default().push(i);
-    }
-    let mut parent: Vec<usize> = (0..points.len()).collect();
-    fn root(parent: &mut [usize], mut i: usize) -> usize {
-        while let Some(&p) = parent.get(i).filter(|&&p| p != i) {
-            let up = parent.get(p).copied().unwrap_or(p);
-            if let Some(slot) = parent.get_mut(i) {
-                *slot = up;
-            }
-            i = p;
+    let mut cells = buckets(points, distance);
+    let mut near = Near {
+        distance,
+        linked,
+        parent: (0..points.len()).collect(),
+        taken: points.iter().map(|p| !p.is_finite()).collect(),
+        queue: vec![],
+        #[cfg(test)]
+        work: 0,
+    };
+    for (start, &p) in points.iter().enumerate() {
+        if near.taken.get(start) != Some(&false) {
+            continue;
         }
-        i
-    }
-    for (&(cx, cy), members) in &grid {
-        for (dx, dy) in [(0, 0), (1, 0), (-1, 1), (0, 1), (1, 1)] {
-            let Some(others) = grid.get(&(cx.saturating_add(dx), cy.saturating_add(dy))) else { continue };
-            for &i in members {
-                for &j in others.iter().filter(|&&j| (dx, dy) != (0, 0) || j > i) {
-                    let (Some(a), Some(b)) = (points.get(i), points.get(j)) else { continue };
-                    let d = a.distance(*b);
-                    if d <= distance && d > 0.0 && !linked(i.min(j), i.max(j)) {
-                        let (ri, rj) = (root(&mut parent, i), root(&mut parent, j));
-                        if let Some(slot) = parent.get_mut(ri.max(rj)) {
-                            *slot = ri.min(rj);
-                        }
-                    }
-                }
-            }
+        if let Some(t) = near.taken.get_mut(start) {
+            *t = true;
+        }
+        let mut at = Some((start, p));
+        while let Some((i, p)) = at {
+            near.reach(&mut cells, i, p);
+            at = near.queue.pop().and_then(|j| Some((j, *points.get(j)?)));
         }
     }
+    #[cfg(test)]
+    NEAR_WORK.with(|w| w.set(near.work));
     // Only pairs that aren't at the same place were joined, so no group is all at one place.
     let mut groups: HashMap<usize, Vec<usize>> = HashMap::new();
     for i in 0..points.len() {
-        groups.entry(root(&mut parent, i)).or_default().push(i);
+        groups.entry(root(&mut near.parent, i)).or_default().push(i);
     }
     let mut out: Vec<Vec<usize>> = groups.into_values().filter(|g| g.len() > 1).collect();
     out.sort_unstable();
     out
 }
 
+#[cfg(test)]
+thread_local! {
+    /// Buckets, places and points looked at by the last [`near_point_groups`] on this thread.
+    static NEAR_WORK: std::cell::Cell<usize> = const { std::cell::Cell::new(0) };
+}
+
+/// The finite points of one cell of a grid that lie in one square a quarter of a cell wide, with
+/// the box around them.
+struct Bucket {
+    lo: Point,
+    hi: Point,
+    /// Each place (bit for bit) and the points there that no group has taken yet, though some of
+    /// them may have been taken since.
+    places: Vec<(Point, Vec<usize>)>,
+}
+
+impl Bucket {
+    /// Whether [`Point::distance`] puts every point here further than `distance` from `p`: one
+    /// coordinate alone differs by more. The square of a `distance` of 1e-150 or more does not
+    /// underflow, so the measured distance is at least that difference.
+    fn beyond(&self, p: Point, distance: f64) -> bool {
+        distance >= 1e-150 && [self.lo.x - p.x, p.x - self.hi.x, self.lo.y - p.y, p.y - self.hi.y].into_iter().any(|g| g > distance)
+    }
+}
+
+/// The finite `points` in buckets, by grid cell.
+fn buckets(points: &[Point], distance: f64) -> HashMap<(i64, i64), Vec<Bucket>> {
+    let quarter = distance / 4.0;
+    let mut cells: HashMap<(i64, i64), Vec<Bucket>> = HashMap::new();
+    // Each bucket by its cell and square, and each place by its cell, bucket and coordinates.
+    let mut bucket_at: HashMap<[i64; 4], usize> = HashMap::new();
+    let mut place_at: HashMap<(usize, [i64; 2], [u64; 2]), usize> = HashMap::new();
+    for (i, &p) in points.iter().enumerate().filter(|(_, p)| p.is_finite()) {
+        let c = grid_cell(p, distance);
+        let list = cells.entry(c).or_default();
+        let square = grid_cell(p, quarter);
+        let b = *bucket_at.entry([c.0, c.1, square.0, square.1]).or_insert_with(|| {
+            list.push(Bucket { lo: p, hi: p, places: vec![] });
+            list.len() - 1
+        });
+        let Some(bucket) = list.get_mut(b) else { continue };
+        bucket.lo = Point::new(bucket.lo.x.min(p.x), bucket.lo.y.min(p.y));
+        bucket.hi = Point::new(bucket.hi.x.max(p.x), bucket.hi.y.max(p.y));
+        let at = *place_at.entry((b, [c.0, c.1], [p.x.to_bits(), p.y.to_bits()])).or_insert_with(|| {
+            bucket.places.push((p, vec![]));
+            bucket.places.len() - 1
+        });
+        if let Some((_, here)) = bucket.places.get_mut(at) {
+            here.push(i);
+        }
+    }
+    cells
+}
+
+/// The grouping behind [`near_point_groups`].
+struct Near<F> {
+    distance: f64,
+    linked: F,
+    /// The union-find forest of the groups.
+    parent: Vec<usize>,
+    /// Whether each point is in a group being grown or grown already (non-finite ones never are).
+    taken: Vec<bool>,
+    /// Points taken into the group being grown whose own neighbors are still to be looked at.
+    queue: Vec<usize>,
+    #[cfg(test)]
+    work: usize,
+}
+
+impl<F: Fn(usize, usize) -> bool> Near<F> {
+    /// Take into `i`'s group every point no group has taken yet that lies within the distance of
+    /// `p`, `i`'s place, in the cells around it, at a different place and not linked to `i`.
+    fn reach(&mut self, cells: &mut HashMap<(i64, i64), Vec<Bucket>>, i: usize, p: Point) {
+        for c in around(grid_cell(p, self.distance)) {
+            let Some(buckets) = cells.get_mut(&c) else { continue };
+            buckets.retain_mut(|b| {
+                #[cfg(test)]
+                {
+                    self.work += 1;
+                }
+                if !b.beyond(p, self.distance) {
+                    b.places.retain_mut(|place| {
+                        self.visit(i, p, place);
+                        !place.1.is_empty()
+                    });
+                }
+                !b.places.is_empty()
+            });
+        }
+    }
+
+    /// Take into `i`'s group the points at `place` that `i` joins, and drop the ones taken.
+    fn visit(&mut self, i: usize, p: Point, (at, here): &mut (Point, Vec<usize>)) {
+        #[cfg(test)]
+        {
+            self.work += 1;
+        }
+        let d = p.distance(*at);
+        if !(d <= self.distance && d > 0.0) {
+            return;
+        }
+        here.retain(|&j| {
+            if self.taken.get(j) != Some(&false) {
+                return false;
+            }
+            #[cfg(test)]
+            {
+                self.work += 1;
+            }
+            if (self.linked)(i.min(j), i.max(j)) {
+                return true;
+            }
+            let (ri, rj) = (root(&mut self.parent, i), root(&mut self.parent, j));
+            if let Some(slot) = self.parent.get_mut(ri.max(rj)) {
+                *slot = ri.min(rj);
+            }
+            if let Some(t) = self.taken.get_mut(j) {
+                *t = true;
+            }
+            self.queue.push(j);
+            false
+        });
+    }
+}
+
 /// Open pieces welded into one path by [`rejoin`]: `members` index the pieces in the order they
 /// run along `path`.
 #[derive(Clone, Debug, PartialEq)]
@@ -179,6 +364,12 @@ fn outward(sp: &SubPath, last: bool) -> Option<Vec2> {
 /// Each junction keeps the earlier piece's anchor and the later piece's out handle. Pieces that are
 /// closed, have fewer than two anchors or a non-finite coordinate stay out. → the chains of two or
 /// more pieces, and the single pieces that closed.
+///
+/// Two ends are compared when their cells of a `gap`-wide grid touch, and pairs are taken in order
+/// of the turn, then the distance, then the order of the ends. Each end looks for its first
+/// partner in that order among the ends still free, searching outward from the reverse of its own
+/// direction, and two ends that are each other's first partner meet; these are the pairs that
+/// taking every pair in that order gives.
 pub fn rejoin(pieces: &[SubPath], gap: f64) -> Vec<Rejoined> {
     if !(gap >= 0.0 && gap.is_finite()) {
         return vec![];
@@ -187,7 +378,7 @@ pub fn rejoin(pieces: &[SubPath], gap: f64) -> Vec<Rejoined> {
         !sp.closed && sp.anchors.len() >= 2 && sp.anchors.iter().all(|a| a.p.is_finite() && a.h_in.is_finite() && a.h_out.is_finite())
     };
     // Each usable end: where it is and the way out of its piece there.
-    let ends: Vec<(End, Point, Vec2)> = pieces
+    let ends: Vec<EndAt> = pieces
         .iter()
         .enumerate()
         .filter(|(_, sp)| usable(sp))
@@ -197,49 +388,254 @@ pub fn rejoin(pieces: &[SubPath], gap: f64) -> Vec<Rejoined> {
             Some(((i, last), a.p, outward(sp, last)?))
         })
         .collect();
-    let points: Vec<Point> = ends.iter().map(|e| e.1).collect();
-    // Every pair of ends within the gap, by the turn across the junction.
-    let mut pairs: Vec<(f64, f64, End, End)> = vec![];
-    for (a, b) in candidate_pairs(&points, gap.max(1e-9)) {
-        let (Some(&(ea, pa, oa)), Some(&(eb, pb, ob))) = (ends.get(a), ends.get(b)) else { continue };
-        // A piece's own two ends meet only when closing it leaves a shape (three anchors).
-        if ea.0 == eb.0 && pieces.get(ea.0).is_none_or(|sp| sp.anchors.len() < 3) {
-            continue;
-        }
-        let d = pa.distance(pb);
-        if d <= gap {
-            pairs.push((oa.dot(-ob).clamp(-1.0, 1.0).acos(), d, ea, eb));
-        }
-    }
-    pairs.sort_by(|x, y| x.0.total_cmp(&y.0).then(x.1.total_cmp(&y.1)));
-    let mut link: HashMap<End, End> = HashMap::new();
-    for (_, _, a, b) in pairs {
-        if !link.contains_key(&a) && !link.contains_key(&b) {
-            link.insert(a, b);
-            link.insert(b, a);
-        }
-    }
+    let link = Weld::new(pieces, &ends, gap).link();
     chains(pieces, &link)
 }
 
-/// Pairs (a < b) of `points` that may lie within `gap`: those in neighbouring grid cells.
-fn candidate_pairs(points: &[Point], gap: f64) -> Vec<(usize, usize)> {
-    let cell = |p: Point| ((p.x / gap).floor() as i64, (p.y / gap).floor() as i64);
-    let mut grid: HashMap<(i64, i64), Vec<usize>> = HashMap::new();
-    for (i, p) in points.iter().enumerate() {
-        grid.entry(cell(*p)).or_default().push(i);
-    }
-    let mut out = vec![];
-    for (i, p) in points.iter().enumerate() {
-        let (cx, cy) = cell(*p);
-        for dx in -1..=1 {
-            for dy in -1..=1 {
-                let Some(others) = grid.get(&(cx.saturating_add(dx), cy.saturating_add(dy))) else { continue };
-                out.extend(others.iter().filter(|&&j| j > i).map(|&j| (i, j)));
+/// An end, where it lies and the unit direction pointing out of its piece there.
+type EndAt = (End, Point, Vec2);
+
+/// Where a pair of ends (`a` < `b` in the list of ends) comes in the order [`rejoin`] takes pairs:
+/// the turn across the junction, the distance between the ends, `a`, the place of `b`'s grid cell
+/// among the cells around `a`'s, and `b`.
+type Order = (f64, f64, usize, usize, usize);
+
+/// Whether pair `x` comes before pair `y`.
+fn earlier(x: &Order, y: &Order) -> bool {
+    x.0.total_cmp(&y.0).then(x.1.total_cmp(&y.1)).then((x.2, x.3, x.4).cmp(&(y.2, y.3, y.4))).is_lt()
+}
+
+/// The place of cell `to` among the cells around `from`, listed by x step and then y step from -1
+/// to 1 (cells past the end of the grid repeat their neighbor; the first place counts).
+fn listed_at(from: (i64, i64), to: (i64, i64)) -> usize {
+    let mut steps = [-1, 0, 1].into_iter().flat_map(|dx| [-1, 0, 1].map(|dy| (from.0.saturating_add(dx), from.1.saturating_add(dy))));
+    steps.position(|c| c == to).unwrap_or(9)
+}
+
+/// `angle` as an integer in the same order.
+fn angle_key(angle: f64) -> i64 {
+    let bits = angle.to_bits() as i64;
+    bits ^ (((bits >> 63) as u64) >> 1) as i64
+}
+
+/// The angle between the directions at angles `a` and `b`, from 0 to π.
+fn apart(a: f64, b: f64) -> f64 {
+    let d = (a - b).abs();
+    d.min(std::f64::consts::TAU - d)
+}
+
+/// How far past the turn of the best pair found so far [`Weld::first_partner`] still looks. The
+/// turn [`rejoin`] measures (the arc cosine of a dot product of unit vectors) is within 1e-7 of
+/// the angle between the directions, and the angles the search sorts by are within 1e-14 of it.
+const TURN_SLACK: f64 = 1e-6;
+
+/// Ends in one grid cell at one place with one outward direction (bit for bit): every other end
+/// pairs with each of them at the same turn and distance.
+struct Class {
+    at: Point,
+    angle: f64,
+    /// The ends here no pair has taken yet.
+    free: BTreeSet<usize>,
+}
+
+/// The pairing behind [`rejoin`].
+struct Weld<'a> {
+    pieces: &'a [SubPath],
+    ends: &'a [EndAt],
+    gap: f64,
+    /// The width of the grid's cells.
+    size: f64,
+    classes: Vec<Class>,
+    /// Each end's class.
+    class_of: Vec<usize>,
+    /// Each grid cell's classes with free ends, by the angle of their direction.
+    cells: HashMap<(i64, i64), BTreeSet<(i64, usize)>>,
+    /// Whether each end is still free.
+    free: Vec<bool>,
+    /// Whether each end has looked for its first partner. Another end's first partner in a class
+    /// is one of the class's first three free ends, so the ends after those cannot meet next and
+    /// look only once they move up.
+    awake: Vec<bool>,
+    /// Each awake end's first partner, and where that pair comes in the order.
+    first: Vec<Option<(Order, usize)>>,
+    /// The ends whose first partner each end is (some may have moved on since).
+    wanted_by: Vec<Vec<usize>>,
+    /// Ends to check for a partner whose first partner they are.
+    todo: Vec<usize>,
+    #[cfg(test)]
+    work: std::cell::Cell<usize>,
+}
+
+impl<'a> Weld<'a> {
+    fn new(pieces: &'a [SubPath], ends: &'a [EndAt], gap: f64) -> Self {
+        let size = gap.max(1e-9);
+        let (mut classes, mut class_of, mut cells) = (vec![], vec![], HashMap::new());
+        let mut class_at: HashMap<((i64, i64), [u64; 4]), usize> = HashMap::new();
+        for (i, &(_, p, out)) in ends.iter().enumerate() {
+            let c = grid_cell(p, size);
+            let k = *class_at.entry((c, [p.x, p.y, out.x, out.y].map(f64::to_bits))).or_insert_with(|| {
+                let angle = out.y.atan2(out.x);
+                cells.entry(c).or_insert_with(BTreeSet::new).insert((angle_key(angle), classes.len()));
+                classes.push(Class { at: p, angle, free: BTreeSet::new() });
+                classes.len() - 1
+            });
+            if let Some(class) = classes.get_mut(k) {
+                class.free.insert(i);
             }
+            class_of.push(k);
+        }
+        let n = ends.len();
+        Weld {
+            pieces,
+            ends,
+            gap,
+            size,
+            classes,
+            class_of,
+            cells,
+            free: vec![true; n],
+            awake: vec![false; n],
+            first: vec![None; n],
+            wanted_by: vec![vec![]; n],
+            todo: vec![],
+            #[cfg(test)]
+            work: std::cell::Cell::new(0),
         }
     }
-    out
+
+    /// The links between the ends that meet: two free ends that are each other's first partner
+    /// meet, and the ends that wanted either of them look again, until no such pair is left.
+    fn link(mut self) -> HashMap<End, End> {
+        for k in 0..self.classes.len() {
+            self.wake(k);
+        }
+        let mut link: HashMap<End, End> = HashMap::new();
+        while let Some(v) = self.todo.pop() {
+            let Some(w) = self.partner(v) else { continue };
+            if self.free.get(v) != Some(&true) || self.free.get(w) != Some(&true) || self.partner(w) != Some(v) {
+                continue;
+            }
+            let (Some(&(ev, ..)), Some(&(ew, ..))) = (self.ends.get(v), self.ends.get(w)) else { continue };
+            link.insert(ev, ew);
+            link.insert(ew, ev);
+            self.take(v);
+            self.take(w);
+            for e in [v, w] {
+                if let Some(&k) = self.class_of.get(e) {
+                    self.wake(k);
+                }
+            }
+            let waiting: Vec<usize> = [v, w].iter().flat_map(|&e| self.wanted_by.get_mut(e).map(std::mem::take).unwrap_or_default()).collect();
+            for x in waiting {
+                if self.free.get(x) == Some(&true) && matches!(self.partner(x), Some(p) if p == v || p == w) {
+                    self.look(x);
+                }
+            }
+        }
+        #[cfg(test)]
+        REJOIN_WORK.with(|w| w.set(self.work.get()));
+        link
+    }
+
+    /// End `v`'s first partner, once it has looked.
+    fn partner(&self, v: usize) -> Option<usize> {
+        self.first.get(v).and_then(|f| f.as_ref()).map(|f| f.1)
+    }
+
+    /// Have the first three free ends of class `k` that have not looked for a partner look.
+    fn wake(&mut self, k: usize) {
+        let Some(class) = self.classes.get(k) else { return };
+        let sleeping: Vec<usize> = class.free.iter().take(3).copied().filter(|&e| self.awake.get(e) == Some(&false)).collect();
+        for e in sleeping {
+            if let Some(a) = self.awake.get_mut(e) {
+                *a = true;
+            }
+            self.look(e);
+        }
+    }
+
+    /// Find free end `e`'s first partner, and queue `e` to check whether that end wants it too.
+    fn look(&mut self, e: usize) {
+        let next = self.first_partner(e);
+        if let Some(list) = next.as_ref().and_then(|f| self.wanted_by.get_mut(f.1)) {
+            list.push(e);
+        }
+        if let Some(f) = self.first.get_mut(e) {
+            *f = next;
+        }
+        self.todo.push(e);
+    }
+
+    /// End `e` has met another: it is no longer free.
+    fn take(&mut self, e: usize) {
+        if let Some(f) = self.free.get_mut(e) {
+            *f = false;
+        }
+        let Some(&k) = self.class_of.get(e) else { return };
+        let Some(class) = self.classes.get_mut(k) else { return };
+        class.free.remove(&e);
+        if class.free.is_empty()
+            && let Some(&(_, p, _)) = self.ends.get(e)
+            && let Some(set) = self.cells.get_mut(&grid_cell(p, self.size))
+        {
+            set.remove(&(angle_key(class.angle), k));
+        }
+    }
+
+    /// The free end `v` meets first, and where that pair comes in the order: the classes in the
+    /// cells around `v`'s are looked at by how far their direction is from the reverse of `v`'s,
+    /// until that is more than [`TURN_SLACK`] past the best turn found.
+    fn first_partner(&self, v: usize) -> Option<(Order, usize)> {
+        let &((piece, _), p, out) = self.ends.get(v)?;
+        let toward = (-out.y).atan2(-out.x);
+        let split = (angle_key(toward), 0);
+        let mut best: Option<(Order, usize)> = None;
+        for c in around(grid_cell(p, self.size)) {
+            let Some(set) = self.cells.get(&c) else { continue };
+            // Both ways round from `toward`; between them they meet every class once.
+            let mut up = set.range(split..).chain(set.range(..split)).peekable();
+            let mut down = set.range(..split).rev().chain(set.range(split..).rev()).peekable();
+            let turn = |e: Option<&&(i64, usize)>| e.and_then(|e| self.classes.get(e.1)).map_or(f64::INFINITY, |k| apart(k.angle, toward));
+            for _ in 0..set.len() {
+                let (tu, td) = (turn(up.peek()), turn(down.peek()));
+                if best.as_ref().is_some_and(|b| tu.min(td) > b.0.0 + TURN_SLACK) {
+                    break;
+                }
+                let Some(&(_, k)) = (if tu <= td { up.next() } else { down.next() }) else { break };
+                if let Some(pair) = self.pair_with(v, p, piece, k)
+                    && best.as_ref().is_none_or(|b| earlier(&pair.0, &b.0))
+                {
+                    best = Some(pair);
+                }
+            }
+        }
+        best
+    }
+
+    /// The first free end of class `k` that free end `v` (at `p`, on `piece`) may meet, and where
+    /// that pair comes in the order. Within a class the pairs differ only in the ends' order, so
+    /// the end listed first comes first.
+    fn pair_with(&self, v: usize, p: Point, piece: usize, k: usize) -> Option<(Order, usize)> {
+        #[cfg(test)]
+        self.work.set(self.work.get() + 1);
+        let class = self.classes.get(k)?;
+        if p.distance(class.at) > self.gap {
+            return None;
+        }
+        // A piece's own two ends meet only when closing it leaves a shape (three anchors).
+        let short = self.pieces.get(piece).is_none_or(|sp| sp.anchors.len() < 3);
+        let u = *class.free.iter().find(|&&u| u != v && !(short && self.ends.get(u).is_some_and(|e| e.0.0 == piece)))?;
+        let (a, b) = (v.min(u), v.max(u));
+        let (&(_, pa, oa), &(_, pb, ob)) = (self.ends.get(a)?, self.ends.get(b)?);
+        let at = listed_at(grid_cell(pa, self.size), grid_cell(pb, self.size));
+        Some(((oa.dot(-ob).clamp(-1.0, 1.0).acos(), pa.distance(pb), a, at, b), u))
+    }
+}
+
+#[cfg(test)]
+thread_local! {
+    /// Classes of ends looked at by the last [`rejoin`] on this thread.
+    static REJOIN_WORK: std::cell::Cell<usize> = const { std::cell::Cell::new(0) };
 }
 
 /// Walk the links between piece ends into chains.
diff --git a/crates/pathops/src/addons/cleanup/tests.rs b/crates/pathops/src/addons/cleanup/tests.rs
index 3e1c58d42b..d262128234 100644
--- a/crates/pathops/src/addons/cleanup/tests.rs
+++ b/crates/pathops/src/addons/cleanup/tests.rs
@@ -1,7 +1,12 @@
-use kurbo::{Point, Vec2};
-use vectorcraft_geom::{Anchor, AnchorKind, SubPath};
+use std::collections::HashMap;
 
-use super::{close_nearly_closed, merge_duplicate_anchors, near_point_groups, rejoin, remove_collinear_anchors, retract_tiny_handles};
+use kurbo::{Point, Vec2};
+use vectorcraft_geom::{Anchor, AnchorKind, EPS, SubPath};
+
+use super::{
+    End, Rejoined, chains, close_nearly_closed, merge_duplicate_anchors, near_point_groups, outward, rejoin, remove_collinear_anchors,
+    retract_tiny_handles, root,
+};
 
 fn line(pts: &[(f64, f64)], closed: bool) -> SubPath {
     SubPath::polyline(&pts.iter().map(|&p| p.into()).collect::<Vec<Point>>(), closed)
@@ -30,6 +35,174 @@ fn collinear_anchors_go_and_corners_stay() {
     assert_eq!(remove_collinear_anchors(&mut bend, 0.005), 1);
 }
 
+/// The collinear-anchor removal done one `Vec::remove` at a time, testing the anchor that moves
+/// into the gap next: the result [`remove_collinear_anchors`] gives.
+fn one_at_a_time(sp: &mut SubPath, tolerance: f64) -> usize {
+    let min = if sp.closed { 3 } else { 2 };
+    let mut removed = 0;
+    let mut i = 0;
+    while i < sp.anchors.len() && sp.anchors.len() > min {
+        let n = sp.anchors.len();
+        let inner = sp.closed || (i > 0 && i + 1 < n);
+        let (prev, cur, next) = (sp.anchors[(i + n - 1) % n], sp.anchors[i], sp.anchors[(i + 1) % n]);
+        let straight = !prev.has_out() && !cur.has_in() && !cur.has_out() && !next.has_in();
+        let (a, b) = (cur.p - prev.p, next.p - cur.p);
+        let len = a.hypot() * b.hypot();
+        if inner && straight && len > EPS && a.dot(b) > 0.0 && a.cross(b).abs() / len <= tolerance {
+            sp.anchors.remove(i);
+            removed += 1;
+        } else {
+            i += 1;
+        }
+    }
+    removed
+}
+
+/// A subpath's closed flag and every anchor's coordinates bit for bit, with its kind.
+fn bits(sp: &SubPath) -> (bool, Vec<([u64; 6], AnchorKind)>) {
+    let anchors = sp.anchors.iter().map(|a| ([a.p.x, a.p.y, a.h_in.x, a.h_in.y, a.h_out.x, a.h_out.y].map(f64::to_bits), a.kind));
+    (sp.closed, anchors.collect())
+}
+
+/// Check that [`remove_collinear_anchors`] gives `sp` the count and anchors [`one_at_a_time`] does.
+fn same_as_one_at_a_time(sp: &SubPath, tolerance: f64) {
+    let (mut fast, mut slow) = (sp.clone(), sp.clone());
+    let got = remove_collinear_anchors(&mut fast, tolerance);
+    let want = one_at_a_time(&mut slow, tolerance);
+    assert_eq!(got, want, "count at tolerance {tolerance} for {sp:?}");
+    assert_eq!(bits(&fast), bits(&slow), "anchors at tolerance {tolerance} for {sp:?}");
+}
+
+/// xorshift64: the same numbers on every run.
+struct Rng(u64);
+
+impl Rng {
+    fn next(&mut self) -> u64 {
+        self.0 ^= self.0 << 13;
+        self.0 ^= self.0 >> 7;
+        self.0 ^= self.0 << 17;
+        self.0
+    }
+
+    fn pick<T: Copy>(&mut self, xs: &[T]) -> T {
+        xs[(self.next() % xs.len() as u64) as usize]
+    }
+}
+
+/// A subpath of `n` anchors stepping mostly straight along one direction, with repeated points,
+/// small and large bends, reversals, short and long handles and non-finite or huge coordinates.
+fn generated(r: &mut Rng, n: usize, closed: bool) -> SubPath {
+    let dir = r.pick(&[Vec2::new(1.0, 0.0), Vec2::new(0.0, -1.0), Vec2::new(0.6, 0.8), Vec2::new(-1.0, 1.0)]);
+    let side = Vec2::new(-dir.y, dir.x);
+    let mut at = Point::new(r.pick(&[0.0, 3.0, -7.5]), r.pick(&[0.0, 2.0]));
+    let anchors = (0..n)
+        .map(|_| {
+            at += dir * r.pick(&[1.0, 1.0, 1.0, 2.0, 0.5, 0.0, -1.0]);
+            let bend = r.pick(&[0.0, 0.0, 0.0, 0.0, 1e-3, 0.004, 0.005, 0.006, 1.0]);
+            let mut p = at + side * bend;
+            match r.next() % 40 {
+                0 => p.x = f64::NAN,
+                1 => p.y = f64::INFINITY,
+                2 => p.x = f64::NEG_INFINITY,
+                3 => p = Point::new(p.x * 1e300, p.y * 1e300),
+                4 => p = Point::new(1e-300, -1e-300),
+                _ => {}
+            }
+            let handle = |r: &mut Rng| p + Vec2::new(r.pick(&[0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1e-10, 1.0]), 0.0);
+            let (h_in, h_out) = (handle(r), handle(r));
+            Anchor { p, h_in, h_out, kind: r.pick(&[AnchorKind::Corner, AnchorKind::Corner, AnchorKind::Smooth]) }
+        })
+        .collect();
+    SubPath::new(anchors, closed)
+}
+
+const TOLERANCES: [f64; 12] = [-1.0, -0.0, 0.0, 1e-9, 0.001, 0.005, 0.5, 1.0, 2.0, f64::NAN, f64::INFINITY, f64::NEG_INFINITY];
+
+#[test]
+fn collinear_anchors_go_as_when_removed_one_at_a_time() {
+    let mut r = Rng(0x9e37_79b9_7f4a_7c15);
+    for n in (0..40).chain([64, 100, 257]) {
+        for _ in 0..40 {
+            for closed in [false, true] {
+                let sp = generated(&mut r, n, closed);
+                for tolerance in TOLERANCES {
+                    same_as_one_at_a_time(&sp, tolerance);
+                }
+            }
+        }
+    }
+}
+
+#[test]
+fn collinear_anchors_keep_their_neighbors_order_and_ends() {
+    let ring = |pts: &[(f64, f64)]| line(pts, true);
+    let cases = [
+        // A closed run whose first anchor goes: the last anchor's next neighbor is then the second.
+        ring(&[(5.0, 0.0), (10.0, 0.0), (10.0, 10.0), (0.0, 10.0), (0.0, 0.0)]),
+        ring(&[(0.0, 0.0), (1.0, 0.0), (2.0, 0.0), (3.0, 0.0), (4.0, 0.0), (5.0, 0.0), (2.0, 4.0)]),
+        // Every anchor on one line: the path stops at three (closed) or two (open) anchors.
+        ring(&[(0.0, 0.0), (1.0, 0.0), (2.0, 0.0), (3.0, 0.0), (4.0, 0.0)]),
+        line(&[(0.0, 0.0), (1.0, 0.0), (2.0, 0.0), (3.0, 0.0), (4.0, 0.0)], false),
+        // The last anchor lies between the first one kept and the one before it.
+        ring(&[(0.0, 0.0), (10.0, 0.0), (10.0, 10.0), (0.0, 10.0), (0.0, 5.0)]),
+        // Repeated and zero-length points.
+        line(&[(0.0, 0.0), (0.0, 0.0), (1.0, 0.0), (1.0, 0.0), (2.0, 0.0)], false),
+        ring(&[(0.0, 0.0), (0.0, 0.0), (0.0, 0.0), (0.0, 0.0)]),
+        // Products of huge coordinates overflow; tiny ones underflow.
+        line(&[(0.0, 0.0), (1e200, 1e200), (2e200, 2e200), (f64::MAX, f64::MAX)], false),
+        line(&[(0.0, 0.0), (1e-200, 0.0), (2e-200, 0.0)], false),
+        line(&[(0.0, 0.0), (f64::NAN, 0.0), (2.0, 0.0), (f64::INFINITY, 0.0), (5.0, 0.0)], false),
+    ];
+    for sp in &cases {
+        for tolerance in TOLERANCES {
+            same_as_one_at_a_time(sp, tolerance);
+        }
+    }
+    let mut first = cases[0].clone();
+    assert_eq!(remove_collinear_anchors(&mut first, 0.0), 1);
+    assert_eq!(points(&first), [(10.0, 0.0), (10.0, 10.0), (0.0, 10.0), (0.0, 0.0)]);
+    let (mut closed, mut open) = (cases[2].clone(), cases[3].clone());
+    assert_eq!((remove_collinear_anchors(&mut closed, 0.0), closed.anchors.len()), (2, 3));
+    assert_eq!(remove_collinear_anchors(&mut open, 0.0), 3);
+    assert_eq!(points(&open), [(0.0, 0.0), (4.0, 0.0)], "open ends stay");
+    // A handle on either side of either segment keeps the anchor between them.
+    for side in 0..4 {
+        let mut sp = line(&[(0.0, 0.0), (1.0, 0.0), (2.0, 0.0)], false);
+        let (i, h) = [(0, 1), (1, 0), (1, 1), (2, 0)][side];
+        for kind in [AnchorKind::Corner, AnchorKind::Smooth] {
+            let a = &mut sp.anchors[i];
+            a.kind = kind;
+            *(if h == 0 { &mut a.h_in } else { &mut a.h_out }) = a.p + Vec2::new(0.0, 1.0);
+            same_as_one_at_a_time(&sp, 0.005);
+            assert_eq!(remove_collinear_anchors(&mut sp.clone(), 0.005), 0, "handle {h} on anchor {i}");
+        }
+    }
+}
+
+#[test]
+fn collinear_anchors_at_the_tolerance_go_and_just_past_it_stay() {
+    let sp = line(&[(0.0, 0.0), (100.0, 0.1), (200.0, 0.0)], false);
+    let (a, b) = (sp.anchors[1].p - sp.anchors[0].p, sp.anchors[2].p - sp.anchors[1].p);
+    let sine = a.cross(b).abs() / (a.hypot() * b.hypot());
+    for (tolerance, removed) in [(sine, 1), (sine.next_down(), 0), (sine.next_up(), 1)] {
+        same_as_one_at_a_time(&sp, tolerance);
+        assert_eq!(remove_collinear_anchors(&mut sp.clone(), tolerance), removed, "{tolerance}");
+    }
+}
+
+#[test]
+fn collinear_anchors_go_in_linear_work() {
+    for (n, closed) in [(10_000, false), (20_000, false), (20_000, true)] {
+        let pts: Vec<Point> = (0..n).map(|i| Point::new(i as f64, 0.0)).collect();
+        let mut sp = SubPath::polyline(&pts, closed);
+        let removed = remove_collinear_anchors(&mut sp, 0.005);
+        // A closed straight run stops at three anchors; an open one keeps its two ends.
+        assert_eq!((removed, sp.anchors.len()), if closed { (n - 3, 3) } else { (n - 2, 2) });
+        let work = super::COLLINEAR_WORK.with(std::cell::Cell::get);
+        assert!(work <= 2 * n, "{n} anchors took {work} steps");
+    }
+}
+
 #[test]
 fn duplicate_anchors_merge_keeping_the_outer_handles() {
     let mut sp = SubPath::new(
@@ -89,6 +262,244 @@ fn near_points_group_but_linked_and_equal_ones_do_not() {
     assert!(near_point_groups(&far, 1e-300, |_, _| false).is_empty());
 }
 
+/// Nearly Coincident Anchors' grouping comparing every pair of points in touching grid cells: the
+/// groups [`near_point_groups`] gives.
+fn near_groups_pair_by_pair(points: &[Point], distance: f64, linked: impl Fn(usize, usize) -> bool) -> Vec<Vec<usize>> {
+    if !(distance > 0.0 && distance.is_finite()) {
+        return vec![];
+    }
+    let cell = |p: Point| ((p.x / distance).floor() as i64, (p.y / distance).floor() as i64);
+    let mut grid: HashMap<(i64, i64), Vec<usize>> = HashMap::new();
+    for (i, p) in points.iter().enumerate().filter(|(_, p)| p.is_finite()) {
+        grid.entry(cell(*p)).or_default().push(i);
+    }
+    let mut parent: Vec<usize> = (0..points.len()).collect();
+    for (&(cx, cy), members) in &grid {
+        for (dx, dy) in [(0, 0), (1, 0), (-1, 1), (0, 1), (1, 1)] {
+            let Some(others) = grid.get(&(cx.saturating_add(dx), cy.saturating_add(dy))) else { continue };
+            for &i in members {
+                for &j in others.iter().filter(|&&j| (dx, dy) != (0, 0) || j > i) {
+                    let d = points[i].distance(points[j]);
+                    if d <= distance && d > 0.0 && !linked(i.min(j), i.max(j)) {
+                        let (ri, rj) = (root(&mut parent, i), root(&mut parent, j));
+                        parent[ri.max(rj)] = ri.min(rj);
+                    }
+                }
+            }
+        }
+    }
+    let mut groups: HashMap<usize, Vec<usize>> = HashMap::new();
+    for i in 0..points.len() {
+        groups.entry(root(&mut parent, i)).or_default().push(i);
+    }
+    let mut out: Vec<Vec<usize>> = groups.into_values().filter(|g| g.len() > 1).collect();
+    out.sort_unstable();
+    out
+}
+
+/// Points around a few spots on a lattice, many at the same place, with signed zeros, huge, tiny
+/// and non-finite coordinates.
+fn generated_points(r: &mut Rng, n: usize, distance: f64) -> Vec<Point> {
+    let spots: Vec<Point> = (0..1 + r.next() % 6).map(|_| Point::new((r.next() % 5) as f64 * distance, (r.next() % 5) as f64 * distance)).collect();
+    (0..n)
+        .map(|_| {
+            let jitter = |r: &mut Rng| r.pick(&[0.0, 0.0, -0.0, 0.1, 0.3, 0.5, 0.7, 0.99, 1.0, 1.01, 2.0, -0.5, -1.0]) * distance;
+            let at = r.pick(&spots);
+            let p = Point::new(at.x + jitter(r), at.y + jitter(r));
+            match r.next() % 50 {
+                0 => Point::new(f64::NAN, p.y),
+                1 => Point::new(p.x, f64::INFINITY),
+                2 => Point::new(-0.0, -0.0),
+                3 => Point::new(p.x * 1e300, p.y),
+                4 => Point::new(1e-300, 1e-300),
+                _ => p,
+            }
+        })
+        .collect()
+}
+
+#[test]
+fn near_points_group_as_when_every_pair_is_compared() {
+    let mut r = Rng(0x2545_f491_4f6c_dd1d);
+    let linked: [fn(usize, usize) -> bool; 5] =
+        [|_, _| false, |a, b| b == a + 1, |a, b| (a * 31 + b * 17) % 5 == 0, |a, b| (a + b) % 2 == 0, |_, _| true];
+    for n in (0..30).chain([64, 150, 300]) {
+        for _ in 0..12 {
+            for distance in [0.01, 0.5, 1.0, 3.0, 1e-200, 1e300] {
+                let pts = generated_points(&mut r, n, distance);
+                for (k, linked) in linked.iter().enumerate() {
+                    let want = near_groups_pair_by_pair(&pts, distance, linked);
+                    assert_eq!(near_point_groups(&pts, distance, linked), want, "linked {k} at {distance} for {pts:?}");
+                }
+            }
+        }
+    }
+}
+
+/// A tight cluster of `n` points around `at`, no two at the same place.
+fn cluster(at: Point, radius: f64, n: usize) -> Vec<Point> {
+    (0..n).map(|i| at + Vec2::from_angle(i as f64 * 2.399_963) * (radius * ((i + 1) as f64 / n as f64).sqrt())).collect()
+}
+
+#[test]
+fn near_points_in_a_dense_cell_group_in_linear_work() {
+    let n = 4_000;
+    let neighbours = |a: usize, b: usize| b == a + 1;
+    let one_cluster = cluster(Point::new(0.005, 0.005), 0.004, n);
+    let mut two_clusters = cluster(Point::new(0.001, 0.001), 0.0009, n / 2);
+    two_clusters.extend(cluster(Point::new(0.019, 0.019), 0.0009, n / 2));
+    let same_place = vec![Point::new(0.005, 0.005); n];
+    let cases = [(&one_cluster, 1, n), (&two_clusters, 2, n / 2), (&same_place, 0, 0)];
+    for (pts, count, size) in cases {
+        let groups = near_point_groups(pts, 0.01, neighbours);
+        let work = super::NEAR_WORK.with(std::cell::Cell::get);
+        assert_eq!(groups, near_groups_pair_by_pair(pts, 0.01, neighbours));
+        assert_eq!(groups.len(), count);
+        assert!(groups.iter().all(|g| g.len() == size), "{:?}", groups.iter().map(Vec::len).collect::<Vec<_>>());
+        assert!(work <= 4 * n, "{n} points took {work} steps");
+    }
+}
+
+/// One pair of ends (`a` < `b` in the list of ends) for [`rejoin_sorting_every_pair`]: the turn,
+/// the distance, and the two ends.
+type Pair = (f64, f64, End, End);
+
+/// Rejoin listing every pair of ends in touching grid cells and sorting them by turn and then
+/// distance: the joins [`rejoin`] makes.
+fn rejoin_sorting_every_pair(pieces: &[SubPath], gap: f64) -> Vec<Rejoined> {
+    if !(gap >= 0.0 && gap.is_finite()) {
+        return vec![];
+    }
+    let usable = |sp: &&SubPath| {
+        !sp.closed && sp.anchors.len() >= 2 && sp.anchors.iter().all(|a| a.p.is_finite() && a.h_in.is_finite() && a.h_out.is_finite())
+    };
+    let ends: Vec<(End, Point, Vec2)> = pieces
+        .iter()
+        .enumerate()
+        .filter(|(_, sp)| usable(sp))
+        .flat_map(|(i, sp)| [(i, sp, false), (i, sp, true)])
+        .filter_map(|(i, sp, last)| {
+            let a = if last { sp.anchors.last() } else { sp.anchors.first() }?;
+            Some(((i, last), a.p, outward(sp, last)?))
+        })
+        .collect();
+    let size = gap.max(1e-9);
+    let cell = |p: Point| ((p.x / size).floor() as i64, (p.y / size).floor() as i64);
+    let mut grid: HashMap<(i64, i64), Vec<usize>> = HashMap::new();
+    for (i, e) in ends.iter().enumerate() {
+        grid.entry(cell(e.1)).or_default().push(i);
+    }
+    let mut pairs: Vec<Pair> = vec![];
+    for (a, &(ea, pa, oa)) in ends.iter().enumerate() {
+        let (cx, cy) = cell(pa);
+        for dx in -1..=1 {
+            for dy in -1..=1 {
+                let Some(others) = grid.get(&(cx.saturating_add(dx), cy.saturating_add(dy))) else { continue };
+                for &b in others.iter().filter(|&&b| b > a) {
+                    let (eb, pb, ob) = ends[b];
+                    if ea.0 == eb.0 && pieces[ea.0].anchors.len() < 3 {
+                        continue;
+                    }
+                    let d = pa.distance(pb);
+                    if d <= gap {
+                        pairs.push((oa.dot(-ob).clamp(-1.0, 1.0).acos(), d, ea, eb));
+                    }
+                }
+            }
+        }
+    }
+    pairs.sort_by(|x, y| x.0.total_cmp(&y.0).then(x.1.total_cmp(&y.1)));
+    let mut link: HashMap<End, End> = HashMap::new();
+    for (_, _, a, b) in pairs {
+        if !link.contains_key(&a) && !link.contains_key(&b) {
+            link.insert(a, b);
+            link.insert(b, a);
+        }
+    }
+    chains(pieces, &link)
+}
+
+/// Check that [`rejoin`] joins `pieces` as [`rejoin_sorting_every_pair`] does, bit for bit.
+fn same_as_sorting_every_pair(pieces: &[SubPath], gap: f64) {
+    let (got, want) = (rejoin(pieces, gap), rejoin_sorting_every_pair(pieces, gap));
+    let shape = |out: &[Rejoined]| out.iter().map(|j| (j.members.clone(), bits(&j.path))).collect::<Vec<_>>();
+    assert_eq!(shape(&got), shape(&want), "gap {gap} for {pieces:?}");
+}
+
+/// Pieces of two to four anchors whose ends sit on or near a few lattice points, often at the
+/// same place, some curving, closed or non-finite.
+fn generated_pieces(r: &mut Rng, n: usize, gap: f64) -> Vec<SubPath> {
+    let spot = |r: &mut Rng| {
+        let near = r.pick(&[0.0, 0.0, 0.0, 0.25, 0.5, 1.0, 1.5]) * gap;
+        Point::new((r.next() % 3) as f64 + near, (r.next() % 3) as f64 + r.pick(&[0.0, 0.0, near, -near]))
+    };
+    (0..n)
+        .map(|_| {
+            let count = 2 + (r.next() % 3) as usize;
+            let mut pts: Vec<Point> = (0..count).map(|_| spot(r)).collect();
+            if r.next().is_multiple_of(30) {
+                pts[0].x = f64::NAN;
+            }
+            let mut sp = SubPath::polyline(&pts, r.next().is_multiple_of(15));
+            if r.next().is_multiple_of(4) {
+                let a = &mut sp.anchors[0];
+                a.h_out = a.p + r.pick(&[Vec2::new(0.3, 0.1), Vec2::new(0.0, -0.2), Vec2::new(-0.1, -0.1)]);
+            }
+            sp
+        })
+        .collect()
+}
+
+#[test]
+fn rejoin_joins_as_when_every_pair_is_sorted() {
+    let mut r = Rng(0xd1b5_4a32_d192_ed03);
+    for n in (0..25).chain([40, 80, 160]) {
+        for _ in 0..15 {
+            for gap in [0.0, 0.05, 0.1, 0.5, 1.0, 3.0] {
+                same_as_sorting_every_pair(&generated_pieces(&mut r, n, gap), gap);
+            }
+        }
+    }
+    // Ends meeting from directions that tie: a plus sign, a doubled star and stacked copies.
+    let plus: Vec<SubPath> = [(1.0, 0.0), (0.0, 1.0), (-1.0, 0.0), (0.0, -1.0)].map(|(x, y)| line(&[(0.0, 0.0), (x, y)], false)).to_vec();
+    let mut doubled = plus.clone();
+    doubled.extend(plus.iter().map(|sp| {
+        let mut sp = sp.clone();
+        sp.reverse();
+        sp
+    }));
+    let stacked = vec![line(&[(0.0, 0.0), (1.0, 0.0), (2.0, 1.0)], false); 7];
+    for pieces in [&plus, &doubled, &stacked] {
+        for gap in [0.0, 0.1, 5.0] {
+            same_as_sorting_every_pair(pieces, gap);
+        }
+    }
+}
+
+#[test]
+fn rejoin_pairs_coincident_ends_in_near_linear_work() {
+    let n = 2_000;
+    let ray = |angle: f64| {
+        let tip = Vec2::from_angle(angle) * 100.0;
+        line(&[(0.0, 0.0), (tip.x, tip.y)], false)
+    };
+    let star: Vec<SubPath> = (0..n).map(|i| ray(i as f64 * std::f64::consts::TAU / n as f64)).collect();
+    let mut r = Rng(0x9e37_79b9_7f4a_7c15);
+    let scattered: Vec<SubPath> = (0..n).map(|_| ray((r.next() % 1_000_000) as f64 * 1e-6 * std::f64::consts::TAU)).collect();
+    let stacked = vec![line(&[(0.0, 0.0), (10.0, 0.0)], false); n];
+    for (name, pieces) in [("star", &star), ("scattered", &scattered), ("stacked", &stacked)] {
+        let out = rejoin(pieces, 0.1);
+        let work = super::REJOIN_WORK.with(std::cell::Cell::get);
+        let want = rejoin_sorting_every_pair(pieces, 0.1);
+        assert_eq!(
+            out.iter().map(|j| (j.members.clone(), bits(&j.path))).collect::<Vec<_>>(),
+            want.iter().map(|j| (j.members.clone(), bits(&j.path))).collect::<Vec<_>>(),
+            "{name}"
+        );
+        assert!(work <= 20 * n, "{name}: {n} pieces took {work} steps");
+    }
+}
+
 #[test]
 fn rejoin_welds_pieces_within_the_gap() {
     // Three pieces of one zig-zag, the middle one drawn backwards, plus one too far away.