THE DERIVATION OF ALGORITHMS TO COMPUTE ELLIPTIC INTEGRALS OF THE FIRST AND SECOND KIND BY LANDEN-TRANSFORMATION

This paper deals with elliptic integrals of first and second kind and their solution by Landen-transformation . It is explained how the Landen-transformation works and how the different algorithms can be used in geodetic problem solving. The algorithms are recursive and thus easy to implement and stable. Further on the testing is quite simple. The numerical results of the Landen-transformation are compared to this computed with Maple V to show the exactness of the algorithms. For two algorithms the source-code is depicted.


INTRODUCTION
In 1771 and 1775 the British mathematician John Landen (1719-1790) published an appropriate (which means applicable at this time) method to solve elliptic integrals (see LANDEN (1771), LANDEN (1775)).This approach is now called the Landen-transformation or Landen-substitution respectively.Remarkably the above mentioned transformation was published before Joseph Liouville  proved that there exists no anti derivative for the elliptic integrals.However, the definite integral can be solved.
In addition to the Landen-transformation there are still other approaches to solve the elliptic integrals.These are for example the series expansion or the numerical integration.Both topics comprise a number of methods.Also in the group of the transformation/substitution methods another kind of transformation is known, the so-called Gauss-transformation.However, these approaches are beyond the scope of this paper.
In fact there are different ways to depict the elliptic integrals but most commonly the notation of Legendre is used.According to this notation the definition of the incomplete elliptic integrals is: The terms k and φ are denoted as modulus and amplitude respectively.
In mathematical and cartographic geodesy there are different problems where these integrals are found.Only a few examples of these should be shown in this introduction.First the computation of the elliptic arc length S should be mentioned which leads to (see HECK, (1995)) (a = semi-major axis of the ellipsoid of revolution; e = first eccentricity, φ 1 = latitude) where the first term of the formula includes an elliptic integral of second kind with the modulus e. Hopfner explained in (HOPFNER,(1938)) how the arc length can be computed by different methods.In this context, he also mentioned the Landen-transformation.
Especially formula 3 reveals that the integral has to be computed with a high accuracy, because the result of E(e, φ 1 ) is multiplied with the semi-major axis of the ellipsoid of revolution.So the error of the computation of the elliptic integral should be smaller than 1.10 -11 to get the arc length with an accuracy of 1 mm.This is the reason why for all further computations a threshold ε for the recursive processes of ε < 1.10 -15 is introduced.
Also the transformation between isometric and geographic coordinates could be computed with elliptic integrals of the second kind (HECK, 1995, p. 262).Thus, their solution is still of interest and the transformation of isometric coordinates is their main area of application.
Nevertheless there are other geodetic problems which lead to elliptic integrals.For example the so-called indirect geodetic problem.In this context the length of the geodesic s can be calculated by:

THE COMPUTATION OF ELLIPTIC INTEGRALS OF THE FIRST KIND
First of all have a look at the definite elliptic integral of first kind given by: (5) The substitution of φ by ω leads to 1 : 1 In this section as well as in section 3 all explanations concerning the derivation of the Landentransformation are based on (SCHLÖMILCH, 1866).
Bol. Ciênc.Geod., sec.Artigos, Curitiba, v. 17, n o 1, p.03-22, jan-mar, 2011.and the derivate: The division of ( 7) and ( 8) combined with the substitution of cos(2ω) leads to ( 9) Now the above shown elliptic integral of first kind can be written as: Both sides of the formula are elliptic integrals of the first kind.But whereas the left integral is original the right integral is of different modulus and amplitude.This is abbreviated by: Bol. Ciênc.Geod., sec.Artigos, Curitiba, v. 17, n o 1, p.03-22, jan-mar, 2011.11) the modulus of the left integral is smaller whereas the amplitude is bigger.This can be repeated and as a result increasing module and decreasing amplitudes are obtained.Written as a chain of elliptic integrals of the first kind with increasing module we get (κ = κ 0 ) and (φ = φ 0 ): (18)

