Numerically approximate ellipse perimeter

[tomcur]

This implements the (truncated) infinite series for the ellipse perimeter as described by Kummer (1837) and rediscovered by Linderholm and Segal (1995), and the iterative arithmetic-geometric mean method.

In the common case of ellipses that are only moderately eccentric and for "normal" accuracy (say, 0.001), the infinite series converges quickly. The series is known at compile-time, truncated to the 6th power. For a given ellipse and accuracy, a quick check is performed at runtime, to determine whether the truncated series' approximation is within the desired accuracy. If so, the truncated series is evaluated.

If the check determines the approximation is not good enough, the problem is handed to the iterative arithmetic-geometric mean method. This method converges quadratically for all cases.

[tomcur]

This solves the general case up to arbitrary precision.

For ellipses that are only moderately eccentric, the series converges very quickly, such that only a few terms are needed for normal accuracy. Even for practically infinite accuracy (i.e., for Kurbo's case: within machine precision of an f64), only a few terms are needed when ellipses are just slightly eccentric.

Wikipedia states terms up to and including the 4th power are enough for eccentricities of 0.5 or less to reach the limits of f64 precision. With a back-of-the-envelope calculation that seems about right: the next term would be 49 * 0.005^5 / 65536, log2 of that is -48, so 48 bits away from the first term of 1.

We could expand the first few terms explicitly, only starting iteration when necessary. Though we'd have to benchmark whether that really is faster.

[tomcur]

The method I first added (kummer_elliptic_perimeter, based on a series with binomial coefficients in n=1/2) only converges quadratically for ellipses that are not too eccentric. I've added an iterative method based on the arithmetic-geometric mean that actually does converge quadratically for all cases, requiring 7 iterations even in the case of an ellipse with radii (0.000_000_1, 1.) and desired accuracy 0.000_001. The binomial series required 355 terms, with the difference increasing rapidly with higher eccentricity or better accuracy (with an accuracy of 0.000_000_01, it's 8 iterations vs. 3534 terms).

The method based on binomial coefficients is still nice, in the sense that the coefficients can be known at compile-time. We can expand that series up to some number of terms, evaluating that in the case the approximation is good enough for the given eccentricity and accuracy. Wikipedia states terms up to and including the 4th power are enough for eccentricities of 0.5 or less to reach the limits of f64 precision. With a back-of-the-envelope calculation that seems about right: the next term would be 49 * 0.005^5 / 65536, log2 of that is -48, so 48 bits away from the first term of 1.

That allows a loopless approximation for the most common cases, using the quadratically converging iterative arithmetic-geometric mean method for very high eccentricity or very precise accuracy.

[tomcur]

That allows a loopless approximation for the most common cases, using the quadratically converging iterative arithmetic-geometric mean method for very high eccentricity or very precise accuracy.

Done, and PR text updated. Perhaps we should benchmark whether this is actually enough of a performance improvement for common ellipses, to justify the additional lines of code.

[raphlinus]

Very cool, from what I can see the math checks out.

I think the AGM algorithm is probably very good, but for one other possibility that seems to work well in the high eccentricity case, see "Cayley" in The Perimeter of an Ellipse. That formula requires logarithms, though, so it's not clear it will be better than AGM. Also, I suspect the accuracy of table 1, as perimeter is not monotonic.

[tomcur]

one other possibility that seems to work well in the high eccentricity case, see "Cayley" in The Perimeter of an Ellipse. That formula requires logarithms, though, so it's not clear it will be better than AGM.

Very nice, that does seem to work well. It only requires evaluating ln(4/k), so the cost is constant. At the very least it's a good candidate to handle the same way as kummer_elliptic_perimeter, as the coefficients are known at compile-time, requiring just a few additions and multiplications per term.

I wonder if the overlap between Kummer and Cayley is enough to guarantee machine precision for all eccentricities, without too many terms.

Also, I suspect the accuracy of table 1, as perimeter is not monotonic.

Well-spotted. Having checked a few, it seems from b=0.002 onwards the reported perimeter decimals are accidentally shifted left by one. The decimals are otherwise correct.

[tomcur]

I think the Cayley formula is worth exploring, but will leave that for another PR. I haven't looked into it yet. In the meantime this is an improvement and is useful also for #381.