diff --git a/Cargo.toml b/Cargo.toml
index 28c2eb7..964bbeb 100644
--- a/Cargo.toml
+++ b/Cargo.toml
@@ -63,6 +63,10 @@ name = "synthetic"
 harness = false
 required-features = ["generators"]
 
+[[bench]]
+name = "painted_dreams"
+harness = false
+
 [[example]]
 name = "boolean_op"
 required-features = ["generators"]
diff --git a/benches/painted_dreams.rs b/benches/painted_dreams.rs
new file mode 100644
index 0000000..aa96cc1
--- /dev/null
+++ b/benches/painted_dreams.rs
@@ -0,0 +1,18 @@
+use criterion::{criterion_group, criterion_main, Criterion};
+use kurbo::BezPath;
+
+use linesweeper::{binary_op_with_eps, BinaryOp, FillRule};
+
+fn painted_dreams(c: &mut Criterion) {
+    let path_a = BezPath::from_svg("M0,340C161.737914,383.575765 107.564182,490.730587 273,476 C419,463 481.741198,514.692273 481.333333,768 C481.333333,768 -0,768 -0,768 C-0,768 0,340 0,340 Z").unwrap();
+    let path_b = BezPath::from_svg(
+		"M458.370270,572.165771C428.525848,486.720093 368.618805,467.485992 273,476 C107.564178,490.730591 161.737915,383.575775 0,340 C0,340 0,689 0,689 C56,700 106.513901,779.342590 188,694.666687 C306.607422,571.416260 372.033966,552.205139 458.370270,572.165771 Z",
+	).unwrap();
+
+    c.bench_function("painted dreams", |b| {
+        b.iter(|| binary_op_with_eps(&path_a, &path_b, FillRule::EvenOdd, BinaryOp::Union, 1e-5))
+    });
+}
+
+criterion_group!(benches, painted_dreams);
+criterion_main!(benches);
diff --git a/profile.json.gz b/profile.json.gz
new file mode 100644
index 0000000..a183754
Binary files /dev/null and b/profile.json.gz differ
diff --git a/src/curve/mod.rs b/src/curve/mod.rs
index d8a2d64..4c8b2cb 100644
--- a/src/curve/mod.rs
+++ b/src/curve/mod.rs
@@ -1,7 +1,7 @@
 //! Curve comparison utilities.
 use arrayvec::ArrayVec;
 use kurbo::{
-    common::solve_cubic, Affine, CubicBez, Line, ParamCurve, PathSeg, QuadBez, Shape, Vec2,
+    common::solve_cubic, Affine, CubicBez, Line, ParamCurve, PathSeg, Point, QuadBez, Shape, Vec2,
 };
 use polycool::Cubic;
 
@@ -967,6 +967,17 @@ fn with_rotation(
     true
 }
 
