Who discovered the "AGM-like" approximation of arctangent?

Question

MathWorld (43-47) gives the following algorithm to compute the arctangent:

To find $\tan^{-1}x$ numerically, the following arithmetic-geometric mean-like algorithm can be used. Let

$$a_0 = (1+x^2)^{-1/2}$$

$$b_0 = 1$$

Then compute

$$a_{i+1} = (a_i+b_i)/2$$

$$b_{i+1} = \sqrt{a_{i+1}b_i}$$

and the inverse tangent is given by

$$\tan^{-1}x=\lim_{n\to\infty}\frac{x}{\sqrt{1+x^2}a_n} $$

(Acton 1990).

NB that it is not the AGM series. See the $a_{i+1}$ in the right-rand side for $b_{i+1}$.

Acton 1990 is Acton, F. S. "The Arctangent." In Numerical Methods that Work, upd. and rev. Washington, DC: Math. Assoc. Amer., pp. 6-10, 1990; originally published in 1970.

However, a much earlier document exists demonstrating a variation of the algorithm. Below is a page from Edsgar Dijkstra's notes on his Algol-60 compiler, describing the implementation of the standard library function arctan. There $a_0 = 1, b_0 = \sqrt{1+x^2}$, and the expression under the limit would be simply $x/b_n$.

For the purpose of my question, disregard the weighted sum. Instead of iterating until the desired precision of $10^{-12}$, the program runs 4 iterations, for which they have found weights of $b_0 ... b_4$ plus a constant term achieving the closest approximation.

My question is, if the "AGM-like" algorithm to compute the arctangent was already known by 1961, who is its true discoverer?

Answer 1

The AGM-like iteration used here is sometimes called Borchardt's algorithm, as it appeared in a publication by the German mathematician Carl Wilhelm Borchardt:

C.W. Borchardt, "Sur deux algorithmes analogues à celui de la moyenne arithmético-géométrique de deux éléments." In L. Cremona and E. Beltrami (eds.), In Memoriam Dominici Chelini, Milan 1881, pp. 206-212.

Translated title: "On two algorithms analogous to that of the arithmetic-geometric average of two elements". Borchardt's findings already appeared in correspondence between J.F. Pfaff and Gauss dating to 1800: Gauss, Werke, Vol. 10, 1917, p. 234. The iterative formula itself also appeared in part 1 of Jacques Schwab, Élémens de géométrie. Nancy: C.-J. Hissette 1813, p. 104.

There is no mention of $\arctan$ in Borchardt's publication, instead it is stated (page 210) that for $m \lt n$, the iteration converges to $$\omega = \frac{\sqrt{n^2-m^2}}{\arccos\left(\frac{m}{n}\right)}$$ By straightforward substitution, using $m=1, n=\sqrt{x^{2}+1}$, one arrives at $\omega = \frac{x}{\arctan x}$. This iteration converges quite slowly and is therefore of little practical use unless one can find a way to accelerate the convergence.

At first I was not able to determine what fourth-order convergence acceleration scheme is used in the note reproduced above. I found that the British mathematician Peter Wynn worked at the MC from 1960 to 1964 where he had some involvement with libraries for the Electrologica X1. In his career he was mostly focused on approximation theory. Several of his publications from the 1950s and 1960s deal with convergence acceleration. This makes it likely that he devised the convergence acceleration scheme used for the note reproduced above.

B. C. Carlson, "An Algorithm for Computing Logarithms and Arctangents." Mathematics of Computation, Vol. 26, No. 118, April 1972, pp. 543-549, also computes $\arctan$ with Borchardt's algorithm, and provides convergence acceleration using a recursion scheme that basically interpolates between the $a_i$.

Expanding the recursion in Carlson's scheme gives:

$$ d(5,5) = \frac{1073741824a_{5} -357564416a_{4}+23744512a_{3}-371008a_{2}+1364a_{1}-a_{0}}{739552275}$$

Based on Borchardt's algorithm one can trivially express the $a_{i}$ in terms of the $b_{i}$ as follows: $a_1=(1+b_{0})/2$, $a_{2}=(1+b_{0}+2b_{1})/4$, $ \ \ldots \ $, $a_{5}(1+b_{0}+2b_{1}+4b_{2}+8b_{3}+16b_{4})/32$.

Plugging this into the interpolation formula for $d(5,5)$ above results in:

$$\frac{536870912b_{4}+89653248b_{3}+56698880b_{2}+28163936b_{1}+14082650b_{0}+14082649}{739552275}$$

These are precisely the interpolation weights used in the handwritten note from the question. Whoever devised this interpolation scheme thus anticipated Carlson's work a decade prior to its publication. Interpolation utilizing the $b_{i}$ instead of the $a_{i}$ may have an advantage in fixed-point computation as all coefficients are strictly less than unity and thus pure fractions.

I note that the coefficients used in the convergence acceleration scheme from the note in the question are not "optimal". Using these, the maximum relative error of the $\arctan$ approximation is about $6.6 \cdot 10^{-13}$, while the maximum relative error can be reduced to $\lt 1.92 \cdot 10^{-13}$ by using the following set of coefficients:

$0.725940450929954960$, $0.121226383896385010$, $0.076666493926230372,$
$0.038082414121090931$, $0.019042129238796116$, $0.019042127887638049$

Answer 2

FYI, AGM like formula is equivalent to arctan "half-angle" formula.

$x = \tan(θ) = \frac{2\,\tan(\frac{θ}{2})} { 1 \;-\; \tan^2(\frac{θ}{2}) }$

$ \tan(\frac{θ}{2}) = \frac{x}{\sqrt{x^2+1}\;+1}$

$\displaystyle θ = \tan^{-1}(x) = 2\;\tan^{-1} \left(\frac{x}{\sqrt{x^2+1}\;+\;1} \right)$

If we estimate RHS, $\;\tan^{-1}(ε) ≈ ε\;$, it is equivalent to AGM-like iterations.
We can prove this by induction. Here, I just shown they are numerically equivalent.

lua> k, x = 1, 1                -- k * atan(x) = pi/4
lua> a, b = 1/sqrt(1+x*x), 1
lua> c = x*a
lua> for i=1, 10 do
:    a=(a+b)/2; b=sqrt(a*b)     -- AGM like iteration
:    k=k+k; x=x/(sqrt(x*x+1)+1) -- half-angle formula
:    print(c/a, k*x)            -- estimate atan(x)≈x
:    end

0.8284271247461901      0.8284271247461902
0.795649469518632       0.7956494695186321
0.787931226857314       0.787931226857314
0.7860295963114761      0.7860295963114762
0.7855559074856141      0.7855559074856143
0.7854375922922415      0.7854375922922417
0.7854080201757955      0.7854080201757957
0.7854006275642023      0.7854006275642026
0.7853987794373972      0.7853987794373977
0.7853983174073268      0.7853983174073271