Fix math

src/ellipse.rs

@@ -401,7 +401,7 @@ fn agm_elliptic_perimeter(accuracy: f64, radii: Vec2) -> f64 {
     let mut c = (1. - g.powi(2)).sqrt();
     let mut mul = 0.5;
 
-    loop {
+    for m in 0.. {
         let c2 = c.powi(2);
         let term = mul * c2;
         sum -= term;
@@ -412,19 +412,22 @@ fn agm_elliptic_perimeter(accuracy: f64, radii: Vec2) -> f64 {
         // <https://web.archive.org/web/20241002140956/https://www.maths.lancs.ac.uk/jameson/ellagm.pdf>)
         //
         // Therefore
-        // ∑ 2^(i-1) c_i^2 from i = 0 ≤ ∑ 2^(i-1) ((1/2)^i c_0)^2 from i = 0
-        //                            = ∑ 2^-(i+1) c_0^2          from i = 0
-        //                            = c_0^2
+        // ∑ 2^(i-1) c_i^2 from i = 1 < ∑ 2^(i-1) ((1/2)^i c_0)^2 from i = 1
+        //                            = ∑ 2^-(i+1) c_0^2          from i = 1
+        //                            = 1/2 c_0^2
         //
-        // or, for arbitrary starting point n > 0:
-        // ∑ 2^(i-1) c_i^2 from i = n ≤ ∑ 2^(i-1) ((1/2)^(i-n) c_n)^2 from i = n
-        //                            = ∑ 2^(2n - i - 1) c_n^2        from i = n
-        //                            = 2^n ∑ 2^(-(i+1)) c_n^2        from i = n
-        //                            = 2^n 2^(2-n) c_n^2
-        //                            = 4 c_n^2
+        // or, for arbitrary starting point i = n+1:
+        // ∑ 2^(i-1) c_i^2 from i = n+1 < ∑ 2^(i-1) ((1/2)^(i-n) c_n)^2 from i = n+1
+        //                              = ∑ 2^(2n - i - 1) c_n^2        from i = n+1
+        //                              = 2^(2n) ∑ 2^(-(i+1)) c_n^2     from i = n+1
+        //                              = 2^(2n) 2^(-(n+1)) c_n^2
+        //                              = 2^(n-1) c_n^2
         //
-        // Therefore, the remainder of the series sums to less than 4 c_n^2.
-        if 4. * c2 <= accuracy {
+        // Therefore, the remainder of the series sums to less than 2^(n-1) c_n^2.
+        //
+        // TODO: is there a way to directly specify a power-of-two float, i.e., directly setting
+        // the exponent part?
+        if c2 * 2f64.powi(m-1) <= accuracy {
             // `sum` currently overestimates the true value - subtract the upper bound of the
             // remaining series. We will then underestimate the true value, but by no more than
             // `accuracy`.