L
With this equation we gained a practical formula to compute the values of an elliptical integral of first kind.
But there is still another way to read equation ( 11).It follows directly: ( which means that with a given q and ω an elliptical integral with a smaller modulus is found on the right.The relations between p, φ and q, ω are: ( Written as a suitable formula the above could be combined to (with κ 0 = q and φ 0 = ω) κ ω Geod., sec.Artigos, Curitiba, v. 17, n o 1, p.03-22, jan-mar, 2011.

Solution with Increasing Modulus
In this chapter numerical aspects of an algorithm based on (18) are outlined.This formula could be easily implemented as recursive algorithm.At first, the convergence of this method should be examined.Therefore the module of each iteration loop is considered.The results are depicted in table 1 (the values are computed with φ = 40 °).Table (1) shows, that even for k 0 = 0.1 the algorithm needs only 5 steps to achieve the break criterion |κ n-1 -κ n | = ε < 1.10 -15 (the results in all the tables are calculated with ε = 1. 10 -15 ).This gives an idea for the convergence of the Landen-transformation.
k 0 = 0.1 k 0 = 0.5 k 0 = 0.9 k 1 0.574959574576069 0.942809041582063 0.998613997947909 k 2 0.962895683726437 0.999566630206701 0.999999759541600 k 3 0.999821325230250 0.999999976513650 0.999999999999993 k 4 0.999999996008703 1.0 1.0 k 5 1.0 Furthermore it is remarkable that for big module the break criterion of the iteration is fulfilled faster (as expected).For κ = 0.1 e. g. there are 5 steps needed, whereas for κ = 0.9 4 steps are enough.If the break criterion ε would be changed to ε = 1. 10 -14 , than for κ = 0.9 even 3 steps instead of 4 are sufficient.The following table shows the relation between the modulus and its correlated amplitude.As an example the modulus κ 0 = 0.1 is used (table 2).Whereas the modulus in column 2 increases from 0.1 to 1.0 within 5 steps, the amplitude (see column 3) decreases from 40.0 degrees to about 17 degrees.Table 3 gives an overview about the numerical results calculated by increasing module.These values could be compared with those of table 6.The implementation of the algorithm based on ( 18) is shown below.The source code is written in C. Like all recursive processes the function to compute F (k, φ) is rather compact.Further on the example is intended to show how easy equation ( 18) can be implemented.The quantities modul and amplitude in the call of the function are selfexplanatory.eps denotes the break criterion for the recrusive process and can be any arbitrary value.By changing this quantity the user can control the processing time and the accuracy of the algorithm.The quantity EllInt_PI is a symbolic constant and simply stands for the circle constant π.

Solution with Decreasing Modulus
In this chapter an iteration based on (25) is outlined.As already shown before, this formula could be used to compute module converging to 0. To calculate the different values for the module and the amplitudes the formulas (20) and ( 21) could be applied.
According to what we said before, in this chapter the results of the tables 1 -3 should now be computed with an algorithm, which works with decreasing module.
As column two in table 4 points out, this algorithm needs for κ = 0.1 only four steps to reach its break criterion, whereas for κ = 0.5 and κ = 0.9 five steps are needed.In chapter 2.1 this was contrariwise.The method using increasing module needed five steps for κ = 0.1 and for κ = 0.5 and κ = 0.9 four are sufficient.This result is exactly what is expected.Because the goal of the algorithm working with increasing module is to reach a modulus of 1, whereas the other one should reach a modulus of 0.
Bol. Ciênc.Geod., sec.Artigos, Curitiba, v. 17, n o 1, p.03-22, jan-mar, 2011.Now have a look at table 5, which depicts the decreasing module and the increasing amplitudes.In the last row of column three the angle is about two times bigger than that one of the last but one iteration.Besides the angle of this iteration is almost twice as big as the angle computed in the step before.Of cause this is not incidental.The connection between φ n and φ n+1 is given by ( 26) Hence, for small module the angle will double from one iteration to the next.The derivation of algorithms to compute elliptic integrals of the first and...
Finally the results of table 3 and table 6 should be compared.In seven out of twenty cases the results are identical.Nine times the difference is about 1.10 -15 .In three cases the difference is 2.10 -15 and only once it is 11.10 -15 .This means that the results confirm each other.

Solution with Decreasing Modulus and Series Expansion
Using decreasing module there is still another way to find an algorithm based on formula (22).As table 4 shows, the module decrease rapidly.Even for κ = 0.9 the modulus of step three is less than 5.10 -4 .This seems to be small enough to try a series expansion.
The term (1 -x) -½ (|x| < 1 but it works better if |x| << 1 could be written as ( 28) Thus ( 22) will become to (30) This formula could be also implemented as a recurrence.The break criterion would then be a modulus k which is smaller than 1.10 -3 .An approximation leads us to the following results.If κ = 1.10 -3 (worst case scenario) we get for the 4 th term R in equation ( 30). (31) The effect of this term is marginal.So it is sufficient to take the first three terms of the series in formula (30).Table 7 depicts the results computed with the first three terms.

Comparison of the Results
Obviously it is difficult to decide which of the results is right.Therefore all computations are compared with those gained by Maple V, which are treated as a reference.The results are illustrated in table 8. Comparing table 8 with the tables 3 -7 shows that only two out of twenty values have a difference of more than 1.10 -15 .One of these differences is 2.10 -15 and the other is 8.10 -15 , computed for F (0.999, 90 °) with increasing module.At first sight, the deviation d = 8.10 -15 seems to be a contradiction to the kind of method which is applied.Because increasing module should lead to rather good results, as the number of iterations should be small.If we look at the module, which are κ 0 = 0.999, κ 1 = 0.999999874874898, κ 2 = 0.999999999999998 and κ 3 = 1.0, we realize, that which means that the break criterion is failed only by 1.10 -15 .I. e. even the last but one iteration κ 2 almost fulfills the break criterion.Due to numerical aspects the next iteration delivers no further contribution to the result.So the derivation between the reference and the computed solution is relatively big.Nevertheless the Landentransformation is a useful method to compute elliptic integrals of the first kind.Even the biggest difference of 8.10 -15 between the two computed results easily fulfills the geodetic requirements mentioned in section 1.

THE COMPUTATION OF ELLIPTIC INTEGRALS OF SECOND KIND
Also for the elliptic integral of the second kind (see 2) there exist explicit solutions for two specific modules.The first is κ = 1 and the second is κ = 0.
Bol. Ciênc.Geod., sec.Artigos, Curitiba, v. 17, n o 1, p.03-22, jan-mar, 2011.(40) or (41) Changing the indexes to κ = κ 0 and q = κ 1 as well as φ = φ 0 and ω = φ 1 we obtain In this formula we have on the left the original elliptic integral of the second kind with its modulus κ 0 and its amplitude φ 0 , whereas on the right we obtain an elliptic integral of the second kind with a smaller modulus κ 1 combined with an elliptic integral of first kind whereas the amplitude is modified in both cases to φ 1 .

Solution with Increasing Modulus
In fact there are two different ways to read formula ( 42).The first one is that by subsequent conversion of the term E(κ n , φ n ) until (43) and thus is obtained.Obviously this can be implemented as a recursive algorithm too.It must be considered that this formula encompasses two recursive processes because both, the elliptic integral of second kind and that of the first kind have to be calculated that way.But for the computation of the latter we already found various solutions in section 2.
In this formula the elliptic integral of the first kind has the same modulus and the same amplitude as the elliptic integral of the second kind.This conversion could be repeated in the same way as we have done before.As a result we get The above shown equation seems to be not so time consuming than (46) because the elliptic integral of the first kind has to be computed only once.Due to numerical aspects it is not possible to implement it as a recursive algorithm.In step n the modulus is approximately 1, because this is the break criterion.But now, beginning with this modulus, the decreasing chain of the moduli has to be computed.This of cause is also a recursive process (at least in theory).In fact we get as a result of Bol. Ciênc.Geod., sec.Artigos, Curitiba, v. 17, n o 1, p.03-22, jan-mar, 2011.
As a consequence the recursive process will have no end.This problem could only be solved with a quasi -recursive algorithm, which puts the computed moduli of step i in a stack from which the product of the moduli is derived (see mark '!' in equation ( 46)).If done so, the numerical values are exactly the same as the results computed with (44) (see table 9 together with table 10).As the results are nearly the same, the algorithm based non ( 46) is about 35% faster than this based on (44) 2 .However the implementation of ( 46) is much more complicated, because it is not strictly recursive, as already mentioned.
An implementation of an algorithm based on ( 44) is depicted below.F_inc(n_mod,n_amp,eps)-modul*sin(amplitude); } } The example shows further that the recursive process invoked by the call of E_inc 3 comprises a second one invoked by F_inc (see section 2.1).I. e. the second is nested within the first.Thus the algorithm is extremely compact.

Solution with Decreasing Modulus
An algorithm to compute E(p, φ) with decreasing module based on ( 42) is obtained by (48)

E
In this formula we assume that E(κ 1 , φ 1 ) has to be computed and E(κ 0 , φ 0 ) is known.Therefore it makes sense to change the indexes and we get (49) The terms κ 1 and φ 1 are computed with (20) and ( 21) respectively.So the modulus κ 1 is smaller than κ 0 whereas φ 1 is larger than φ 0 .E(κ 1, φ 1 ) could be treated in the same way as E(κ 0, φ 0 ) until we get κ < ε, which means in practice κ ≈ 0. With (50) a break criterion is given.Thus a recursive process to compute E(p, φ) is found.The numerical results are shown in table 11.

Solution with Decreasing Modulus and Series Expansion
In the same way we did in section 2.3, we can combine the Landentransformation with the series expansion.In this context this seems to be even more advantageous because the computation of an elliptic integral of the second kind is even more time-consuming.
The series expansion of As above mentioned formula (51) will deliver the best results if |x| << 1and thus κ << 1. Hence we decrease the modulus by Landen-transformation to apply the series expansion to the transformed integral subsequently.
If we look at the 4 th term of the series in formula (53) and integrate it assuming k = 0.001 and the maximum amplitude, the result r is r < 2.7.10 -26 .So it is sufficient to take only the first three terms of the series.Thus the fourth and all further terms are omitted.So formula (49) in combination with (53) yields   It makes sense to compute the elliptic integral of first kind in equation ( 54) also by decreasing module in combination with the series expansion, because in each step of the recursion the modulus decreases and thus the break criterion is reached faster.

Comparison of theRresults
The comparison between the tables 9 -12 shows that most of the values are identical with those of table 13, which is treated as a reference.The biggest difference concerns E(0.999, 90 °), which leads to 2.10 -15 between table 12 and 13.Beside the accuracy there is still another criterion which should not be neglected.This criterion is the processing time.All algorithms of section 3 comprise two recursive processes.So it should be considered to take the right combination of both, to get the best result.Whereas in this case best means fastest.E. g. it seems to be different if we solve E(κ, φ) by decreasing module, then it makes sense to apply decreasing module too, in oder to compute the appropriate elliptic integral of the first kind.
This can easily be proved using E(0.5, 90°) as an example.The algorithm based on 54 in combination with that based on 30 is about 25% faster than that based on 54 together with 25.The same test applied on E(0.1, 90 °) delivers still an Bol.Ciênc.Geod., sec.Artigos, Curitiba, v. 17, n o 1, p. 03-22, jan-mar, 2011.advantage of 20% between the two approaches.So the decision between the different methods is mainly influenced by the processing time, because the accuracy is the same.

CONCLUSIONS
As an example the result of the arc length for φ = 48 ° is computed.With (3) we get S = 5317885.233m (φ denotes the amplitude; further on the parameters of the ellipsoid of Bessel are chosen) which is exactly the same result as published by (SCHÖDLBAUER, 1981).The result remains the same regardless which algorithm is applied.Even the combination of decreasing and increasing module (e.g.formula 49 together with 18 or 42 together with 25) leads to the same result.Therefore it is almost needless to explain that as a consequence also the transformation from geographical coordinates to conformal coordinates (e. g. of type Gauss-Krueger or UTM) yields the same results.
The advantages of the algorithms based on the Landen-transformation are obvious.They are flexible which means that they can be applied to different problems and the accuracy can be fitted to the requirements.Further on they are simple to implement as a recursive process and thus the test of the implementation is easy as well.As a consequence the presented solutions are an easy to use alternative to the known approaches.
have done in section 3.1.Now the elliptic integral of the first kind needs to be computed only once.
depicts the last term of the recursion.At this stage the modulus κ is κ < 10 -3 .
Thus we got an easy way to use formula to compute ω.With

Table 2 -
The moduli and the related amplitudes for k 0 = 0.1.

Table 3 -
F (k i , φ i ) computed by increasing moduli.

Table 4 -
The changing moduli for three different iterations.

Table 5 -
The moduli and the appropriate amplitudes for k 0 = 0.9.

Table 7 -
F (k i , φ i ) computed with decreasing moduli and series expansion.

Table 8 -
F (k i , φ i ) computed with Maple V.