diff --git a/crates/color/src/gradient.rs b/crates/color/src/gradient.rs
index 402cf3ae2..8441ebcc5 100644
--- a/crates/color/src/gradient.rs
+++ b/crates/color/src/gradient.rs
@@ -175,6 +175,40 @@ fn between(a: &GradientStop, b: &GradientStop, t: f32, mix: &impl Fn(&Color, &Co
     (mix(&a.color, &b.color, v), a.opacity + (b.opacity - a.opacity) * v)
 }
 
+/// For each pair of neighbouring stops (`k`, `k + 1`), the largest offset among stops `1..=k + 1`.
+/// [`Gradient::sample_with`] picks the first pair whose end offset `t` doesn't pass; that is the
+/// first pair whose reach `t` doesn't pass, and the reaches never decrease, so the pair is found by
+/// binary search whether or not the stops are in order.
+struct StopReach(Vec<f32>);
+
+impl StopReach {
+    fn new(stops: &[GradientStop]) -> Self {
+        // `f32::max` skips a NaN offset, which `t <= offset` never chooses either.
+        let reach = stops.iter().skip(1).scan(f32::NEG_INFINITY, |m, s| {
+            *m = m.max(s.offset);
+            Some(*m)
+        });
+        Self(reach.collect())
+    }
+
+    /// What [`Gradient::sample_with`] gives at `t` for `stops`, the stops this was built from.
+    fn sample_with(&self, stops: &[GradientStop], t: f32, mix: &impl Fn(&Color, &Color, f32) -> Color) -> (Color, f32) {
+        let Some(first) = stops.first() else { return (Color::BLACK, 1.0) };
+        if t <= first.offset {
+            return (first.color, first.opacity);
+        }
+        let last = || stops.last().map_or((Color::BLACK, 1.0), |l| (l.color, l.opacity));
+        if t.is_nan() {
+            return last();
+        }
+        let k = self.0.partition_point(|r| *r < t);
+        match (stops.get(k), stops.get(k + 1)) {
+            (Some(a), Some(b)) => between(a, b, t, mix),
+            _ => last(),
+        }
+    }
+}
+
 fn one() -> f32 {
     1.0
 }
@@ -245,9 +279,11 @@ impl Gradient {
     /// the one before it stopped, so the work is the samples plus the stops, not their product.
     pub fn sample_increasing_with(&self, ts: impl IntoIterator<Item = f32>, mix: impl Fn(&Color, &Color, f32) -> Color) -> Vec<(Color, f32)> {
         let stops = &self.stops;
-        // The shared index assumes ordered stops; others (a file's, unsorted) sample one by one.
+        // The shared cursor assumes ordered stops; others (a file's, unsorted) look each sample up
+        // in a [`StopReach`], so they cost a binary search each rather than a scan of the stops.
         if stops.windows(2).any(|w| w[1].offset < w[0].offset) {
-            return ts.into_iter().map(|t| self.sample_with(t, &mix)).collect();
+            let reach = StopReach::new(stops);
+            return ts.into_iter().map(|t| reach.sample_with(stops, t, &mix)).collect();
         }
         let mut j = 0;
         ts.into_iter()
@@ -284,6 +320,9 @@ impl Gradient {
     fn expanded_with_stops(&self) -> impl Iterator<Item = (f32, Color, f32, Option<f32>, Option<usize>)> + '_ {
         // Direct interpolation matches `sample` only for strictly increasing offsets.
         let increasing = self.stops.windows(2).all(|w| w[0].offset < w[1].offset);
+        // Otherwise (repeated or unsorted offsets) each midpoint is looked up as `sample` would
+        // find it, but by binary search: a scan of all the stops per midpoint cost their square.
+        let reach = if increasing || !self.method.is_classic() { None } else { Some(StopReach::new(&self.stops)) };
         self.stops.iter().enumerate().flat_map(move |(i, s)| {
             let next = self.stops.get(i + 1);
             let off_centre = (s.midpoint - 0.5).abs() > 1e-3;
@@ -293,7 +332,10 @@ impl Gradient {
                 // Mixed from the two stops it lies between, as `sample` would: searching all the
                 // stops again for each midpoint made many stops cost their square.
                 let mix = |a: &Color, b: &Color, v| self.method.mix(a, b, v);
-                let (c, o) = if increasing { between(s, n, t, &mix) } else { self.sample_with(t, mix) };
+                let (c, o) = match &reach {
+                    Some(reach) => reach.sample_with(&self.stops, t, &mix),
+                    None => between(s, n, t, &mix),
+                };
                 (t, c, o, Some(s.midpoint), None)
             });
             // Perceptual: stops along the mix (`v`) where a straight line would stray from it, the
diff --git a/crates/color/src/gradient/tests.rs b/crates/color/src/gradient/tests.rs
index 224803075..da9a3f2cc 100644
--- a/crates/color/src/gradient/tests.rs
+++ b/crates/color/src/gradient/tests.rs
@@ -23,6 +23,49 @@ fn increasing_samples_match_single_samples() {
     assert_eq!(g.sample_increasing_with(ts, Color::lerp_keep_model), one_by_one);
 }
 
+/// Looking stops up by their reach finds what scanning them does, for any order of offsets,
+/// repeated, out-of-order and NaN ones included, and for positions outside and between them.
+#[test]
+fn stop_reach_matches_scanning_the_stops() {
+    let colors = ["#ff0000", "#00ff00", "#0000ff", "#ffff00", "#00ffff", "#000000"];
+    let offsets: [&[f32]; 5] = [
+        &[0.1, 0.4, 0.4, 0.2, 0.9, 0.5],
+        &[0.5, 0.5, 0.5, 0.5],
+        &[0.9, 0.1, 0.8, 0.2, 0.7, 0.3],
+        &[0.2, f32::NAN, 0.6, 0.1, f32::NAN, 1.0],
+        &[f32::NAN, 0.3, 0.3, 0.8],
+    ];
+    for offsets in offsets {
+        let stops: Vec<GradientStop> = offsets
+            .iter()
+            .zip(colors.iter().cycle())
+            .map(|(&o, hex)| GradientStop { midpoint: 0.3, ..GradientStop::new(o, Color::from_hex(hex).unwrap()) })
+            .collect();
+        let g = Gradient::new(GradientKind::Linear, stops);
+        let reach = StopReach::new(&g.stops);
+        for t in (-4..=44).map(|k| k as f32 / 40.0).chain([f32::NAN, f32::INFINITY, f32::NEG_INFINITY]) {
+            // Compared as printed: a NaN offset mixes NaN colours, which are equal only so.
+            let (got, want) = (reach.sample_with(&g.stops, t, &Color::lerp_keep_model), g.sample_with(t, Color::lerp_keep_model));
+            assert_eq!(format!("{got:?}"), format!("{want:?}"), "{offsets:?} at {t}");
+        }
+    }
+}
+
+/// Repeated offsets turn off the direct midpoint mix, and the midpoints still come out as
+/// `sample` gives them, with work that stays near linear in the stops (thousands of stops at one
+/// offset used to scan them all for each midpoint).
+#[test]
+fn expanded_with_repeated_offsets_samples_midpoints() {
+    let stop = |offset: f32, hex: &str| GradientStop { midpoint: 0.25, ..GradientStop::new(offset, Color::from_hex(hex).unwrap()) };
+    let g = Gradient::new(GradientKind::Linear, vec![stop(0.0, "#ff0000"), stop(0.5, "#00ff00"), stop(0.5, "#0000ff"), stop(1.0, "#000000")]);
+    for (t, c, o, mid) in g.expanded().filter(|e| e.3.is_some()) {
+        assert_eq!(mid, Some(0.25));
+        assert_eq!((c, o), g.sample(t), "midpoint at {t}");
+    }
+    let many = Gradient::new(GradientKind::Linear, (0..200_000).map(|i| stop(if i == 0 { 0.0 } else { 0.5 }, "#808080")).collect());
+    assert_eq!(many.expanded().count(), 2 * 200_000 - 1);
+}
+
 #[test]
 fn midpoint_shifts() {
     let mut g = Gradient::default();