Calculate incomplete integrals to absolute accuracy

Author: tomcur

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})"
             );
         }