Fix #3437: Spiral segments follow the spiral

The Spiral tool places its anchors on a logarithmic spiral, with the
radius shrinking by the decay each quarter turn, but gave every handle
the direction of a circle's tangent and a length from an averaged
radius. The curve bulged off the spiral between anchors and turned at
each anchor, so the spiral looked squared: up to 2.8% off its radius at
the default 80% decay and 23% at 30%.

The handles now follow the spiral's own tangent, with lengths that
differ at each segment's outer and inner end. The curve stays within
0.03% of the spiral at 80% decay and 0.8% at 20%. Anchor positions and
counts are unchanged.

crates/geom/src/shapes.rs

@@ -210,21 +210,27 @@ pub fn arc(a: Point, b: Point, slope: f64, closed: bool, _base_along_x: bool) ->
     PathData::single(SubPath::new(anchors, closed))
 }
 
-/// Archimedean-style spiral approximating Illustrator's Spiral tool (decay-based, `segments` quarter turns).
+/// Logarithmic spiral approximating the reference app's Spiral tool (decay-based, `segments` quarter turns).
 pub fn spiral(center: Point, radius: f64, decay_percent: f64, segments: u32, clockwise: bool) -> PathData {
     let decay = (decay_percent / 100.0).clamp(0.05, 0.9999);
     let dir = if clockwise { 1.0 } else { -1.0 };
+    // The radius shrinks by `decay` each quarter turn: r(φ) = radius · e^(kφ) after turning φ.
+    let k = decay.ln() / (PI / 2.0);
+    // Handle lengths per unit of the derivative below, at a segment's outer and inner end: √r
+    // changes evenly with the curve's parameter along each quarter turn. Both are KAPPA when the
+    // radius does not shrink, and the curve stays within 0.8% of the spiral's radius at 20% decay.
+    let outer = 2.0 * KAPPA * (decay.sqrt() - 1.0) / decay.ln();
+    let inner = outer / decay.sqrt();
     let mut anchors = Vec::new();
     let mut r = radius;
     for i in 0..=segments {
         let a = dir * (i as f64) * PI / 2.0;
         let p = Point::new(center.x + r * a.cos(), center.y + r * a.sin());
-        // Tangent direction for a quarter arc of radius r.
-        let t = Vec2::new(-a.sin(), a.cos()) * dir;
-        let r_prev = r / decay;
-        let r_next = r * decay;
-        let h_in = p - t * (KAPPA * (r + r_prev) / 2.0);
-        let h_out = p + t * (KAPPA * (r + r_next) / 2.0);
+        // The spiral's derivative by the angle turned: the radius changes along the radial
+        // direction while the point turns along the circle's tangent.
+        let d = Vec2::new(a.cos(), a.sin()) * (r * k) + Vec2::new(-a.sin(), a.cos()) * (r * dir);
+        let h_in = p - d * inner;
+        let h_out = p + d * outer;
         anchors.push(Anchor { p, h_in: if i == 0 { p } else { h_in }, h_out: if i == segments { p } else { h_out }, kind: AnchorKind::Smooth });
         r *= decay;
     }
@@ -460,6 +466,34 @@ mod tests {
         assert!(!p.subpaths[0].closed);
     }
 
+    /// #3437: every segment stays on the spiral through its anchors, whose radius shrinks by
+    /// `decay` each quarter turn, for strong and weak decay and both directions.
+    #[test]
+    fn spiral_segments_follow_the_spiral() {
+        let center = Point::new(40.0, -25.0);
+        for decay in [20.0, 30.0, 50.0, 80.0, 95.0] {
+            for clockwise in [true, false] {
+                let p = spiral(center, 100.0, decay, 8, clockwise);
+                let anchors = &p.subpaths[0].anchors;
+                let k = (decay / 100.0_f64).ln() / (PI / 2.0);
+                let mut worst = 0.0_f64;
+                for pair in anchors.windows(2) {
+                    let c = kurbo::CubicBez::new(pair[0].p, pair[0].h_out, pair[1].h_in, pair[1].p);
+                    let r0 = (pair[0].p - center).hypot();
+                    let a0 = (pair[0].p - center).atan2();
+                    for step in 1..20 {
+                        let q = kurbo::ParamCurve::eval(&c, f64::from(step) / 20.0) - center;
+                        // The angle turned from the segment's start, toward the outer end.
+                        let turned = (q.atan2() - a0 + 3.0 * PI).rem_euclid(TAU) - PI;
+                        let want = r0 * (-k * turned.abs()).exp();
+                        worst = worst.max((q.hypot() - want).abs() / want);
+                    }
+                }
+                assert!(worst < 0.008, "decay {decay} clockwise {clockwise}: radius off by {:.2}%", worst * 100.0);
+            }
+        }
+    }
+
     #[test]
     fn grids() {
         let r = Rect::new(0.0, 0.0, 10.0, 10.0);