From e0a54ebf68a1d26eb409fe68e9b1b810fc597945 Mon Sep 17 00:00:00 2001
From: Bruce Mitchener <bruce.mitchener@gmail.com>
Date: Sat, 14 Sep 2024 23:43:34 +0700
Subject: [PATCH 01/30] Impl `ParamCurve`, `ParamCurveArclen` for `Arc`

This uses an approximation of the arc length using beziers as the
analytic solution is quite involved.
---
 CHANGELOG.md |  2 ++
 src/arc.rs   | 40 ++++++++++++++++++++++++++++++++++++++--
 2 files changed, 40 insertions(+), 2 deletions(-)

diff --git a/CHANGELOG.md b/CHANGELOG.md
index f133ff48..128cc582 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -20,6 +20,7 @@ This release has an [MSRV][] of 1.65.
 - `Stroke` is now `PartialEq`, `StrokeOpts` is now `Clone`, `Copy`, `Debug`, `Eq`, `PartialEq`. ([#379][] by [@waywardmonkeys][])
 - Implement `Sum` for `Vec2`. ([#399][] by [@Philipp-M][])
 - Add triangle shape. ([#350][] by [@juliapaci][])
+- `Arc` now implements `ParamCurve` and `ParamCurveArclen`. ([#378] by [@waywardmonkeys])
 
 ### Changed
 
@@ -90,6 +91,7 @@ Note: A changelog was not kept for or before this release
 [#370]: https://github.com/linebender/kurbo/pull/370
 [#375]: https://github.com/linebender/kurbo/pull/375
 [#376]: https://github.com/linebender/kurbo/pull/376
+[#378]: https://github.com/linebender/kurbo/pull/378
 [#379]: https://github.com/linebender/kurbo/pull/379
 [#383]: https://github.com/linebender/kurbo/pull/383
 [#388]: https://github.com/linebender/kurbo/pull/388
diff --git a/src/arc.rs b/src/arc.rs
index 9ac642be..d87ec5ad 100644
--- a/src/arc.rs
+++ b/src/arc.rs
@@ -3,11 +3,11 @@
 
 //! An ellipse arc.
 
-use crate::{Affine, Ellipse, PathEl, Point, Rect, Shape, Vec2};
+use crate::{Affine, Ellipse, ParamCurve, ParamCurveArclen, PathEl, Point, Rect, Shape, Vec2};
 use core::{
     f64::consts::{FRAC_PI_2, PI},
     iter,
-    ops::Mul,
+    ops::{Mul, Range},
 };
 
 #[cfg(not(feature = "std"))]
@@ -171,6 +171,42 @@ fn rotate_pt(pt: Vec2, angle: f64) -> Vec2 {
     )
 }
 
+impl ParamCurve for Arc {
+    fn eval(&self, t: f64) -> Point {
+        let angle = self.start_angle + (self.sweep_angle * t);
+        sample_ellipse(self.radii, self.x_rotation, angle).to_point()
+    }
+
+    fn subsegment(&self, range: Range<f64>) -> Self {
+        Self {
+            center: self.center,
+            radii: self.radii,
+            start_angle: self.start_angle + (self.sweep_angle * range.start),
+            sweep_angle: self.sweep_angle - (self.sweep_angle * (range.end - range.start)),
+            x_rotation: self.x_rotation,
+        }
+    }
+
+    fn start(&self) -> Point {
+        sample_ellipse(self.radii, self.x_rotation, self.start_angle).to_point()
+    }
+
+    fn end(&self) -> Point {
+        sample_ellipse(
+            self.radii,
+            self.x_rotation,
+            self.start_angle + self.sweep_angle,
+        )
+        .to_point()
+    }
+}
+
+impl ParamCurveArclen for Arc {
+    fn arclen(&self, accuracy: f64) -> f64 {
+        self.path_segments(0.1).perimeter(accuracy)
+    }
+}
+
 impl Shape for Arc {
     type PathElementsIter<'iter> = iter::Chain<iter::Once<PathEl>, ArcAppendIter>;
 

From 2d20d64184dacc983ede374591560141d3c85a5e Mon Sep 17 00:00:00 2001
From: Thomas Churchman <thomas@kepow.org>
Date: Fri, 27 Sep 2024 22:19:47 +0200
Subject: [PATCH 02/30] Numerically approximate ParamCurveArclen

---
 src/arc.rs | 220 ++++++++++++++++++++++++++++++++++++++++++++++++++++-
 1 file changed, 218 insertions(+), 2 deletions(-)

diff --git a/src/arc.rs b/src/arc.rs
index d87ec5ad..814be58b 100644
--- a/src/arc.rs
+++ b/src/arc.rs
@@ -202,8 +202,111 @@ impl ParamCurve for Arc {
 }
 
 impl ParamCurveArclen for Arc {
-    fn arclen(&self, accuracy: f64) -> f64 {
-        self.path_segments(0.1).perimeter(accuracy)
+    fn arclen(&self, _accuracy: f64) -> f64 {
+        // TODO: wire up accuracy. The Carlson numerical approximation provides a bound on the relative
+        // error
+        let relative_error = 1e-20;
+
+        /// Numerically approximate the incomplete elliptic integral of the second kind
+        /// parameterized by `phi` and `m = k^2` in Legendre's trigonometric form.
+        fn incomplete_elliptic_integral_second_kind(relative_error: f64, phi: f64, m: f64) -> f64 {
+            // Approximate the incomplete elliptic integral through Carlson symmetric forms:
+            // https://en.wikipedia.org/w/index.php?title=Carlson_symmetric_form&oldid=1223277638#Incomplete_elliptic_integrals
+            phi.sin()
+                * carlson_rf(
+                    relative_error,
+                    phi.cos().powi(2),
+                    1. - m * phi.sin().powi(2),
+                    1.,
+                )
+                - 1. / 3.
+                    * m
+                    * phi.sin().powi(3)
+                    * carlson_rd(
+                        relative_error,
+                        phi.cos().powi(2),
+                        1. - m * phi.sin().powi(2),
+                        1.,
+                    )
+        }
+
+        /// Calculate the length of an arc along an ellipse defined by `radii`, from `start_angle`
+        /// to `end_angle`, with the angles being inside the first quadrant, i.e.,
+        /// 0 <= start_angle <= end_angle <= PI / 2
+        fn quadrant_ellipse_arc_length(
+            relative_error: f64,
+            radii: Vec2,
+            start_angle: f64,
+            end_angle: f64,
+        ) -> f64 {
+            debug_assert!(radii.y >= radii.x);
+            debug_assert!(start_angle >= 0.);
+            debug_assert!(end_angle >= start_angle);
+            debug_assert!(end_angle <= PI / 2.);
+
+            // Ellipse arc length calculated through the incomplete elliptic integral of the second
+            // kind:
+            // https://en.wikipedia.org/w/index.php?title=Ellipse&oldid=1248023575#Arc_length
+            let x0 = sample_ellipse(radii, 0., start_angle).x;
+            let x1 = sample_ellipse(radii, 0., end_angle).x;
+
+            let z0 = (x0 / radii.x).acos();
+            let z1 = (x1 / radii.x).acos();
+
+            let m = 1. - radii.x.powi(2) / radii.y.powi(2);
+
+            radii.y
+                * (incomplete_elliptic_integral_second_kind(relative_error, z1, m)
+                    - incomplete_elliptic_integral_second_kind(relative_error, z0, m))
+        }
+
+        // Normalize ellipse to have radius y >= radius x
+        let (radii, mut start_angle) = if self.radii.y > self.radii.x {
+            (self.radii, self.start_angle)
+        } else {
+            (
+                Vec2::new(self.radii.y, self.radii.x),
+                self.start_angle + PI / 2.,
+            )
+        };
+
+        let mut sweep_angle = self.sweep_angle;
+        // Normalize sweep angle to be non-negative
+        if sweep_angle < 0. {
+            start_angle -= sweep_angle;
+            sweep_angle = -sweep_angle;
+        }
+
+        // Normalize start_angle to be non-negative and to be on the first quadrant of the ellipse
+        start_angle = start_angle.rem_euclid(PI / 2.);
+
+        if start_angle + sweep_angle > PI / 2. {
+            // The arc crosses from the first into the second quadrant, first calculate from start
+            // to end of first quadrant
+            let mut arclen = 0.;
+            arclen += quadrant_ellipse_arc_length(relative_error, radii, start_angle, PI / 2.);
+
+            sweep_angle -= PI / 2. - start_angle;
+            if sweep_angle >= PI / 2. {
+                // The sweep angle may twist around ellipse quadrants multiple times
+                let quarter_turns = (sweep_angle / (PI / 2.)).floor();
+                // Calculating quarter of a complete ellipse arc length, this could be
+                // special-cased as it's a complete elliptic integral
+                arclen +=
+                    quarter_turns * quadrant_ellipse_arc_length(relative_error, radii, 0., PI / 2.)
+            }
+
+            sweep_angle = sweep_angle.rem_euclid(PI / 2.);
+            arclen += quadrant_ellipse_arc_length(relative_error, radii, 0., sweep_angle);
+            arclen
+        } else {
+            quadrant_ellipse_arc_length(
+                relative_error,
+                radii,
+                start_angle,
+                start_angle + sweep_angle,
+            )
+        }
     }
 }
 
@@ -259,6 +362,92 @@ impl Mul<Arc> for Affine {
     }
 }
 
+/// Approximation of Carlson RF function as defined in "Numerical computation of real or complex
+/// elliptic integrals" (Carlson, Bille C.): https://arxiv.org/abs/math/9409227v1
+fn carlson_rf(relative_error: f64, x: f64, y: f64, z: f64) -> f64 {
+    let mut x = x;
+    let mut y = y;
+    let mut z = z;
+
+    let a0 = (x + y + z) / 3.;
+    let q = (3. * relative_error).powf(-1. / 6.)
+        * (a0 - x).abs().max((a0 - y).abs()).max((a0 - z).abs());
+
+    let mut a = a0;
+    let mut m = 0;
+    loop {
+        if 4f64.powi(-m) * q <= a.abs() {
+            break;
+        }
+
+        let lambda = x.sqrt() * y.sqrt() + x.sqrt() * z.sqrt() + y.sqrt() * z.sqrt();
+
+        a = (a + lambda) / 4.;
+        x = (x + lambda) / 4.;
+        y = (y + lambda) / 4.;
+        z = (z + lambda) / 4.;
+
+        m += 1;
+    }
+
+    let x = (a0 - x) / 4f64.powi(m) * a;
+    let y = (a0 - y) / 4f64.powi(m) * a;
+    let z = -x - y;
+
+    let e2 = x * y - z.powi(2);
+    let e3 = x * y * z;
+
+    a.powf(-1. / 2.)
+        * (1. - 1. / 10. * e2 + 1. / 14. * e3 + 1. / 24. * e2.powi(2) - 3. / 44. * e2 * e3)
+}
+
+/// Approximation of Carlson RD function as defined in "Numerical computation of real or complex
+/// elliptic integrals" (Carlson, Bille C.): https://arxiv.org/abs/math/9409227v1
+fn carlson_rd(relative_error: f64, x: f64, y: f64, z: f64) -> f64 {
+    let mut x = x;
+    let mut y = y;
+    let mut z = z;
+
+    let a0 = (x + y + 3. * z) / 5.;
+    let q = (relative_error / 4.).powf(-1. / 6.)
+        * (a0 - x).abs().max((a0 - y).abs()).max((a0 - z).abs());
+
+    let mut sum = 0.;
+    let mut a = a0;
+    let mut m = 0;
+    loop {
+        if 4f64.powi(-m) * q <= a.abs() {
+            break;
+        }
+
+        let lambda = x.sqrt() * y.sqrt() + x.sqrt() * z.sqrt() + y.sqrt() * z.sqrt();
+        sum += 4f64.powi(-m) / (z.sqrt() * (z + lambda));
+        a = (a + lambda) / 4.;
+        x = (x + lambda) / 4.;
+        y = (y + lambda) / 4.;
+        z = (z + lambda) / 4.;
+
+        m += 1;
+    }
+
+    let x = (a0 - x) / (4f64.powi(4) * a);
+    let y = (a0 - y) / (4f64.powi(4) * a);
+    let z = -(x + y) / 3.;
+
+    let e2 = x * y - 6. * z.powi(2);
+    let e3 = (3. * x * y - 8. * z.powi(2)) * z;
+    let e4 = 3. * x * y - z.powi(2) * z.powi(2);
+    let e5 = x * y * z.powi(3);
+
+    4f64.powi(-m)
+        * a.powf(-3. / 2.)
+        * (1. - 3. / 14. * e2 + 1. / 6. * e3 + 9. / 88. * e2.powi(2)
+            - 3. / 22. * e4
+            - 9. / 52. * e2 * e3
+            + 3. / 26. * e5)
+        + 3. * sum
+}
+
 #[cfg(test)]
 mod tests {
     use super::*;
@@ -278,4 +467,31 @@ mod tests {
         // Reversing it again should result in the original arc
         assert_eq!(a, f.reversed());
     }
+
+    #[test]
+    fn length() {
+        // TODO: when arclen actually uses specified accuracy, update EPSILON and the accuracy
+        // params
+        const EPSILON: f64 = 1e-6;
+
+        let a = Arc::new((0., 0.), (1., 1.), 0., PI * 4., 0.);
+        assert!((dbg!(a.arclen(0.1)) - PI * 4.).abs() <= EPSILON);
+
+        // TODO: add some more known cases (perhaps calculated by other implementations)
+    }
+
+    #[test]
+    fn carlson_numerical_checks() {
+        // TODO: relative bound on error doesn't seem to be quite correct yet, use a large epsilon
+        // for now
+        const EPSILON: f64 = 1e-6;
+
+        // Numerical checks from section 3 of "Numerical computation of real or complex elliptic
+        // integrals" (Carlson, Bille C.): https://arxiv.org/abs/math/9409227v1 (real-valued calls)
+        assert!((carlson_rf(1e-20, 1., 2., 0.) - 1.3110_28777_1461).abs() <= EPSILON);
+        assert!((carlson_rf(1e-20, 2., 3., 4.) - 0.58408_28416_7715).abs() <= EPSILON);
+
+        assert!((carlson_rd(1e-20, 0., 2., 1.) - 1.7972_10352_1034).abs() <= EPSILON);
+        assert!((carlson_rd(1e-20, 2., 3., 4.) - 0.16510_52729_4261).abs() <= EPSILON);
+    }
 }

From 7bd3bf5d552e05bc444a233a1146e79a18a3c724 Mon Sep 17 00:00:00 2001
From: Thomas Churchman <thomas@kepow.org>
Date: Sat, 28 Sep 2024 12:18:08 +0200
Subject: [PATCH 03/30] Calculate for ellipse upper half, simplify arccos(a
 cos(x)/a)

Also add some more tests
---
 src/arc.rs | 187 ++++++++++++++++++++++++++++++-----------------------
 1 file changed, 105 insertions(+), 82 deletions(-)

diff --git a/src/arc.rs b/src/arc.rs
index 814be58b..002ff860 100644
--- a/src/arc.rs
+++ b/src/arc.rs
@@ -207,61 +207,8 @@ impl ParamCurveArclen for Arc {
         // error
         let relative_error = 1e-20;
 
-        /// Numerically approximate the incomplete elliptic integral of the second kind
-        /// parameterized by `phi` and `m = k^2` in Legendre's trigonometric form.
-        fn incomplete_elliptic_integral_second_kind(relative_error: f64, phi: f64, m: f64) -> f64 {
-            // Approximate the incomplete elliptic integral through Carlson symmetric forms:
-            // https://en.wikipedia.org/w/index.php?title=Carlson_symmetric_form&oldid=1223277638#Incomplete_elliptic_integrals
-            phi.sin()
-                * carlson_rf(
-                    relative_error,
-                    phi.cos().powi(2),
-                    1. - m * phi.sin().powi(2),
-                    1.,
-                )
-                - 1. / 3.
-                    * m
-                    * phi.sin().powi(3)
-                    * carlson_rd(
-                        relative_error,
-                        phi.cos().powi(2),
-                        1. - m * phi.sin().powi(2),
-                        1.,
-                    )
-        }
-
-        /// Calculate the length of an arc along an ellipse defined by `radii`, from `start_angle`
-        /// to `end_angle`, with the angles being inside the first quadrant, i.e.,
-        /// 0 <= start_angle <= end_angle <= PI / 2
-        fn quadrant_ellipse_arc_length(
-            relative_error: f64,
-            radii: Vec2,
-            start_angle: f64,
-            end_angle: f64,
-        ) -> f64 {
-            debug_assert!(radii.y >= radii.x);
-            debug_assert!(start_angle >= 0.);
-            debug_assert!(end_angle >= start_angle);
-            debug_assert!(end_angle <= PI / 2.);
-
-            // Ellipse arc length calculated through the incomplete elliptic integral of the second
-            // kind:
-            // https://en.wikipedia.org/w/index.php?title=Ellipse&oldid=1248023575#Arc_length
-            let x0 = sample_ellipse(radii, 0., start_angle).x;
-            let x1 = sample_ellipse(radii, 0., end_angle).x;
-
-            let z0 = (x0 / radii.x).acos();
-            let z1 = (x1 / radii.x).acos();
-
-            let m = 1. - radii.x.powi(2) / radii.y.powi(2);
-
-            radii.y
-                * (incomplete_elliptic_integral_second_kind(relative_error, z1, m)
-                    - incomplete_elliptic_integral_second_kind(relative_error, z0, m))
-        }
-
-        // Normalize ellipse to have radius y >= radius x
-        let (radii, mut start_angle) = if self.radii.y > self.radii.x {
+        // Normalize ellipse to have radius x >= radius y
+        let (radii, mut start_angle) = if self.radii.x >= self.radii.y {
             (self.radii, self.start_angle)
         } else {
             (
@@ -270,43 +217,36 @@ impl ParamCurveArclen for Arc {
             )
         };
 
-        let mut sweep_angle = self.sweep_angle;
         // Normalize sweep angle to be non-negative
+        let mut sweep_angle = self.sweep_angle;
         if sweep_angle < 0. {
             start_angle -= sweep_angle;
             sweep_angle = -sweep_angle;
         }
 
-        // Normalize start_angle to be non-negative and to be on the first quadrant of the ellipse
-        start_angle = start_angle.rem_euclid(PI / 2.);
-
-        if start_angle + sweep_angle > PI / 2. {
-            // The arc crosses from the first into the second quadrant, first calculate from start
-            // to end of first quadrant
-            let mut arclen = 0.;
-            arclen += quadrant_ellipse_arc_length(relative_error, radii, start_angle, PI / 2.);
-
-            sweep_angle -= PI / 2. - start_angle;
-            if sweep_angle >= PI / 2. {
-                // The sweep angle may twist around ellipse quadrants multiple times
-                let quarter_turns = (sweep_angle / (PI / 2.)).floor();
-                // Calculating quarter of a complete ellipse arc length, this could be
-                // special-cased as it's a complete elliptic integral
-                arclen +=
-                    quarter_turns * quadrant_ellipse_arc_length(relative_error, radii, 0., PI / 2.)
-            }
-
-            sweep_angle = sweep_angle.rem_euclid(PI / 2.);
-            arclen += quadrant_ellipse_arc_length(relative_error, radii, 0., sweep_angle);
-            arclen
+        start_angle = start_angle.rem_euclid(PI);
+
+        let mut arclen = 0.;
+        let half_turns = (sweep_angle / PI).floor();
+        if half_turns > 0. {
+            // Half of the ellipse circumference is a complete elliptic integral and could be special-cased
+            arclen += half_turns * half_ellipse_arc_length(relative_error, radii, 0., PI);
+            sweep_angle = sweep_angle.rem_euclid(PI);
+        }
+
+        if start_angle + sweep_angle > PI {
+            arclen += half_ellipse_arc_length(relative_error, radii, start_angle, PI);
+            arclen += half_ellipse_arc_length(relative_error, radii, 0., PI - sweep_angle);
         } else {
-            quadrant_ellipse_arc_length(
+            arclen += half_ellipse_arc_length(
                 relative_error,
                 radii,
                 start_angle,
                 start_angle + sweep_angle,
-            )
+            );
         }
+
+        arclen
     }
 }
 
@@ -448,6 +388,66 @@ fn carlson_rd(relative_error: f64, x: f64, y: f64, z: f64) -> f64 {
         + 3. * sum
 }
 
+/// Numerically approximate the incomplete elliptic integral of the second kind
+/// parameterized by `phi` and `m = k^2` in Legendre's trigonometric form.
+fn incomplete_elliptic_integral_second_kind(relative_error: f64, phi: f64, m: f64) -> f64 {
+    // Approximate the incomplete elliptic integral through Carlson symmetric forms:
+    // https://en.wikipedia.org/w/index.php?title=Carlson_symmetric_form&oldid=1223277638#Incomplete_elliptic_integrals
+
+    debug_assert!(phi >= -PI / 2.);
+    debug_assert!(phi <= PI / 2.);
+    debug_assert!(m * phi.sin().powi(2) >= 0.);
+    debug_assert!(m * phi.sin().powi(2) <= 1.);
+
+    phi.sin()
+        * carlson_rf(
+            relative_error,
+            phi.cos().powi(2),
+            1. - m * phi.sin().powi(2),
+            1.,
+        )
+        - 1. / 3.
+            * m
+            * phi.sin().powi(3)
+            * carlson_rd(
+                relative_error,
+                phi.cos().powi(2),
+                1. - m * phi.sin().powi(2),
+                1.,
+            )
+}
+
+/// Calculate the length of an arc along an ellipse defined by `radii`, from `start_angle`
+/// to `end_angle`, with the angles being inside the first two quadrants, i.e.,
+/// 0 <= start_angle <= end_angle <= PI.
+///
+/// This assumes radii.x >= radii.y
+fn half_ellipse_arc_length(
+    relative_error: f64,
+    radii: Vec2,
+    start_angle: f64,
+    end_angle: f64,
+) -> f64 {
+    debug_assert!(radii.x >= radii.y);
+    debug_assert!(start_angle >= 0.);
+    debug_assert!(end_angle >= start_angle);
+    debug_assert!(end_angle <= PI);
+
+    let radii = Vec2::new(radii.y, radii.x);
+    let start_angle = start_angle - PI / 2.;
+    let end_angle = end_angle - PI / 2.;
+
+    // Ellipse arc length calculated through the incomplete elliptic integral of the second
+    // kind:
+    // https://en.wikipedia.org/w/index.php?title=Ellipse&oldid=1248023575#Arc_length
+
+    let m = 1. - (radii.x / radii.y).powi(2);
+
+    radii.y
+        * (incomplete_elliptic_integral_second_kind(relative_error, end_angle, m)
+            - incomplete_elliptic_integral_second_kind(relative_error, start_angle, m))
+}
+
 #[cfg(test)]
 mod tests {
     use super::*;
@@ -475,9 +475,13 @@ mod tests {
         const EPSILON: f64 = 1e-6;
 
         let a = Arc::new((0., 0.), (1., 1.), 0., PI * 4., 0.);
-        assert!((dbg!(a.arclen(0.1)) - PI * 4.).abs() <= EPSILON);
+        assert!((a.arclen(0.000_1) - PI * 4.).abs() <= EPSILON);
+
+        let a = Arc::new((0., 0.), (2.23, 3.05), 0., 0.2, 0.);
+        assert!((a.arclen(0.000_1) - 0.60811714277).abs() <= EPSILON);
 
-        // TODO: add some more known cases (perhaps calculated by other implementations)
+        let a = Arc::new((0., 0.), (3.05, 2.23), 0., 0.2, 0.);
+        assert!((a.arclen(0.000_1) - 0.448555).abs() <= EPSILON);
     }
 
     #[test]
@@ -494,4 +498,23 @@ mod tests {
         assert!((carlson_rd(1e-20, 0., 2., 1.) - 1.7972_10352_1034).abs() <= EPSILON);
         assert!((carlson_rd(1e-20, 2., 3., 4.) - 0.16510_52729_4261).abs() <= EPSILON);
     }
+
+    #[test]
+    fn elliptic_e_numerical_checks() {
+        const EPSILON: f64 = 1e-6;
+
+        for (phi, m, elliptic_e) in [
+            (0.0, 0.0, 0.0),
+            (0.5, 0.0, 0.5),
+            (1.0, 0.0, 1.0),
+            (0.0, 1.0, 0.0),
+            (1.0, 1.0, 0.84147098),
+        ] {
+            let elliptic_e_approx = incomplete_elliptic_integral_second_kind(1e-20, phi, m);
+            assert!(
+                (elliptic_e_approx - elliptic_e).abs() < EPSILON,
+                "Approximated elliptic e {elliptic_e_approx} does not match known value {elliptic_e} for E({phi}|{m})"
+            );
+        }
+    }
 }

From 2f8606e9fb3c0d0fc20b39ed846f3c1f0d2f2290 Mon Sep 17 00:00:00 2001
From: Thomas Churchman <thomas@kepow.org>
Date: Sat, 28 Sep 2024 13:15:57 +0200
Subject: [PATCH 04/30] Correctly wrap around

---
 src/arc.rs | 19 ++++++++++++++++++-
 1 file changed, 18 insertions(+), 1 deletion(-)

diff --git a/src/arc.rs b/src/arc.rs
index 002ff860..60ffd6e7 100644
--- a/src/arc.rs
+++ b/src/arc.rs
@@ -236,7 +236,8 @@ impl ParamCurveArclen for Arc {
 
         if start_angle + sweep_angle > PI {
             arclen += half_ellipse_arc_length(relative_error, radii, start_angle, PI);
-            arclen += half_ellipse_arc_length(relative_error, radii, 0., PI - sweep_angle);
+            arclen +=
+                half_ellipse_arc_length(relative_error, radii, 0., start_angle + sweep_angle - PI);
         } else {
             arclen += half_ellipse_arc_length(
                 relative_error,
@@ -474,6 +475,22 @@ mod tests {
         // params
         const EPSILON: f64 = 1e-6;
 
+        // Circular checks:
+        for (start_angle, sweep_angle, length) in [
+            (0., 1., 1.),
+            (0., 2., 2.),
+            (0., 5., 5.),
+            (1.0, 3., 3.),
+            (1.5, 10., 10.),
+        ] {
+            let a = Arc::new((0., 0.), (1., 1.), start_angle, sweep_angle, 0.);
+            let arc_length = a.arclen(0.000_1);
+            assert!(
+                (arc_length - length).abs() <= EPSILON,
+                "Got arc length {arc_length}, expected {length} for circular arc {a:?}"
+            );
+        }
+
         let a = Arc::new((0., 0.), (1., 1.), 0., PI * 4., 0.);
         assert!((a.arclen(0.000_1) - PI * 4.).abs() <= EPSILON);
 

From 6537992d880a303a38e0422f74b0b735a9e7a8de Mon Sep 17 00:00:00 2001
From: Thomas Churchman <thomas@kepow.org>
Date: Sat, 28 Sep 2024 13:18:35 +0200
Subject: [PATCH 05/30] Fix negative sweep angle normalization

---
 src/arc.rs | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/arc.rs b/src/arc.rs
index 60ffd6e7..f35987ae 100644
--- a/src/arc.rs
+++ b/src/arc.rs
@@ -220,7 +220,7 @@ impl ParamCurveArclen for Arc {
         // Normalize sweep angle to be non-negative
         let mut sweep_angle = self.sweep_angle;
         if sweep_angle < 0. {
-            start_angle -= sweep_angle;
+            start_angle = PI - start_angle;
             sweep_angle = -sweep_angle;
         }
 

From 50250662f6e4aa7e1a0e2a7a3dad6e8bcef60f72 Mon Sep 17 00:00:00 2001
From: Thomas Churchman <thomas@kepow.org>
Date: Sat, 28 Sep 2024 13:20:47 +0200
Subject: [PATCH 06/30] Add check comparing with bezier perimeter result

---
 src/arc.rs | 21 +++++++++++++++++++++
 1 file changed, 21 insertions(+)

diff --git a/src/arc.rs b/src/arc.rs
index f35987ae..b87d2c20 100644
--- a/src/arc.rs
+++ b/src/arc.rs
@@ -501,6 +501,27 @@ mod tests {
         assert!((a.arclen(0.000_1) - 0.448555).abs() <= EPSILON);
     }
 
+    #[test]
+    fn length_compare_with_bez_length() {
+        const EPSILON: f64 = 1e-3;
+
+        for radii in [(1., 1.), (0.5, 1.), (2., 1.)] {
+            for start_angle in [0., 0.5, 1., 2., PI, -1.] {
+                for sweep_angle in [0., 0.5, 1., 2., PI, -1.] {
+                    let a = Arc::new((0., 0.), radii, start_angle, sweep_angle, 0.);
+
+                    let arc_length = a.arclen(0.000_1);
+                    let bez_length = a.path_segments(0.000_1).perimeter(0.000_1);
+
+                    assert!(
+                         (arc_length - bez_length).abs() < EPSILON,
+                         "Numerically approximated arc length ({arc_length}) does not match bezier segment perimeter length ({bez_length}) for arc {a:?}"
+                     );
+                }
+            }
+        }
+    }
+
     #[test]
     fn carlson_numerical_checks() {
         // TODO: relative bound on error doesn't seem to be quite correct yet, use a large epsilon

From 65b548481a97330c2d314cfda17cf728e1f1f5d2 Mon Sep 17 00:00:00 2001
From: Thomas Churchman <thomas@kepow.org>
Date: Sat, 28 Sep 2024 13:37:59 +0200
Subject: [PATCH 07/30] Improve docs

---
 src/arc.rs | 12 ++++++++----
 1 file changed, 8 insertions(+), 4 deletions(-)

diff --git a/src/arc.rs b/src/arc.rs
index b87d2c20..53d24766 100644
--- a/src/arc.rs
+++ b/src/arc.rs
@@ -418,8 +418,8 @@ fn incomplete_elliptic_integral_second_kind(relative_error: f64, phi: f64, m: f6
             )
 }
 
-/// Calculate the length of an arc along an ellipse defined by `radii`, from `start_angle`
-/// to `end_angle`, with the angles being inside the first two quadrants, i.e.,
+/// Calculate the length of an arc along an ellipse defined by `radii`, from `start_angle` to
+/// `end_angle`, with the angles being inside the upper half of the ellipse, i.e.,
 /// 0 <= start_angle <= end_angle <= PI.
 ///
 /// This assumes radii.x >= radii.y
@@ -434,12 +434,16 @@ fn half_ellipse_arc_length(
     debug_assert!(end_angle >= start_angle);
     debug_assert!(end_angle <= PI);
 
+    // This function takes angles from 0 to pi for convenience of the caller, but the numerical
+    // methods used here operate from -1/2 pi to 1/2 pi. Rotate the ellipse by 1/2 pi radians (a
+    // quarter turn):
     let radii = Vec2::new(radii.y, radii.x);
     let start_angle = start_angle - PI / 2.;
     let end_angle = end_angle - PI / 2.;
 
-    // Ellipse arc length calculated through the incomplete elliptic integral of the second
-    // kind:
+    // The arc length is equal to radii.y * (E(end_angle | m) - E(start_angle | m)) with E the
+    // incomplete elliptic integral of the second kind and parameter
+    // m = 1 - (radii.x / radii.y)^2 = k^2. See also:
     // https://en.wikipedia.org/w/index.php?title=Ellipse&oldid=1248023575#Arc_length
 
     let m = 1. - (radii.x / radii.y).powi(2);

From 361783e38f16c4c40cf2674075c115230313d8fe Mon Sep 17 00:00:00 2001
From: Thomas Churchman <thomas@kepow.org>
Date: Sat, 28 Sep 2024 13:40:02 +0200
Subject: [PATCH 08/30] Improve readability

---
 src/arc.rs | 23 +++++++----------------
 1 file changed, 7 insertions(+), 16 deletions(-)

diff --git a/src/arc.rs b/src/arc.rs
index 53d24766..58bbb73a 100644
--- a/src/arc.rs
+++ b/src/arc.rs
@@ -400,22 +400,13 @@ fn incomplete_elliptic_integral_second_kind(relative_error: f64, phi: f64, m: f6
     debug_assert!(m * phi.sin().powi(2) >= 0.);
     debug_assert!(m * phi.sin().powi(2) <= 1.);
 
-    phi.sin()
-        * carlson_rf(
-            relative_error,
-            phi.cos().powi(2),
-            1. - m * phi.sin().powi(2),
-            1.,
-        )
-        - 1. / 3.
-            * m
-            * phi.sin().powi(3)
-            * carlson_rd(
-                relative_error,
-                phi.cos().powi(2),
-                1. - m * phi.sin().powi(2),
-                1.,
-            )
+    let (sin, cos) = phi.sin_cos();
+    let sin2 = sin.powi(2);
+    let sin3 = sin.powi(3);
+    let cos2 = cos.powi(2);
+
+    sin * carlson_rf(relative_error, cos2, 1. - m * sin2, 1.)
+        - 1. / 3. * m * sin3 * carlson_rd(relative_error, cos2, 1. - m * sin2, 1.)
 }
 
 /// Calculate the length of an arc along an ellipse defined by `radii`, from `start_angle` to

From 8ff7424c92cc6adbd19910ca192334a0db0d3048 Mon Sep 17 00:00:00 2001
From: Thomas Churchman <thomas@kepow.org>
Date: Sat, 28 Sep 2024 13:42:21 +0200
Subject: [PATCH 09/30] Clippy

---
 src/arc.rs | 8 ++++----
 1 file changed, 4 insertions(+), 4 deletions(-)

diff --git a/src/arc.rs b/src/arc.rs
index 58bbb73a..cffbde48 100644
--- a/src/arc.rs
+++ b/src/arc.rs
@@ -304,7 +304,7 @@ impl Mul<Arc> for Affine {
 }
 
 /// Approximation of Carlson RF function as defined in "Numerical computation of real or complex
-/// elliptic integrals" (Carlson, Bille C.): https://arxiv.org/abs/math/9409227v1
+/// elliptic integrals" (Carlson, Bille C.): <https://arxiv.org/abs/math/9409227v1>
 fn carlson_rf(relative_error: f64, x: f64, y: f64, z: f64) -> f64 {
     let mut x = x;
     let mut y = y;
@@ -343,7 +343,7 @@ fn carlson_rf(relative_error: f64, x: f64, y: f64, z: f64) -> f64 {
 }
 
 /// Approximation of Carlson RD function as defined in "Numerical computation of real or complex
-/// elliptic integrals" (Carlson, Bille C.): https://arxiv.org/abs/math/9409227v1
+/// elliptic integrals" (Carlson, Bille C.): <https://arxiv.org/abs/math/9409227v1>
 fn carlson_rd(relative_error: f64, x: f64, y: f64, z: f64) -> f64 {
     let mut x = x;
     let mut y = y;
@@ -411,9 +411,9 @@ fn incomplete_elliptic_integral_second_kind(relative_error: f64, phi: f64, m: f6
 
 /// Calculate the length of an arc along an ellipse defined by `radii`, from `start_angle` to
 /// `end_angle`, with the angles being inside the upper half of the ellipse, i.e.,
-/// 0 <= start_angle <= end_angle <= PI.
+/// 0 <= `start_angle` <= `end_angle` <= PI.
 ///
-/// This assumes radii.x >= radii.y
+/// This assumes `radii.x` >= `radii.y`
 fn half_ellipse_arc_length(
     relative_error: f64,
     radii: Vec2,

From ea71123f43f04c08a484843c1f36194864ec4301 Mon Sep 17 00:00:00 2001
From: Thomas Churchman <thomas@kepow.org>
Date: Sat, 28 Sep 2024 13:59:08 +0200
Subject: [PATCH 10/30] Improve docs

---
 src/arc.rs | 10 +++++++---
 1 file changed, 7 insertions(+), 3 deletions(-)

diff --git a/src/arc.rs b/src/arc.rs
index cffbde48..708524c9 100644
--- a/src/arc.rs
+++ b/src/arc.rs
@@ -303,8 +303,10 @@ impl Mul<Arc> for Affine {
     }
 }
 
-/// Approximation of Carlson RF function as defined in "Numerical computation of real or complex
+/// Approximation of the Carlson RF function as defined in "Numerical computation of real or complex
 /// elliptic integrals" (Carlson, Bille C.): <https://arxiv.org/abs/math/9409227v1>
+///
+/// RF = 1/2 ∫ 1 / ( sqrt(t+x) sqrt(t+y) sqrt(t+z) ) dt from 0 to inf
 fn carlson_rf(relative_error: f64, x: f64, y: f64, z: f64) -> f64 {
     let mut x = x;
     let mut y = y;
@@ -342,8 +344,10 @@ fn carlson_rf(relative_error: f64, x: f64, y: f64, z: f64) -> f64 {
         * (1. - 1. / 10. * e2 + 1. / 14. * e3 + 1. / 24. * e2.powi(2) - 3. / 44. * e2 * e3)
 }
 
-/// Approximation of Carlson RD function as defined in "Numerical computation of real or complex
-/// elliptic integrals" (Carlson, Bille C.): <https://arxiv.org/abs/math/9409227v1>
+/// Approximation of the Carlson RD function as defined in "Numerical computation of real or
+/// complex elliptic integrals" (Carlson, Bille C.): <https://arxiv.org/abs/math/9409227v1>
+///
+/// RD = 3/2 ∫ 1 / ( sqrt(t+x) sqrt(t+y) (t+z)^(3/2) ) dt from 0 to inf
 fn carlson_rd(relative_error: f64, x: f64, y: f64, z: f64) -> f64 {
     let mut x = x;
     let mut y = y;

From 1cdabf4143fc11e1099f49f915f5f025b1147777 Mon Sep 17 00:00:00 2001
From: Thomas Churchman <thomas@kepow.org>
Date: Sat, 28 Sep 2024 14:31:35 +0200
Subject: [PATCH 11/30] Probable performance improvements

---
 src/arc.rs | 5 ++---
 1 file changed, 2 insertions(+), 3 deletions(-)

diff --git a/src/arc.rs b/src/arc.rs
index 708524c9..37e5ffd2 100644
--- a/src/arc.rs
+++ b/src/arc.rs
@@ -340,7 +340,7 @@ fn carlson_rf(relative_error: f64, x: f64, y: f64, z: f64) -> f64 {
     let e2 = x * y - z.powi(2);
     let e3 = x * y * z;
 
-    a.powf(-1. / 2.)
+    1. / a.sqrt()
         * (1. - 1. / 10. * e2 + 1. / 14. * e3 + 1. / 24. * e2.powi(2) - 3. / 44. * e2 * e3)
 }
 
@@ -384,8 +384,7 @@ fn carlson_rd(relative_error: f64, x: f64, y: f64, z: f64) -> f64 {
     let e4 = 3. * x * y - z.powi(2) * z.powi(2);
     let e5 = x * y * z.powi(3);
 
-    4f64.powi(-m)
-        * a.powf(-3. / 2.)
+    4f64.powi(-m) * 1. / (a * a.sqrt())
         * (1. - 3. / 14. * e2 + 1. / 6. * e3 + 9. / 88. * e2.powi(2)
             - 3. / 22. * e4
             - 9. / 52. * e2 * e3

From f91abce9f90b8878be5f3decd9267d90dc6eced2 Mon Sep 17 00:00:00 2001
From: Thomas Churchman <thomas@kepow.org>
Date: Sat, 28 Sep 2024 14:44:09 +0200
Subject: [PATCH 12/30] Fix wording

---
 src/arc.rs | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/src/arc.rs b/src/arc.rs
index 37e5ffd2..cf782242 100644
--- a/src/arc.rs
+++ b/src/arc.rs
@@ -392,8 +392,8 @@ fn carlson_rd(relative_error: f64, x: f64, y: f64, z: f64) -> f64 {
         + 3. * sum
 }
 
-/// Numerically approximate the incomplete elliptic integral of the second kind
-/// parameterized by `phi` and `m = k^2` in Legendre's trigonometric form.
+/// Numerically approximate the incomplete elliptic integral of the second kind to `phi`
+/// parameterized by `m = k^2` in Legendre's trigonometric form.
 fn incomplete_elliptic_integral_second_kind(relative_error: f64, phi: f64, m: f64) -> f64 {
     // Approximate the incomplete elliptic integral through Carlson symmetric forms:
     // https://en.wikipedia.org/w/index.php?title=Carlson_symmetric_form&oldid=1223277638#Incomplete_elliptic_integrals

From abb02be758bf1bbf59a088fe39db0f0fa15d57d3 Mon Sep 17 00:00:00 2001
From: Thomas Churchman <thomas@kepow.org>
Date: Sun, 29 Sep 2024 17:13:36 +0200
Subject: [PATCH 13/30] Division instead of `powi`

---
 src/arc.rs | 11 ++++++-----
 1 file changed, 6 insertions(+), 5 deletions(-)

diff --git a/src/arc.rs b/src/arc.rs
index cf782242..c60bc9f8 100644
--- a/src/arc.rs
+++ b/src/arc.rs
@@ -313,23 +313,23 @@ fn carlson_rf(relative_error: f64, x: f64, y: f64, z: f64) -> f64 {
     let mut z = z;
 
     let a0 = (x + y + z) / 3.;
-    let q = (3. * relative_error).powf(-1. / 6.)
+    let mut q = (3. * relative_error).powf(-1. / 6.)
         * (a0 - x).abs().max((a0 - y).abs()).max((a0 - z).abs());
 
     let mut a = a0;
     let mut m = 0;
     loop {
-        if 4f64.powi(-m) * q <= a.abs() {
+        if q <= a.abs() {
             break;
         }
 
         let lambda = x.sqrt() * y.sqrt() + x.sqrt() * z.sqrt() + y.sqrt() * z.sqrt();
-
         a = (a + lambda) / 4.;
         x = (x + lambda) / 4.;
         y = (y + lambda) / 4.;
         z = (z + lambda) / 4.;
 
+        q /= 4.;
         m += 1;
     }
 
@@ -354,14 +354,14 @@ fn carlson_rd(relative_error: f64, x: f64, y: f64, z: f64) -> f64 {
     let mut z = z;
 
     let a0 = (x + y + 3. * z) / 5.;
-    let q = (relative_error / 4.).powf(-1. / 6.)
+    let mut q = (relative_error / 4.).powf(-1. / 6.)
         * (a0 - x).abs().max((a0 - y).abs()).max((a0 - z).abs());
 
     let mut sum = 0.;
     let mut a = a0;
     let mut m = 0;
     loop {
-        if 4f64.powi(-m) * q <= a.abs() {
+        if q <= a.abs() {
             break;
         }
 
@@ -372,6 +372,7 @@ fn carlson_rd(relative_error: f64, x: f64, y: f64, z: f64) -> f64 {
         y = (y + lambda) / 4.;
         z = (z + lambda) / 4.;
 
+        q /= 4.;
         m += 1;
     }
 

From f3e12f86f8f768e73787e2c76f0180e37794494f Mon Sep 17 00:00:00 2001
From: Thomas Churchman <thomas@kepow.org>
Date: Mon, 30 Sep 2024 10:17:21 +0200
Subject: [PATCH 14/30] Start angle is modulo half turns

---
 src/arc.rs | 3 ++-
 1 file changed, 2 insertions(+), 1 deletion(-)

diff --git a/src/arc.rs b/src/arc.rs
index c60bc9f8..75a50743 100644
--- a/src/arc.rs
+++ b/src/arc.rs
@@ -220,10 +220,11 @@ impl ParamCurveArclen for Arc {
         // Normalize sweep angle to be non-negative
         let mut sweep_angle = self.sweep_angle;
         if sweep_angle < 0. {
-            start_angle = PI - start_angle;
+            start_angle = -start_angle;
             sweep_angle = -sweep_angle;
         }
 
+        // Normalize start angle to be on the upper half of the ellipse
         start_angle = start_angle.rem_euclid(PI);
 
         let mut arclen = 0.;

From 6ab359209e7b3c78cdad5da8daf8b5c7035ab045 Mon Sep 17 00:00:00 2001
From: Thomas Churchman <thomas@kepow.org>
Date: Mon, 30 Sep 2024 10:33:14 +0200
Subject: [PATCH 15/30] trunc is faster than floor

---
 src/arc.rs | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/arc.rs b/src/arc.rs
index 75a50743..e86e3208 100644
--- a/src/arc.rs
+++ b/src/arc.rs
@@ -228,7 +228,7 @@ impl ParamCurveArclen for Arc {
         start_angle = start_angle.rem_euclid(PI);
 
         let mut arclen = 0.;
-        let half_turns = (sweep_angle / PI).floor();
+        let half_turns = (sweep_angle / PI).trunc();
         if half_turns > 0. {
             // Half of the ellipse circumference is a complete elliptic integral and could be special-cased
             arclen += half_turns * half_ellipse_arc_length(relative_error, radii, 0., PI);

From 3d8fd7f65b00c108bb6db08deb9997c7f89991ca Mon Sep 17 00:00:00 2001
From: Thomas Churchman <thomas@kepow.org>
Date: Mon, 30 Sep 2024 13:57:33 +0200
Subject: [PATCH 16/30] Reduce number of calls to integral calculation

---
 src/arc.rs | 105 ++++++++++++++++++++++++-----------------------------
 1 file changed, 47 insertions(+), 58 deletions(-)

diff --git a/src/arc.rs b/src/arc.rs
index e86e3208..c33d1766 100644
--- a/src/arc.rs
+++ b/src/arc.rs
@@ -207,8 +207,9 @@ impl ParamCurveArclen for Arc {
         // error
         let relative_error = 1e-20;
 
-        // Normalize ellipse to have radius x >= radius y
-        let (radii, mut start_angle) = if self.radii.x >= self.radii.y {
+        // Normalize ellipse to have radius y >= radius x, required for the parameter assumptions
+        // of `incomplete_elliptic_integral_second_kind`.
+        let (radii, mut start_angle) = if self.radii.y >= self.radii.x {
             (self.radii, self.start_angle)
         } else {
             (
@@ -216,6 +217,7 @@ impl ParamCurveArclen for Arc {
                 self.start_angle + PI / 2.,
             )
         };
+        let m = 1. - (radii.x / radii.y).powi(2);
 
         // Normalize sweep angle to be non-negative
         let mut sweep_angle = self.sweep_angle;
@@ -225,30 +227,44 @@ impl ParamCurveArclen for Arc {
         }
 
         // Normalize start angle to be on the upper half of the ellipse
-        start_angle = start_angle.rem_euclid(PI);
-
+        let start_angle = start_angle.rem_euclid(PI);
+
+        let end_angle = start_angle + sweep_angle;
+
+        let mut quarter_turns = (end_angle / PI * 2.).trunc() - (start_angle / PI * 2.).trunc();
+        let end_angle = end_angle % PI;
+
+        // The elliptic arc length is equal to radii.y * (E(end_angle | m) - E(start_angle | m))
+        // with E the incomplete elliptic integral of the second kind and parameter
+        // m = 1 - (radii.x / radii.y)^2 = k^2.
+        //
+        // See also:
+        // https://en.wikipedia.org/w/index.php?title=Ellipse&oldid=1248023575#Arc_length
+        //
+        // The implementation here allows calculating the incomplete elliptic integral in the range
+        // 0 <= phi <= 1/2 pi (the first elliptic quadrant), so split the arc into segments in
+        // that range.
         let mut arclen = 0.;
-        let half_turns = (sweep_angle / PI).trunc();
-        if half_turns > 0. {
-            // Half of the ellipse circumference is a complete elliptic integral and could be special-cased
-            arclen += half_turns * half_ellipse_arc_length(relative_error, radii, 0., PI);
-            sweep_angle = sweep_angle.rem_euclid(PI);
+
+        if start_angle >= PI / 2. {
+            arclen += incomplete_elliptic_integral_second_kind(relative_error, PI - start_angle, m);
+            quarter_turns -= 1.;
+        } else {
+            arclen -= incomplete_elliptic_integral_second_kind(relative_error, start_angle, m);
         }
 
-        if start_angle + sweep_angle > PI {
-            arclen += half_ellipse_arc_length(relative_error, radii, start_angle, PI);
-            arclen +=
-                half_ellipse_arc_length(relative_error, radii, 0., start_angle + sweep_angle - PI);
+        if end_angle >= PI / 2. {
+            arclen -= incomplete_elliptic_integral_second_kind(relative_error, PI - end_angle, m);
+            quarter_turns += 1.;
         } else {
-            arclen += half_ellipse_arc_length(
-                relative_error,
-                radii,
-                start_angle,
-                start_angle + sweep_angle,
-            );
+            arclen += incomplete_elliptic_integral_second_kind(relative_error, end_angle, m);
         }
 
-        arclen
+        // Note: this uses the complete elliptic integral, which can be special-cased.
+        arclen +=
+            quarter_turns * incomplete_elliptic_integral_second_kind(relative_error, PI / 2., m);
+
+        radii.y * arclen
     }
 }
 
@@ -396,6 +412,10 @@ fn carlson_rd(relative_error: f64, x: f64, y: f64, z: f64) -> f64 {
 
 /// Numerically approximate the incomplete elliptic integral of the second kind to `phi`
 /// parameterized by `m = k^2` in Legendre's trigonometric form.
+///
+/// Assumes:
+/// 0 <= phi <= pi / 2
+/// and 0 <= m sin^2(phi) <= 1
 fn incomplete_elliptic_integral_second_kind(relative_error: f64, phi: f64, m: f64) -> f64 {
     // Approximate the incomplete elliptic integral through Carlson symmetric forms:
     // https://en.wikipedia.org/w/index.php?title=Carlson_symmetric_form&oldid=1223277638#Incomplete_elliptic_integrals
@@ -410,45 +430,13 @@ fn incomplete_elliptic_integral_second_kind(relative_error: f64, phi: f64, m: f6
     let sin3 = sin.powi(3);
     let cos2 = cos.powi(2);
 
+    // note: this actually allows calculating from -1/2 pi <= phi <= 1/2 pi, but there are some
+    // alternative translations from the Legendre form that are potentially better, that do
+    // restrict the domain to 0 <= phi <= 1/2 pi.
     sin * carlson_rf(relative_error, cos2, 1. - m * sin2, 1.)
         - 1. / 3. * m * sin3 * carlson_rd(relative_error, cos2, 1. - m * sin2, 1.)
 }
 
-/// Calculate the length of an arc along an ellipse defined by `radii`, from `start_angle` to
-/// `end_angle`, with the angles being inside the upper half of the ellipse, i.e.,
-/// 0 <= `start_angle` <= `end_angle` <= PI.
-///
-/// This assumes `radii.x` >= `radii.y`
-fn half_ellipse_arc_length(
-    relative_error: f64,
-    radii: Vec2,
-    start_angle: f64,
-    end_angle: f64,
-) -> f64 {
-    debug_assert!(radii.x >= radii.y);
-    debug_assert!(start_angle >= 0.);
-    debug_assert!(end_angle >= start_angle);
-    debug_assert!(end_angle <= PI);
-
-    // This function takes angles from 0 to pi for convenience of the caller, but the numerical
-    // methods used here operate from -1/2 pi to 1/2 pi. Rotate the ellipse by 1/2 pi radians (a
-    // quarter turn):
-    let radii = Vec2::new(radii.y, radii.x);
-    let start_angle = start_angle - PI / 2.;
-    let end_angle = end_angle - PI / 2.;
-
-    // The arc length is equal to radii.y * (E(end_angle | m) - E(start_angle | m)) with E the
-    // incomplete elliptic integral of the second kind and parameter
-    // m = 1 - (radii.x / radii.y)^2 = k^2. See also:
-    // https://en.wikipedia.org/w/index.php?title=Ellipse&oldid=1248023575#Arc_length
-
-    let m = 1. - (radii.x / radii.y).powi(2);
-
-    radii.y
-        * (incomplete_elliptic_integral_second_kind(relative_error, end_angle, m)
-            - incomplete_elliptic_integral_second_kind(relative_error, start_angle, m))
-}
-
 #[cfg(test)]
 mod tests {
     use super::*;
@@ -482,6 +470,7 @@ mod tests {
             (0., 5., 5.),
             (1.0, 3., 3.),
             (1.5, 10., 10.),
+            (2.5, 10., 10.),
         ] {
             let a = Arc::new((0., 0.), (1., 1.), start_angle, sweep_angle, 0.);
             let arc_length = a.arclen(0.000_1);
@@ -514,9 +503,9 @@ mod tests {
                     let bez_length = a.path_segments(0.000_1).perimeter(0.000_1);
 
                     assert!(
-                         (arc_length - bez_length).abs() < EPSILON,
-                         "Numerically approximated arc length ({arc_length}) does not match bezier segment perimeter length ({bez_length}) for arc {a:?}"
-                     );
+                        (arc_length - bez_length).abs() < EPSILON,
+                        "Numerically approximated arc length ({arc_length}) does not match bezier segment perimeter length ({bez_length}) for arc {a:?}"
+                        );
                 }
             }
         }

From ef1dddd7051d3504d91fff15ae1beb09a3e26457 Mon Sep 17 00:00:00 2001
From: Thomas Churchman <thomas@kepow.org>
Date: Tue, 29 Oct 2024 14:24:20 +0100
Subject: [PATCH 17/30] Provide rem_euclid in FloatFuncs

---
 src/common.rs | 23 +++++++++++++++++++++++
 1 file changed, 23 insertions(+)

diff --git a/src/common.rs b/src/common.rs
index 30036869..1a89802b 100644
--- a/src/common.rs
+++ b/src/common.rs
@@ -34,6 +34,9 @@ macro_rules! define_float_funcs {
             /// Special implementation for signum, because libm doesn't have it.
             fn signum(self) -> Self;
 
+            /// Special implementation for `rem_euclid`, because libm doesn't have it.
+            fn rem_euclid(self, rhs: Self) -> Self;
+
             $(fn $name(self $(,$arg: $arg_ty)*) -> $ret;)+
         }
 
@@ -51,6 +54,16 @@ macro_rules! define_float_funcs {
                 }
             }
 
+            #[inline]
+            fn rem_euclid(self, rhs: f32) -> f32 {
+                let r = self % rhs;
+                if r < 0.0 {
+                    r + rhs.abs()
+                } else {
+                    r
+                }
+            }
+
             $(fn $name(self $(,$arg: $arg_ty)*) -> $ret {
                 #[cfg(feature = "libm")]
                 return libm::$lfname(self $(,$arg as _)*);
@@ -73,6 +86,16 @@ macro_rules! define_float_funcs {
                 }
             }
 
+            #[inline]
+            fn rem_euclid(self, rhs: f64) -> f64 {
+                let r = self % rhs;
+                if r < 0.0 {
+                    r + rhs.abs()
+                } else {
+                    r
+                }
+            }
+
             $(fn $name(self $(,$arg: $arg_ty)*) -> $ret {
                 #[cfg(feature = "libm")]
                 return libm::$lname(self $(,$arg as _)*);

From 3f36dc171279814ae5630f9ee29daf9e5093edc5 Mon Sep 17 00:00:00 2001
From: Thomas Churchman <thomas@kepow.org>
Date: Tue, 29 Oct 2024 14:31:25 +0100
Subject: [PATCH 18/30] clippy: group digits

---
 src/arc.rs | 8 ++++----
 1 file changed, 4 insertions(+), 4 deletions(-)

diff --git a/src/arc.rs b/src/arc.rs
index c33d1766..32c52c67 100644
--- a/src/arc.rs
+++ b/src/arc.rs
@@ -519,11 +519,11 @@ mod tests {
 
         // Numerical checks from section 3 of "Numerical computation of real or complex elliptic
         // integrals" (Carlson, Bille C.): https://arxiv.org/abs/math/9409227v1 (real-valued calls)
-        assert!((carlson_rf(1e-20, 1., 2., 0.) - 1.3110_28777_1461).abs() <= EPSILON);
-        assert!((carlson_rf(1e-20, 2., 3., 4.) - 0.58408_28416_7715).abs() <= EPSILON);
+        assert!((carlson_rf(1e-20, 1., 2., 0.) - 1.311_028_777_146_1).abs() <= EPSILON);
+        assert!((carlson_rf(1e-20, 2., 3., 4.) - 0.584_082_841_677_15).abs() <= EPSILON);
 
-        assert!((carlson_rd(1e-20, 0., 2., 1.) - 1.7972_10352_1034).abs() <= EPSILON);
-        assert!((carlson_rd(1e-20, 2., 3., 4.) - 0.16510_52729_4261).abs() <= EPSILON);
+        assert!((carlson_rd(1e-20, 0., 2., 1.) - 1.797_210_352_103_4).abs() <= EPSILON);
+        assert!((carlson_rd(1e-20, 2., 3., 4.) - 0.165_105_272_942_61).abs() <= EPSILON);
     }
 
     #[test]

From 73ddfff20a3250546a26fc2ec6562d9249fde4e3 Mon Sep 17 00:00:00 2001
From: Thomas Churchman <thomas@kepow.org>
Date: Tue, 29 Oct 2024 18:19:30 +0100
Subject: [PATCH 19/30] Fix erroneous agm_elliptic_perimeter if radius y >
 radius x

This doesn't (or doesn't easily?) occur for ellipses, as they instead
get a rotation when constructed with radius y > radius x.
---
 src/ellipse.rs | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/src/ellipse.rs b/src/ellipse.rs
index f5f2cc1a..ebc743a5 100644
--- a/src/ellipse.rs
+++ b/src/ellipse.rs
@@ -376,7 +376,7 @@ fn agm_elliptic_perimeter(accuracy: f64, radii: Vec2) -> f64 {
         Vec2::new(radii.y, radii.x)
     };
 
-    let accuracy = accuracy / (2. * PI * radii.x);
+    let accuracy = accuracy / (2. * PI * x);
 
     let mut sum = 1.;
     let mut a = 1.;
@@ -427,7 +427,7 @@ fn agm_elliptic_perimeter(accuracy: f64, radii: Vec2) -> f64 {
         a = a_next;
     }
 
-    2. * PI * radii.x / a * sum
+    2. * PI * x / a * sum
 }
 
 #[cfg(test)]

From 0234382b3d121ec69ccd91e237aa1a14232b00de Mon Sep 17 00:00:00 2001
From: Thomas Churchman <thomas@kepow.org>
Date: Tue, 29 Oct 2024 18:24:12 +0100
Subject: [PATCH 20/30] Special-case complete elliptic perimeter

---
 src/arc.rs     | 11 +++++++----
 src/ellipse.rs | 33 ++++++++++++++++++++-------------
 2 files changed, 27 insertions(+), 17 deletions(-)

diff --git a/src/arc.rs b/src/arc.rs
index 32c52c67..9769601a 100644
--- a/src/arc.rs
+++ b/src/arc.rs
@@ -3,7 +3,10 @@
 
 //! An ellipse arc.
 
-use crate::{Affine, Ellipse, ParamCurve, ParamCurveArclen, PathEl, Point, Rect, Shape, Vec2};
+use crate::{
+    ellipse::complete_elliptic_perimeter, Affine, Ellipse, ParamCurve, ParamCurveArclen, PathEl,
+    Point, Rect, Shape, Vec2,
+};
 use core::{
     f64::consts::{FRAC_PI_2, PI},
     iter,
@@ -259,12 +262,12 @@ impl ParamCurveArclen for Arc {
         } else {
             arclen += incomplete_elliptic_integral_second_kind(relative_error, end_angle, m);
         }
+        arclen *= radii.y;
 
         // Note: this uses the complete elliptic integral, which can be special-cased.
-        arclen +=
-            quarter_turns * incomplete_elliptic_integral_second_kind(relative_error, PI / 2., m);
+        arclen += 0.25 * quarter_turns * complete_elliptic_perimeter(self.radii, relative_error);
 
-        radii.y * arclen
+        arclen
     }
 }
 
diff --git a/src/ellipse.rs b/src/ellipse.rs
index ebc743a5..37c0bc96 100644
--- a/src/ellipse.rs
+++ b/src/ellipse.rs
@@ -247,19 +247,7 @@ impl Shape for Ellipse {
             return f64::NAN;
         }
 
-        // Check for the trivial case where the ellipse has one of its radii equal to 0, i.e.,
-        // where it describes a line, as the numerical method used breaks down with this extreme.
-        if radii.x == 0. || radii.y == 0. {
-            return 4. * f64::max(radii.x, radii.y);
-        }
-
-        // Evaluate an approximation based on a truncated infinite series. If it returns a good
-        // enough value, we do not need to iterate.
-        if kummer_elliptic_perimeter_range(radii) <= accuracy {
-            return kummer_elliptic_perimeter(radii);
-        }
-
-        agm_elliptic_perimeter(accuracy, radii)
+        complete_elliptic_perimeter(radii, accuracy)
     }
 
     fn winding(&self, pt: Point) -> i32 {
@@ -304,6 +292,25 @@ impl Shape for Ellipse {
     }
 }
 
+/// See [`Ellipse::perimeter`]. This is extracted to its own function for use by
+/// [`crate::Arc::perimeter`] without having to construct an [`Ellipse`].
+#[inline]
+pub(crate) fn complete_elliptic_perimeter(radii: Vec2, accuracy: f64) -> f64 {
+    // Check for the trivial case where the ellipse has one of its radii equal to 0, i.e.,
+    // where it describes a line, as the numerical method used breaks down with this extreme.
+    if radii.x == 0. || radii.y == 0. {
+        return 4. * f64::max(radii.x, radii.y);
+    }
+
+    // Evaluate an approximation based on a truncated infinite series. If it returns a good
+    // enough value, we do not need to iterate.
+    if kummer_elliptic_perimeter_range(radii) <= accuracy {
+        return kummer_elliptic_perimeter(radii);
+    }
+
+    agm_elliptic_perimeter(accuracy, radii)
+}
+
 /// Calculates circumference C of an ellipse with radii (x, y) as the infinite series
 ///
 /// C = pi (x+y) * ∑ binom(1/2, n)^2 * h^n from n = 0 to inf

From 6e111a86b3706128ec7ed176c94d866f922a2dcc Mon Sep 17 00:00:00 2001
From: Thomas Churchman <thomas@kepow.org>
Date: Tue, 29 Oct 2024 18:24:35 +0100
Subject: [PATCH 21/30] Lift calculation to compile-time

---
 src/arc.rs | 3 +--
 1 file changed, 1 insertion(+), 2 deletions(-)

diff --git a/src/arc.rs b/src/arc.rs
index 9769601a..f439fd3e 100644
--- a/src/arc.rs
+++ b/src/arc.rs
@@ -231,10 +231,9 @@ impl ParamCurveArclen for Arc {
 
         // Normalize start angle to be on the upper half of the ellipse
         let start_angle = start_angle.rem_euclid(PI);
-
         let end_angle = start_angle + sweep_angle;
 
-        let mut quarter_turns = (end_angle / PI * 2.).trunc() - (start_angle / PI * 2.).trunc();
+        let mut quarter_turns = (2. / PI * end_angle).trunc() - (2. / PI * start_angle).trunc();
         let end_angle = end_angle % PI;
 
         // The elliptic arc length is equal to radii.y * (E(end_angle | m) - E(start_angle | m))

From 902b543c5f24e9410198a3faf2e6497ed88e78e5 Mon Sep 17 00:00:00 2001
From: Thomas Churchman <thomas@kepow.org>
Date: Tue, 29 Oct 2024 18:24:56 +0100
Subject: [PATCH 22/30] Don't calculate radii again

---
 src/ellipse.rs | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/ellipse.rs b/src/ellipse.rs
index 37c0bc96..f3ec4691 100644
--- a/src/ellipse.rs
+++ b/src/ellipse.rs
@@ -243,7 +243,7 @@ impl Shape for Ellipse {
         // Note: the radii are calculated from an inner affine map (`self.inner`), and may be NaN.
         // Currently, constructing an ellipse with infinite radii will produce an ellipse whose
         // calculated radii are NaN.
-        if !self.radii().is_finite() {
+        if !radii.is_finite() {
             return f64::NAN;
         }
 

From 56c049a9f308ddf6c087d22b2dbb03be859673b5 Mon Sep 17 00:00:00 2001
From: Thomas Churchman <thomas@kepow.org>
Date: Wed, 30 Oct 2024 14:57:26 +0100
Subject: [PATCH 23/30] Improve wording

---
 src/arc.rs | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/arc.rs b/src/arc.rs
index f439fd3e..60380873 100644
--- a/src/arc.rs
+++ b/src/arc.rs
@@ -412,7 +412,7 @@ fn carlson_rd(relative_error: f64, x: f64, y: f64, z: f64) -> f64 {
         + 3. * sum
 }
 
-/// Numerically approximate the incomplete elliptic integral of the second kind to `phi`
+/// Numerically approximate the incomplete elliptic integral of the second kind from 0 to `phi`
 /// parameterized by `m = k^2` in Legendre's trigonometric form.
 ///
 /// Assumes:

From f57a05344353fd0413293a3535810e1d3975da18 Mon Sep 17 00:00:00 2001
From: Thomas Churchman <thomas@kepow.org>
Date: Fri, 1 Nov 2024 15:39:54 +0100
Subject: [PATCH 24/30] Add debug asserts

---
 src/arc.rs | 7 +++++++
 1 file changed, 7 insertions(+)

diff --git a/src/arc.rs b/src/arc.rs
index 60380873..3a6e2e4b 100644
--- a/src/arc.rs
+++ b/src/arc.rs
@@ -327,6 +327,9 @@ impl Mul<Arc> for Affine {
 ///
 /// RF = 1/2 ∫ 1 / ( sqrt(t+x) sqrt(t+y) sqrt(t+z) ) dt from 0 to inf
 fn carlson_rf(relative_error: f64, x: f64, y: f64, z: f64) -> f64 {
+    // At most one of (x, y, z) may be 0.
+    debug_assert!((x == 0.) as u8 + (y == 0.) as u8 + (z == 0.) as u8 <= 1);
+
     let mut x = x;
     let mut y = y;
     let mut z = z;
@@ -368,6 +371,10 @@ fn carlson_rf(relative_error: f64, x: f64, y: f64, z: f64) -> f64 {
 ///
 /// RD = 3/2 ∫ 1 / ( sqrt(t+x) sqrt(t+y) (t+z)^(3/2) ) dt from 0 to inf
 fn carlson_rd(relative_error: f64, x: f64, y: f64, z: f64) -> f64 {
+    // At most one of (x, y) may be 0, z must be nonzero.
+    debug_assert!(z != 0.);
+    debug_assert!(x != 0. || y != 0.);
+
     let mut x = x;
     let mut y = y;
     let mut z = z;

From 68cc0c107e7d512094bc84df64c582acbf37e004 Mon Sep 17 00:00:00 2001
From: Thomas Churchman <thomas@kepow.org>
Date: Fri, 1 Nov 2024 15:44:35 +0100
Subject: [PATCH 25/30] Merge sqrts

---
 src/arc.rs | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/src/arc.rs b/src/arc.rs
index 3a6e2e4b..b550724b 100644
--- a/src/arc.rs
+++ b/src/arc.rs
@@ -345,7 +345,7 @@ fn carlson_rf(relative_error: f64, x: f64, y: f64, z: f64) -> f64 {
             break;
         }
 
-        let lambda = x.sqrt() * y.sqrt() + x.sqrt() * z.sqrt() + y.sqrt() * z.sqrt();
+        let lambda = (x * y).sqrt() + (x * z).sqrt() + (y * z).sqrt();
         a = (a + lambda) / 4.;
         x = (x + lambda) / 4.;
         y = (y + lambda) / 4.;
@@ -391,7 +391,7 @@ fn carlson_rd(relative_error: f64, x: f64, y: f64, z: f64) -> f64 {
             break;
         }
 
-        let lambda = x.sqrt() * y.sqrt() + x.sqrt() * z.sqrt() + y.sqrt() * z.sqrt();
+        let lambda = (x * y).sqrt() + (x * z).sqrt() + (y * z).sqrt();
         sum += 4f64.powi(-m) / (z.sqrt() * (z + lambda));
         a = (a + lambda) / 4.;
         x = (x + lambda) / 4.;

From 9779d4e968e5646a12310dc44d926c1becf81e9c Mon Sep 17 00:00:00 2001
From: Thomas Churchman <thomas@kepow.org>
Date: Fri, 13 Dec 2024 21:26:57 +0100
Subject: [PATCH 26/30] Calculate incomplete integrals to absolute accuracy

---
 src/arc.rs | 230 ++++++++++++++++++++++++++++++++++++-----------------
 1 file changed, 157 insertions(+), 73 deletions(-)

diff --git a/src/arc.rs b/src/arc.rs
index b550724b..9814f439 100644
--- a/src/arc.rs
+++ b/src/arc.rs
@@ -8,7 +8,10 @@ use crate::{
     Point, Rect, Shape, Vec2,
 };
 use core::{
-    f64::consts::{FRAC_PI_2, PI},
+    f64::{
+        self,
+        consts::{FRAC_PI_2, PI},
+    },
     iter,
     ops::{Mul, Range},
 };
@@ -205,11 +208,7 @@ impl ParamCurve for Arc {
 }
 
 impl ParamCurveArclen for Arc {
-    fn arclen(&self, _accuracy: f64) -> f64 {
-        // TODO: wire up accuracy. The Carlson numerical approximation provides a bound on the relative
-        // error
-        let relative_error = 1e-20;
-
+    fn arclen(&self, accuracy: f64) -> f64 {
         // Normalize ellipse to have radius y >= radius x, required for the parameter assumptions
         // of `incomplete_elliptic_integral_second_kind`.
         let (radii, mut start_angle) = if self.radii.y >= self.radii.x {
@@ -248,23 +247,43 @@ impl ParamCurveArclen for Arc {
         // that range.
         let mut arclen = 0.;
 
+        // The available accuracy (tolerance) is distributed over the calculation of the two
+        // incomplete and one complete elliptic integrals.
+        let accuracy_per_incomplete_integral = 1. / 3. * accuracy / radii.y;
         if start_angle >= PI / 2. {
-            arclen += incomplete_elliptic_integral_second_kind(relative_error, PI - start_angle, m);
+            arclen += incomplete_elliptic_integral_second_kind(
+                accuracy_per_incomplete_integral,
+                PI - start_angle,
+                m,
+            );
             quarter_turns -= 1.;
         } else {
-            arclen -= incomplete_elliptic_integral_second_kind(relative_error, start_angle, m);
+            arclen -= incomplete_elliptic_integral_second_kind(
+                accuracy_per_incomplete_integral,
+                start_angle,
+                m,
+            );
         }
 
         if end_angle >= PI / 2. {
-            arclen -= incomplete_elliptic_integral_second_kind(relative_error, PI - end_angle, m);
+            arclen -= incomplete_elliptic_integral_second_kind(
+                accuracy_per_incomplete_integral,
+                PI - end_angle,
+                m,
+            );
             quarter_turns += 1.;
         } else {
-            arclen += incomplete_elliptic_integral_second_kind(relative_error, end_angle, m);
+            arclen += incomplete_elliptic_integral_second_kind(
+                accuracy_per_incomplete_integral,
+                end_angle,
+                m,
+            );
         }
         arclen *= radii.y;
 
-        // Note: this uses the complete elliptic integral, which can be special-cased.
-        arclen += 0.25 * quarter_turns * complete_elliptic_perimeter(self.radii, relative_error);
+        arclen += 1. / 4.
+            * quarter_turns
+            * complete_elliptic_perimeter(radii, 1. / 4. / 3. * accuracy * quarter_turns.max(1.));
 
         arclen
     }
@@ -326,22 +345,49 @@ impl Mul<Arc> for Affine {
 /// elliptic integrals" (Carlson, Bille C.): <https://arxiv.org/abs/math/9409227v1>
 ///
 /// RF = 1/2 ∫ 1 / ( sqrt(t+x) sqrt(t+y) sqrt(t+z) ) dt from 0 to inf
-fn carlson_rf(relative_error: f64, x: f64, y: f64, z: f64) -> f64 {
+fn carlson_rf(accuracy: f64, x: f64, y: f64, z: f64) -> f64 {
     // At most one of (x, y, z) may be 0.
     debug_assert!((x == 0.) as u8 + (y == 0.) as u8 + (z == 0.) as u8 <= 1);
 
+    // This mostly follows "Numerical computation of real or complex elliptic integrals", but using
+    // an absolute upper error bound rather than a relative one.
+    //
+    // From "Numerical computation of real or complex elliptic integrals" we have
+    //
+    // X_n = (a_0 - x_0) / (4^n a_n)
+    // (and the same for variables (Y,y), (Z,z)).
+    //
+    // From "Computing Elliptic Integrals by Duplication" we have an upper error bound of
+    //
+    // |err_n| < a_n^(-1/2) epsilon_n^6 / (4 (1 - epsilon_n))
+    // with epsilon_n = max(X_n, Y_n, Z_n)
+    //                = max(a_0 - x_0, a_0 - y_0, a_0 - z_0) / (4^n a_n).
+    //
+    // Define e_0 = max(a_0 - x_0, a_0 - y_0, a_0 - z_0). Rewrite for ease of computation,
+    //
+    //    |err_n| < a_n^(-1/2) epsilon_n^6 / (4 (1 - epsilon_n))
+    //            = a_n^(-1/2) e_0^6 / (4^n a_n)^6 / (4 (1 - epsilon_n))
+    // -> |err_n| a_n^(1/2) (4^n a_n)^6  / e_0^6 < 1 / (4 (1 - epsilon_n))
+    // -> |err_n| a_n^(1/2) a_n^6 4^(6n + 1) / e_0^6 < 1 / (1 - epsilon_n)
+    // -> |err_n| a_n^(1/2) a_n^6 4^(6n + 1) / e_0^6 (1 - epsilon_n) < 1.
+    //
+    // To reach an error upper bound of `accuracy`, iterate until
+    // 1 <= accuracy * a_n^(1/2) a_n^6 4^(6n + 1) / e_0^6 (1 - epsilon_n).
+
     let mut x = x;
     let mut y = y;
     let mut z = z;
 
-    let a0 = (x + y + z) / 3.;
-    let mut q = (3. * relative_error).powf(-1. / 6.)
-        * (a0 - x).abs().max((a0 - y).abs()).max((a0 - z).abs());
+    let mut a = (x + y + z) / 3.;
 
-    let mut a = a0;
-    let mut m = 0;
+    // These are partial terms of the inequality derived above. The multiply by (powers of) 4 are
+    // performed per iteration for computational efficiency.
+    let mut e = a - x.min(y).min(z);
+    let mut r = accuracy * 4. * e.powi(-6);
+
+    // let mut q = 1. / (3. * f64::EPSILON).cbrt().sqrt() * (a - x.min(y).min(z));
     loop {
-        if q <= a.abs() {
+        if 1. <= r * a.powi(6) * a.sqrt() * (1. - e / a) {
             break;
         }
 
@@ -351,81 +397,114 @@ fn carlson_rf(relative_error: f64, x: f64, y: f64, z: f64) -> f64 {
         y = (y + lambda) / 4.;
         z = (z + lambda) / 4.;
 
-        q /= 4.;
-        m += 1;
+        r *= 4f64.powi(6);
+        e /= 4.;
     }
 
-    let x = (a0 - x) / 4f64.powi(m) * a;
-    let y = (a0 - y) / 4f64.powi(m) * a;
+    let x = 1. - x / a;
+    let y = 1. - y / a;
     let z = -x - y;
 
     let e2 = x * y - z.powi(2);
     let e3 = x * y * z;
 
-    1. / a.sqrt()
-        * (1. - 1. / 10. * e2 + 1. / 14. * e3 + 1. / 24. * e2.powi(2) - 3. / 44. * e2 * e3)
+    (1. + (-1. / 10. * e2 + 1. / 14. * e3 + 1. / 24. * e2.powi(2) - 3. / 44. * e2 * e3)) / a.sqrt()
 }
 
 /// Approximation of the Carlson RD function as defined in "Numerical computation of real or
 /// complex elliptic integrals" (Carlson, Bille C.): <https://arxiv.org/abs/math/9409227v1>
 ///
 /// RD = 3/2 ∫ 1 / ( sqrt(t+x) sqrt(t+y) (t+z)^(3/2) ) dt from 0 to inf
-fn carlson_rd(relative_error: f64, x: f64, y: f64, z: f64) -> f64 {
+fn carlson_rd(accuracy: f64, x: f64, y: f64, z: f64) -> f64 {
     // At most one of (x, y) may be 0, z must be nonzero.
     debug_assert!(z != 0.);
     debug_assert!(x != 0. || y != 0.);
 
+    // As above for RF, find the absolute upper error bound rather than a relative one, the
+    // derivation of which is along the same lines.
+    //
+    // Again,
+    //
+    // X_n = (a_0 - x_0) / (4^n a_n)
+    // (and the same for variables (Y,y), (Z,z)).
+    //
+    // From "Computing Elliptic Integrals by Duplication" we have
+    //
+    // |err_n| < 4^-n a_n^(-3/2) 3 epsilon_n^6 / (1 - epsilon_n)^(3/2)
+    // with epsilon_n = max(X_n, Y_n, Z_n)
+    //                = max(a_0 - x_0, a_0 - y_0, a_0 - z_0) / (4^n a_n).
+    //
+    // Define e_0 = max(a_0 - x_0, a_0 - y_0, a_0 - z_0). Rewriting for ease of computation,
+    //
+    // |err_n| < 4^-n a_n^(-3/2) 3 epsilon_n^6 / (1 - epsilon_n)^(3/2)
+    //         = 4^-n a_n^(-3/2) 3 e_0^6 / 4^(6n) a_n^6 / (1 - epsilon_n)^(3/2)
+    // -> |err_n| 4^(7n) a_n^(3/2) a_n^6 / (3 e_0^6) < 1 / (1 - epsilon_n)^(3/2)
+    // -> |err_n| 4^(7n) a_n^(3/2) a_n^6 (1/3) / e_0^6 < (1 / 1 - epsilon_n)^(3/2),
+    // raise to the power 2/3,
+    // -> |err_n|^(2/3) 4^(14/3 n) a_n a_n^4 (1/3)^(2/3) / e_0^4 < 1 / (1 - epsilon_n)
+    // -> |err_n|^(2/3) 4^(14/3 n) a_n^5 (1/3)^(2/3) / e_0^4 (1 - epsilon_n) < 1
+    //
+    // That means, to reach an error upper bound of `accuracy`, iterate until
+    // 1 <= accuracy^(2/3) 4^(14/3 n) a_n^5 (1/3)^(2/3) / e_0^4 (1 - epsilon)
+
     let mut x = x;
     let mut y = y;
     let mut z = z;
 
     let a0 = (x + y + 3. * z) / 5.;
-    let mut q = (relative_error / 4.).powf(-1. / 6.)
-        * (a0 - x).abs().max((a0 - y).abs()).max((a0 - z).abs());
+    let mut a = a0;
 
     let mut sum = 0.;
-    let mut a = a0;
-    let mut m = 0;
+    let mut mul = 1.;
+
+    // These are partial terms of the inequality derived above. The multiply by (powers of) 4 are
+    // performed per iteration for computational efficiency.
+    let mut e = a - x.min(y).min(z);
+    let mut r = (accuracy / 3.).powf(2. / 3.) * e.powi(-4);
+
     loop {
-        if q <= a.abs() {
+        if 1. <= r * a.powi(5) * (1. - e / a) {
             break;
         }
 
         let lambda = (x * y).sqrt() + (x * z).sqrt() + (y * z).sqrt();
-        sum += 4f64.powi(-m) / (z.sqrt() * (z + lambda));
+        sum += mul / (z.sqrt() * (z + lambda));
         a = (a + lambda) / 4.;
         x = (x + lambda) / 4.;
         y = (y + lambda) / 4.;
         z = (z + lambda) / 4.;
 
-        q /= 4.;
-        m += 1;
+        r *= 4f64.powf(14. / 3.);
+        e /= 4.;
+        mul /= 4.;
     }
 
-    let x = (a0 - x) / (4f64.powi(4) * a);
-    let y = (a0 - y) / (4f64.powi(4) * a);
-    let z = -(x + y) / 3.;
+    let x = 1. - x / a;
+    let y = 1. - y / a;
+    let z = (-x - y) / 3.;
 
     let e2 = x * y - 6. * z.powi(2);
     let e3 = (3. * x * y - 8. * z.powi(2)) * z;
-    let e4 = 3. * x * y - z.powi(2) * z.powi(2);
+    let e4 = 3. * (x * y - z.powi(2)) * z.powi(2);
     let e5 = x * y * z.powi(3);
 
-    4f64.powi(-m) * 1. / (a * a.sqrt())
-        * (1. - 3. / 14. * e2 + 1. / 6. * e3 + 9. / 88. * e2.powi(2)
-            - 3. / 22. * e4
-            - 9. / 52. * e2 * e3
-            + 3. / 26. * e5)
+    (1. - 3. / 14. * e2 + 1. / 6. * e3 + 9. / 88. * e2.powi(2) - 3. / 22. * e4 - 9. / 52. * e2 * e3
+        + 3. / 26. * e5)
+        * mul
+        / (a * a.sqrt())
         + 3. * sum
 }
 
 /// Numerically approximate the incomplete elliptic integral of the second kind from 0 to `phi`
 /// parameterized by `m = k^2` in Legendre's trigonometric form.
 ///
+/// The absolute error between the calculated integral and the true integral is bounded by
+/// `accuracy` (modulo floating point rounding errors).
+///
 /// Assumes:
 /// 0 <= phi <= pi / 2
 /// and 0 <= m sin^2(phi) <= 1
-fn incomplete_elliptic_integral_second_kind(relative_error: f64, phi: f64, m: f64) -> f64 {
+fn incomplete_elliptic_integral_second_kind(accuracy: f64, phi: f64, m: f64) -> f64 {
     // Approximate the incomplete elliptic integral through Carlson symmetric forms:
     // https://en.wikipedia.org/w/index.php?title=Carlson_symmetric_form&oldid=1223277638#Incomplete_elliptic_integrals
 
@@ -442,8 +521,25 @@ fn incomplete_elliptic_integral_second_kind(relative_error: f64, phi: f64, m: f6
     // note: this actually allows calculating from -1/2 pi <= phi <= 1/2 pi, but there are some
     // alternative translations from the Legendre form that are potentially better, that do
     // restrict the domain to 0 <= phi <= 1/2 pi.
-    sin * carlson_rf(relative_error, cos2, 1. - m * sin2, 1.)
-        - 1. / 3. * m * sin3 * carlson_rd(relative_error, cos2, 1. - m * sin2, 1.)
+    let term1 = if sin == 0. {
+        0.
+    } else {
+        sin * carlson_rf(
+            accuracy / (2. * sin),
+            // 1e-30,
+            cos2,
+            1. - m * sin2,
+            1.,
+        )
+    };
+
+    let term2 = if sin == 0. || m == 0. {
+        0.
+    } else {
+        1. / 3. * m * sin3 * carlson_rd(accuracy * 3. / 2. / (m * sin3), cos2, 1. - m * sin2, 1.)
+    };
+
+    term1 - term2
 }
 
 #[cfg(test)]
@@ -468,10 +564,6 @@ mod tests {
 
     #[test]
     fn length() {
-        // TODO: when arclen actually uses specified accuracy, update EPSILON and the accuracy
-        // params
-        const EPSILON: f64 = 1e-6;
-
         // Circular checks:
         for (start_angle, sweep_angle, length) in [
             (0., 1., 1.),
@@ -482,37 +574,35 @@ mod tests {
             (2.5, 10., 10.),
         ] {
             let a = Arc::new((0., 0.), (1., 1.), start_angle, sweep_angle, 0.);
-            let arc_length = a.arclen(0.000_1);
+            let arc_length = a.arclen(1e-7);
             assert!(
-                (arc_length - length).abs() <= EPSILON,
+                (arc_length - length).abs() <= 1e-6,
                 "Got arc length {arc_length}, expected {length} for circular arc {a:?}"
             );
         }
 
         let a = Arc::new((0., 0.), (1., 1.), 0., PI * 4., 0.);
-        assert!((a.arclen(0.000_1) - PI * 4.).abs() <= EPSILON);
+        assert!((a.arclen(1e-13) - PI * 4.).abs() <= 1e-12);
 
         let a = Arc::new((0., 0.), (2.23, 3.05), 0., 0.2, 0.);
-        assert!((a.arclen(0.000_1) - 0.60811714277).abs() <= EPSILON);
+        assert!((a.arclen(1e-13) - 0.608_117_142_773_153_8).abs() <= 1e-12);
 
         let a = Arc::new((0., 0.), (3.05, 2.23), 0., 0.2, 0.);
-        assert!((a.arclen(0.000_1) - 0.448555).abs() <= EPSILON);
+        assert!((a.arclen(1e-13) - 0.448_554_961_296_305_9).abs() <= 1e-12);
     }
 
     #[test]
     fn length_compare_with_bez_length() {
-        const EPSILON: f64 = 1e-3;
-
         for radii in [(1., 1.), (0.5, 1.), (2., 1.)] {
             for start_angle in [0., 0.5, 1., 2., PI, -1.] {
                 for sweep_angle in [0., 0.5, 1., 2., PI, -1.] {
                     let a = Arc::new((0., 0.), radii, start_angle, sweep_angle, 0.);
 
-                    let arc_length = a.arclen(0.000_1);
-                    let bez_length = a.path_segments(0.000_1).perimeter(0.000_1);
+                    let arc_length = a.arclen(1e-8);
+                    let bez_length = a.path_segments(1e-8).perimeter(1e-8);
 
                     assert!(
-                        (arc_length - bez_length).abs() < EPSILON,
+                        (arc_length - bez_length).abs() < 1e-7,
                         "Numerically approximated arc length ({arc_length}) does not match bezier segment perimeter length ({bez_length}) for arc {a:?}"
                         );
                 }
@@ -522,33 +612,27 @@ mod tests {
 
     #[test]
     fn carlson_numerical_checks() {
-        // TODO: relative bound on error doesn't seem to be quite correct yet, use a large epsilon
-        // for now
-        const EPSILON: f64 = 1e-6;
-
         // Numerical checks from section 3 of "Numerical computation of real or complex elliptic
         // integrals" (Carlson, Bille C.): https://arxiv.org/abs/math/9409227v1 (real-valued calls)
-        assert!((carlson_rf(1e-20, 1., 2., 0.) - 1.311_028_777_146_1).abs() <= EPSILON);
-        assert!((carlson_rf(1e-20, 2., 3., 4.) - 0.584_082_841_677_15).abs() <= EPSILON);
+        assert!((carlson_rf(1e-13, 1., 2., 0.) - 1.311_028_777_146_1).abs() <= 1e-12);
+        assert!((carlson_rf(1e-13, 2., 3., 4.) - 0.584_082_841_677_15).abs() <= 1e-12);
 
-        assert!((carlson_rd(1e-20, 0., 2., 1.) - 1.797_210_352_103_4).abs() <= EPSILON);
-        assert!((carlson_rd(1e-20, 2., 3., 4.) - 0.165_105_272_942_61).abs() <= EPSILON);
+        assert!((carlson_rd(1e-13, 0., 2., 1.) - 1.797_210_352_103_4).abs() <= 1e-12);
+        assert!((carlson_rd(1e-13, 2., 3., 4.) - 0.165_105_272_942_61).abs() <= 1e-12);
     }
 
     #[test]
     fn elliptic_e_numerical_checks() {
-        const EPSILON: f64 = 1e-6;
-
         for (phi, m, elliptic_e) in [
             (0.0, 0.0, 0.0),
             (0.5, 0.0, 0.5),
             (1.0, 0.0, 1.0),
             (0.0, 1.0, 0.0),
-            (1.0, 1.0, 0.84147098),
+            (1.0, 1.0, 0.841_470_984_807_896_5),
         ] {
-            let elliptic_e_approx = incomplete_elliptic_integral_second_kind(1e-20, phi, m);
+            let elliptic_e_approx = incomplete_elliptic_integral_second_kind(1e-13, phi, m);
             assert!(
-                (elliptic_e_approx - elliptic_e).abs() < EPSILON,
+                (elliptic_e_approx - elliptic_e).abs() < 1e-12,
                 "Approximated elliptic e {elliptic_e_approx} does not match known value {elliptic_e} for E({phi}|{m})"
             );
         }

From 6fab9bcbf8c300f8f8496b8e5c157c0d94156cef Mon Sep 17 00:00:00 2001
From: Thomas Churchman <thomas@kepow.org>
Date: Fri, 13 Dec 2024 21:34:47 +0100
Subject: [PATCH 27/30] Add changelog

---
 CHANGELOG.md | 2 ++
 1 file changed, 2 insertions(+)

diff --git a/CHANGELOG.md b/CHANGELOG.md
index 128cc582..fe7b07f2 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -26,6 +26,7 @@ This release has an [MSRV][] of 1.65.
 
 - Reduce number of operations in `Triangle::circumscribed_circle`. ([#390][] by [@tomcur][])
 - Numerically approximate ellipse perimeter. ([#383] by [@tomcur][])
+- Numerically approximate elliptic arc perimeter. ([#381] by [@tomcur][])
 
 ## [0.11.1][] (2024-09-12)
 
@@ -93,6 +94,7 @@ Note: A changelog was not kept for or before this release
 [#376]: https://github.com/linebender/kurbo/pull/376
 [#378]: https://github.com/linebender/kurbo/pull/378
 [#379]: https://github.com/linebender/kurbo/pull/379
+[#381]: https://github.com/linebender/kurbo/pull/381
 [#383]: https://github.com/linebender/kurbo/pull/383
 [#388]: https://github.com/linebender/kurbo/pull/388
 [#390]: https://github.com/linebender/kurbo/pull/390

From 4af43fce604e50db444ce2801f878be5049bd300 Mon Sep 17 00:00:00 2001
From: Thomas Churchman <thomas@kepow.org>
Date: Fri, 13 Dec 2024 21:40:26 +0100
Subject: [PATCH 28/30] Remove comment

---
 src/arc.rs | 1 -
 1 file changed, 1 deletion(-)

diff --git a/src/arc.rs b/src/arc.rs
index 9814f439..75db53ed 100644
--- a/src/arc.rs
+++ b/src/arc.rs
@@ -385,7 +385,6 @@ fn carlson_rf(accuracy: f64, x: f64, y: f64, z: f64) -> f64 {
     let mut e = a - x.min(y).min(z);
     let mut r = accuracy * 4. * e.powi(-6);
 
-    // let mut q = 1. / (3. * f64::EPSILON).cbrt().sqrt() * (a - x.min(y).min(z));
     loop {
         if 1. <= r * a.powi(6) * a.sqrt() * (1. - e / a) {
             break;

From 106be0cc66378e8543957893a78a938d523291c0 Mon Sep 17 00:00:00 2001
From: Thomas Churchman <thomas@kepow.org>
Date: Fri, 13 Dec 2024 22:02:47 +0100
Subject: [PATCH 29/30] Group operations

---
 src/arc.rs | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/arc.rs b/src/arc.rs
index 75db53ed..3bf30d58 100644
--- a/src/arc.rs
+++ b/src/arc.rs
@@ -535,7 +535,7 @@ fn incomplete_elliptic_integral_second_kind(accuracy: f64, phi: f64, m: f64) ->
     let term2 = if sin == 0. || m == 0. {
         0.
     } else {
-        1. / 3. * m * sin3 * carlson_rd(accuracy * 3. / 2. / (m * sin3), cos2, 1. - m * sin2, 1.)
+        1. / 3. * m * sin3 * carlson_rd(accuracy * (3. / 2.) / (m * sin3), cos2, 1. - m * sin2, 1.)
     };
 
     term1 - term2

From ad09dead30924605a1eb19041998aa6dbd6d9e57 Mon Sep 17 00:00:00 2001
From: Thomas Churchman <thomas@kepow.org>
Date: Fri, 13 Dec 2024 22:04:05 +0100
Subject: [PATCH 30/30] Recip

---
 src/arc.rs | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/src/arc.rs b/src/arc.rs
index 3bf30d58..2e4859b5 100644
--- a/src/arc.rs
+++ b/src/arc.rs
@@ -378,7 +378,7 @@ fn carlson_rf(accuracy: f64, x: f64, y: f64, z: f64) -> f64 {
     let mut y = y;
     let mut z = z;
 
-    let mut a = (x + y + z) / 3.;
+    let mut a = (x + y + z) * (1. / 3.);
 
     // These are partial terms of the inequality derived above. The multiply by (powers of) 4 are
     // performed per iteration for computational efficiency.
@@ -450,7 +450,7 @@ fn carlson_rd(accuracy: f64, x: f64, y: f64, z: f64) -> f64 {
     let mut y = y;
     let mut z = z;
 
-    let a0 = (x + y + 3. * z) / 5.;
+    let a0 = (x + y + 3. * z) * (1. / 5.);
     let mut a = a0;
 
     let mut sum = 0.;