diff --git a/src/ellipse.rs b/src/ellipse.rs index e453ca3d..7579a73d 100644 --- a/src/ellipse.rs +++ b/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 { // ) // // 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`.