Introduction
A year ago, in a Review [1] of the approximate formulae for ellipse perimeters, I have launched the following challenge:
Those who wish to pursue this topic further should look for an efficient formula using only the four algebraic operations (if possible, avoiding even squareroot) with a maximum error smaller than 10 ppm. If would be also nice if such a formula were exact for both the circle (y=1) and the degenerate flat ellipse (y=0).
Shortly afterwards, David W.Cantrell published an approximation [2] which meets all the above criteria and has the maximum relative error of just 5.2 ppm. Consequently, in an Addendum [3] to the Review, I have strengthened the challenge calling for a maximum error of less than 1 ppm.
This Note describes a discrete family of approximations indexed by a degree n = 1, 2, 3,..., all of which meet the qualitative criteria of the challenge. Moreover, for n &ge, 6 they also meet the 1 ppm maximum error requirement. It is interesting that three of these approximations (those for n = 1, 3, and 5) correspond to already known formulae (Cantrell's Particularly Fruitful, Sykora 2005, and Cantrell 2006, respectively).
Since the Note adopts the same nomenclature as the Review [1], the reader is kindly requested to consult the latter in case of any doubt regarding symbol definitions and their elementary relations.
General properties of algebraic approximations
We will call algebraic any expression which involves a finite number of elements of an algebraic field and uses just the arithmetic operations of sum (+), subtraction (), product (*) and division (/). It is a matter of elementary considerations to see that any such expression can be written as a rational function (i.e., a ratio of two polynomials) in terms of the input elements.
In the case of approximate algebraic expressions for the perimeter of an ellipse S(a,b) in terms of the semiaxes a and b, we are therefore looking for an expression of the type
(1) ,
where N(a,b) and D(a,b) are polynomials in a,b.
This requirement needs to be complemented by proper dimensional (or scaling) behavior. Since a and b are quantities with the dimension of length, all the terms of each of the two polynomials must have the same total degree and, since S(a,b) has also a dimension of length, the degree of N(a,b) must be higher by 1 than that of D(a,b). In addition to all this, we require the approximations to be invariant with respect to the exchange of a and b, regardless of the numerical values of the adjustable coefficients. This is only possible if, for any admissible pair of exponents, terms such as a^{p}b^{q} and a^{q}b^{p} have the same coefficient.
Consequently, the algebraic approximations can be classified according to the degree of D(a,b)  henceforth called the order of the algebraic approximation  and they must be of the form (1) with
(2) ,
where n is order of the approximation and r_{k}, s_{k} are (n+2) coefficients whose values need to be defined.
Any approximation formula for S(a,b) of this form can be rewritten in terms of the ratio y = b/a and leads also to an approximation for the complete elliptic integral E(x):
(3) ,
where t_{k} = r_{k}/4 and y^{2} = 1x.
Since the value of expression (3) does not change when all the coefficients are divided by a common factor, it is possible to normalize it so as s_{0} = 1. In the following we will assume that this has been done. The number of coefficients to be determined is therefore equal to (n+1).
Types of algebraic approximation formulae
The set of conditions used to define the coefficients determines what kind of algebraic approximations one obtains. Thus, for example, matching the (n+1) derivatives of E(x) at y = 1 from 0th to nth optimizes the formulae for ellipses with low eccentricities and leads to a subselection of Padé approximations [1, section K). A similar approach is not possible for flat ellipses with large eccentricities since S(a,b)/4a does not have a convergent Taylor expansion at y = 0 (the behavior .
Another family of approximations arises applying the requirement that the formulae be exact for both y = 0 (flat ellipses, perimeter equal to 4a) and y = 1 (round ellipses, perimeter equal to 2πa).
The first condition implies t_{0} = 1, while the second one amounts to
(4) ,
a relation which can be used to determine, for example, the highest of the coefficients t_{k}. It is interesting that these same conditions guarantee also a match for the first derivative of S(a,b)/4a at y = 1, while this is not true at y = 0. This explains why perimeters of roundish ellipses with low eccentricities get in any case approximated better than perimeters of flat elipses.
With these two constraints, the number of undefined coefficients becomes (n1) meaning, for example, that there is nothing left to adjust in the case n = 1.
When n > 1, one can proceed in two ways. The first one consists in matching the 2nd to (n+1)th derivatives of S(a,b)/4a at y = 1 (0th and 1st derivatives are already matched). In such a case one would obtain Padélike approximations constrained to the exact value at y = 0. This is reminiscent of the linear Padé combinations introduced in [1] in order to achieve essentially the same goal (the extent to which the two approaches might be equivalent is at this moment not clear).
In the second, numerically more brutal approach, all the (n1) undefined coefficients are adjusted so as to minimize the maximum relative error of the approximation over the whole interval [0,1] of y values. This prescription defines a unique set of coefficient values and leads to what we will call globally best algebraic approximations or GBA's.
One might also consider hybrid algebraic approximations in which one matches m derivatives of S(a,b)/4a at y = 1 and adjusts the remaining (n+1m) coefficients to minimize the maximum relative error. This leads to a large set of approximations classified according to their order n and their derivative index m = 0, 1, 2, ..., (n+1).
In this work, the coefficients of the GBA's (m = 0) were determined for orders n = 1 to 7. The respective error curves are shown in Figure 1 (for graph conventions, consult [1]).
Error curves for the globally best algebraic approximations of orders 1 to 7_{ }
The horizontal axis corresponds to y = b/a and is either linear (left graph) or logarithmic (right graph). The vertical axis corresponds to the relative error. The error curves of the GBA approximations of various order n are distinguished by color. Their global maximum is highest for n = 1 and smallest for n = 7. The line thickness of each trace is used to indicate the sign of the error: bold sections are positive, thin sections negative. The graphs were produced using the enclosed MATLAB program EllipseAlgApprox.m.
Particular cases
Order 1, rel.error < 6313 ppm
n =1 is the only case in which there are no freely adjustable coefficients. It coincides with the Cantrell's Particularly Fruitful approximation (see [3] for a discussion of its origin) and can be written as
(A1) .
Order 2, rel.error < 557.3 ppm
This approximation appears to be new:
(A2) .
Order 3, rel.error < 75.96 ppm
This approximation turns out to be equivalent to Sykora 2005 [1:C5] provided that its coefficients are numerically adjusted as suggested by Cantrell [3]. Using the formalism of this Note, it can be written as
(A3) .
Order 4, rel.error < 15.29 ppm
This approximation appears to be new:
(A4) .
Order 5, rel.error < 3.65 ppm
This approximation is equivalent to Cantrell 2006 [2,3] provided that its coefficients are better adjusted. In his original announcement [2], David Cantrell in fact mentions that he had problems with the numerical minimization and that his coefficients and his relative error bound of 5.2 ppm need not be the best. Using the formalism of this Note, the optimized formula can be written as
(A5) .
Expanding both this and Cantrell's 2006 formulas as plain rational functions and comparing the coefficients, one obtains the relations between our t_{k}'s and s_{k}'s and the Cantrell's coefficients s, p, q, r, t. It follows that the optimal values of the latter should have been
s = 22.29785433367,
t = 33.65830998486,
p = 3.98644261166,
q = 82.24320600580,
r = 84.00260009659.
Order 6, rel.error < 0.988 ppm
This approximation is new and it is the first one with guaranteed maximum error smaller than 1 ppm:
(A6) .
Order 7, rel.error < 0.282 ppm
This, of course, is also a new approximation:
(A7) .
Discussion
This Note intentionally avoids the problem of optimizing the approximations from the algorithmic point of view. As presented here, the expressions maintain the elegance originating from a single guiding principle. Minimizing the number of arithmetic operations in order to optimize computational efficiency is another matter which might well hide the fact that the expressions form a single, internally coherent family.
The numerical values of the coefficients reported here are those obtained from the optimization process. I have tried hard to make sure that the reported maximum errors be correct and I do not believe one can improve them except, perhaps, beyond the fourth significant decimal digit. The coefficient values, however, do have a certain tolerance range and could be rounded to a lower number of significant digits (maybe two more than the maximum relative error). It is yet to be seen, for example, whether there is a set of singleprecision coefficients which would maintain the precision of the 6th and 7th order algebraic approximations. I think that it is probable but I have not investigated the matter in detail.
For approximations of 5th and higher orders, the numeric optimization process is very slow and rather unreliable in the sense the local minima due to limited precision and/or discrete errorcurve sampling all but preclude a straightforward minimization. I had to resort to alternate minimization of the absolute value of maximum relative error and the mean square of the relative deviations. Such an alternation, reminiscent of simulated annealing, helps to nudge the process out of false local minima and drive it towards the absolute one. It is also helpful to use Eq.(4) to compute the highest t_{k} coefficient rather than t_{1}.
Determining the coefficients for still higherorder GBA approximations may be an arduous task (a nice exercise in numeric optimization, though) which I, at least for the moment, refrain from undertaking. Extrapolating the current results, there is no doubt that the 8th order approximation has a maximum relative error smaller than 0.1 ppm (100 ppb).