+fn l_infty_distance(p0: Point, p1: Point) -> f64 {
+    (p0.x - p1.x).abs().max((p0.y - p1.y).abs())
+}
+
+fn max_distance(c0: CubicBez, c1: CubicBez) -> f64 {
+    l_infty_distance(c0.p0, c1.p0)
+        .max(l_infty_distance(c0.p1, c1.p1))
+        .max(l_infty_distance(c0.p2, c1.p2))
+        .max(l_infty_distance(c0.p3, c1.p3))
+}
+
 fn intersect_cubics_rec(
     orig_c0: CubicBez,
     orig_c1: CubicBez,
@@ -982,6 +993,10 @@ fn intersect_cubics_rec(
     // dbg!(c0, orig_c0);
     // dbg!(c1, orig_c1);
 
+    if max_distance(c0, c1) <= tolerance {
+        out.push(y1, Order::Ish);
+        return;
+    }
     if y1 - y0 < accuracy {
         // For very short intervals there's some numerical instability in constructing the
         // approximating quadratics, so we just do a coarser comparison based on bounding
@@ -1238,8 +1253,6 @@ pub mod arbtests {
 
 #[cfg(test)]
 mod test {
-    use crate::Segment;
-
     use super::*;
 
     #[test]
@@ -1431,56 +1444,26 @@ mod test {
         );
     }
 
+    // This test recurses a bit too deeply. The two curves are quite close,
+    // and near the end they're just a bit more than eps apart.
     #[test]
-    fn graphite_example() {
-        // Fiddling with graphite turned up this example, in which c0 and c1 compare "ish",
-        // c1 compares left of c2, and c0 compares right of c2. That's all fine: c1, c2, c0
-        // would be a valid sweep-line order.
-        //
-        // But suppose that we start with c0, c1 and then try to insert c2. Since it's to
-        // the right of c1 even with the bigger thresholds, we think it's ok to put to the
-        // right of c1, and then the bigger threshold comparison stops us from seeing c0
-        // when scanning to the left.
-        //
-        // I think the solution has to be to re-introduce the close_before and close_after
-        // stuff for y-slop.
+    fn painted_dreams_deep_recursion() {
         let eps = 1e-6;
         let tolerance = eps;
         let accuracy = eps / 2.0;
-        let y = -227.53699416;
         let c0 = CubicBez::new(
-            (-4.04445106, -227.53699448),
-            (-4.0443963, -227.53699414000002),
-            (-4.0443415400000005, -227.5369938),
-            (-4.04428678, -227.53699347),
+            (268.2700181987094, 465.9600610772311),
+            (261.5036543284458, 474.04917281766745),
+            (255.01588390924064, 481.16474559246694),
+            (248.76144328725377, 487.4191862144538),
         );
         let c1 = CubicBez::new(
-            (-4.04445106, -227.53699448),
-            (-4.04445141, -227.53699414000002),
-            (-4.04445176, -227.5369938),
-            (-4.0444521, -227.53699347),
+            (268.2700181987094, 465.9600610772311),
+            (261.508835610275, 474.04297901254995),
+            (255.02582000674386, 481.1538483487484),
+            (248.77581162941254, 487.404817872295),
         );
-        let c2 = CubicBez::new(
-            (-4.04442515, -227.53699416),
-            (-4.04443617, -227.53699411),
-            (-4.0444375500000005, -227.53699405),
-            (-4.04442929, -227.536994),
-        );
-        let cmp01 = intersect_cubics(c0, c1, tolerance, accuracy).with_y_slop(tolerance);
-        let cmp02 = intersect_cubics(c0, c2, tolerance, accuracy).with_y_slop(tolerance);
-        let cmp12 = intersect_cubics(c1, c2, tolerance, accuracy).with_y_slop(tolerance);
-
-        dbg!(&cmp01, &cmp12, &cmp02);
-        dbg!(solve_x_for_y(c0, y));
-        dbg!(solve_x_for_y(c1, y));
-        assert!(c2.p0.x >= dbg!(c1.p0.x.max(c1.p1.x).max(c1.p2.x).max(c1.p3.x)) + 2.0 * eps);
-
-        let c0_seg = Segment::from_kurbo_cubic(c0);
-        let c1_seg = Segment::from_kurbo_cubic(c1);
-        let c2_seg = Segment::from_kurbo_cubic(c2);
-
-        dbg!(c0_seg.lower(y, eps));
-        dbg!(c1_seg.lower(y, eps));
-        dbg!(c2_seg.lower(y, eps));
+
+        dbg!(intersect_cubics(c0, c1, tolerance, accuracy));
     }
 }
diff --git a/src/lib.rs b/src/lib.rs
index 8009977..b82a941 100644
--- a/src/lib.rs
+++ b/src/lib.rs
@@ -84,7 +84,12 @@ impl std::fmt::Display for Error {
 
 impl std::error::Error for Error {}
 
-/// Computes a boolean operation between two sets, each of which is described as a collection of closed polylines.
+/// Computes a boolean operation between two sets, each of which is described as
+/// a collection of closed polylines.
+///
+/// This function calls [`binary_op_with_eps`] with an automatically-chosen
+/// accuracy parameter. The automatically-chosen accuracy is fairly
+/// conservative, so you may get better performance by choosing your own.
 pub fn binary_op(
     set_a: &kurbo::BezPath,
     set_b: &kurbo::BezPath,
@@ -110,7 +115,21 @@ pub fn binary_op(
     let eps = eps.max(1e-6);
 
     debug_assert!(eps.is_finite());
+    binary_op_with_eps(set_a, set_b, fill_rule, op, eps)
+}
 
+/// Computes a boolean operation between two sets, each of which is described as
+/// a collection of closed polylines.
+///
+/// The output paths will not exactly agree with the input paths in general.
+/// The accuracy parameter `eps` bounds how much they can differ.
+pub fn binary_op_with_eps(
+    set_a: &kurbo::BezPath,
+    set_b: &kurbo::BezPath,
+    fill_rule: FillRule,
+    op: BinaryOp,
+    eps: f64,
+) -> Result<topology::Contours, Error> {
     let top = Topology::from_paths_binary(set_a, set_b, eps)?;
     #[cfg(feature = "debug-svg")]
     {