Finding the coefficients (weights) for a convergence acceleration scheme

Question

This is a followup question to my answer for a previous question, where I identified the $\arctan$ algorithm used in an 1961 Algol compiler for the Electrologica X1 (subsequently also used in an Algol compiler for the English Electric KDF9 in 1963) as Borchardt's algorithm combined with an unspecified convergence acceleration scheme. Borchardt's algorithm uses an AGM-like iteration

$$a_{i+1} = \frac{a_{i}+b_{i}}{2}$$ $$b_{i+1} = \sqrt{a_{i+1} b_{i}}$$

When $0 \le a_{0} \lt b_{0}$, and with $\lim_{i\to\infty} a_{i} = A, \lim_{i\to\infty} b_{i} = B$, the iteration converges to

$$A = B = \frac{\sqrt{b_{0}^{2} - a_{0}^{2}}}{\arccos \left(\frac{a_{0}}{b_{0}}\right)}$$

Setting $a_{0}=1$ and $b_{0}=\sqrt{x^{2}+1}$ results in $$A = B = \frac{x}{\arctan (x)}$$

For practical finite-precision computation of $\arctan$ one might set $a_{0} = 1$, $t_{0} = x^{2}+1$, and $b_{0} = \sqrt{t_{0}}$, then iterate

$$a_{i+1} = \frac{a_{i}+b_{i}}{2}$$ $$t_{i+1} = \frac{a_{i}b_{i}+t_{i}}{2}$$ $$b_{i+1} = \sqrt{t_{i+1}}$$

and after $n$ iterations, compute $\arctan (x)$ as $\frac{x}{b_{n}}$. With either variant of the iteration, the sequence terms are represented by

$$a_{i} = \frac{x}{2^{i}\tan \left(\arctan \left(\frac{x}{2^{i}} \right)\right)} = \cos\left(2^{-i}\arctan(x)\right)\ b_{0}\prod_{k=1}^{i} \cos\left(2^{-k}\arctan(x)\right)$$ and $$b_{i} = \frac{x}{2^{i}\sin \left(\arctan \left(\frac{x}{2^{i}} \right)\right)} = b_{0}\prod_{k=1}^{i} \cos\left(2^{-k}\arctan(x)\right)$$

Denoting the ratio of corresponding sequence terms as $r_{i} = \frac{a_{i}}{b_{i}} = \cos\left(2^{-i}\arctan(x)\right)$, we have

$$r_{i+1} = \sqrt{\frac{1}{2}\left(1+r_{i}\right)}$$

Borchardt's algorithm converges rather slowly: the relative error of the sequence terms vis-a-vis the limit is reduced by a factor of 4 in each iteration. This makes it unsuitable for practical computation.

The convergence acceleration used for the algorithm in the Electrologica X1 and the English Electric KDF9 consists of a simple weighted sum applied to the $b$ sequence. After performing $n$ iterations of Borchardt's algorithm it computes

$$ s_{n} = c_{0} + \sum_{i = 0}^{n} b_{i}c_{i+1}$$

using the coefficients (weights) $c_{i} > 0$ and $\sum_{i=0}^{n+1} c_{i} \approx 1$. This scheme proves remarkably effective and provides significant convergence acceleration for several small $n$ that I tried, as can be seen from the following data. The data suggests quadratic convergence of the $s_{n}$. Here, I computed the $c_{i}$ numerically, first solving a system of equations, then following up with a heuristic search to provide near-minimax approximations. $$ \begin{array}{|c|c|c|l|} \hline n & \text{maxrelerr } b_{n} & \text{maxrelerr } s_{n} & \text{}c_{n+1}, \ldots, c_{0} \\ \hline 0 &3.64\times10^{-1} & 6.33 \times 10^{-3} & 0.640673255622243090, 0.362501346012709110 \\ \hline 1 &9.97\times10^{-2} & 3.42 \times 10^{-5} & 0.696171720382026130, 0.144352031074684220, 0.159461397445679970 \\ \hline 2 & 2.55\times 10^{-2} & 1.98 \times 10^{-7} & 0.721874544769588230, 0.124637158070154310, 0.076900841592368346, 0.076587554224942539 \\ \hline 3 & 6.42\times 10^{-3} & 7.03 \times 10^{-10} & 0.725231395234937290, 0.121816205684205090, 0.076710658890471975, \ 0.038120187333781988, 0.038121552505571692 \\ \hline 4 & 1.61 \times 10^{-3} & 1.92 \times 10^{-13} & 0.725940450929954300, 0.121226383896385460, 0.076666493926229748, 0.038082414121091576, 0.019042129238796179, 0.019042127887638319 \\ \hline 5 & 4.02 \times 10^{-4} & 8.63 \times 10^{-18} & 0.72611772576565280359, 0.12107871257274047205, 0.07665561239099945804, 0.03807299187837815249, 0.01903747959739412644, 0.00951873889725269016, 0.00951873889758229282 \\ \hline \end{array} $$

I have demonstrated the practical utility of this approach for the low-precision computation of $\arctan$ on modern computing platforms where fast (but possible approximate) square root computation is provided by the hardware.

It seems to me that it may also be useful for $\arctan$ computations in arbitrary-precision computation, and I would like to explore this further. However, this is hampered by the lack of an efficient way to compute the coefficients used in this sequence acceleration scheme. Is there a closed-form formula or a recursion formula for rapidly computing these coefficients? From the data so far, the coefficients appear to be following a pattern, so I suspect such methods could exist.

Because I suspected, for a number of reasons, the involvement of the British mathematician Peter Wynn in the creation of the $\arctan$ algorithm for the Electrologica X1, I first looked into Wynn's $\epsilon$-algorithm. While it provides convergence acceleration, the acceleration is not nearly as pronounced as with the simple weighted sum. I also looked at Wynn's research notes and unpublished manuscripts, but found nothing relevant.

I looked at publications on Borchardt's algorithm and related iteration schemes, both those without ^[1]-[4] and those with ^[5]-[11] an acceleration methods, but did not find anything that matched closely. I tried the Romberg method, and while it provides acceleration, it does so much less than the weighted sum method. The closest match is Carlson's extrapolation scheme ^[9], which, however, works with the terms $a_{i}$ rather than the $b$ sequence and also does not provide the same amount of acceleration. I also studied publications specifically on convergence acceleration and extrapolation methods ^[12]-[16] but could not find anything applicable to the problem at hand. I am not a mathematician and have not worked with convergence acceleration or extrapolation methods before, so it is entirely possible that I overlooked some useful connection to the material in these publications.

I used the inverse symbolic calculator to determine whether the coefficients I found numerically are near any "interesting" values, which does not seem to be the case. The coefficients used in the Electrologica X1 for $n=4$ are all fractions with a denominator of $1479104550$; that did not provide a clue either. I also rationalized the coefficients I found numerically into various fractions and looked for numerator and denominator sequences in the OEIS without success.

In terms of further brainstorming, I have considered equating series expansions for the $b$ and $s$ sequences and possible connections to Padé approximants or Chebyshev expansions, but I do not have the mathematical wherewithal to effectively follow through on those ideas.

---

[1] Bille Chandler Carlson, "Algorithms involving arithmetic and geometric means." The American Mathematical Monthly, Vol. 78, No. 5, 1971, pp. 496-505

[2] Francesco Giacomo Tricomi, "Sugli algoritmi iterativi nell'analisi numerica." In Metodi valutativi nella fisica-mathematica, Rome, December 15-19, 1972 (Quaderno N. 217, Accademia Nationale dei Lincei 1975), pp. 105-117

[3] John Todd, "The algorithms of Gauss, Borchardt, and Carlson." In John Todd, Basic Numerical Mathematics, Vol. 1: Numerical Analysis, Basel: Birkhäuser 1979, pp. 13-23

[4] Isaac Jacob Schoenberg, "On the arithmetic-geometric mean and similar iterative algorithms." In Isaac Jacob Schoenberg, Mathematical Time Exposures, Mathematical Association of America 1982, pp. 149-167

[5] R. M. Milne, "Extension of Huygens' approximation to a circular arc." The Mathematical Gazette, Vol. 2, No. 40, Jul. 1903, pp. 309-311

[6] Karl Kommerell, "Der Begriff des Grenzwertes in der Elementarmathematik." In Karl Kommerell, Das Grenzgebiet der elementaren und höheren Mathematik, Leipzig: Koehler Verlag 1936, pp. 1-102

[7] Heinz Rutishauser, "Ausdehnung des Rombergschen Prinzips." Numerische Mathematik, Vol. 5, 1963, pp. 48-54

[8] Siegfried Filippi, "Die Berechnung einiger elementarer transzendenter Funktionen mit Hilfe des Richardson-Algorithmus." Computing, Vol. 1, 1966, pp. 127-132

[9] Bille Chandler Carlson, "An algorithm for computing logarithms and arctangents." Mathematics of Computation, Vol. 26, No. 118, Apr. 1972, pp. 543-549

[10] G. M. Phillips, "Archimedes the Numerical Analyst." The American Mathematical Monthly, Vol. 88, No. 3, Mar. 1981, pp. 165-169

[11] George M. Phillips, "Archimedes and the complex plane." The American Mathematical Monthly, Vol. 91, No. 2, Feb. 1984, pp. 108-114

[12] Hermann Engels, "Zur Geschichte der Richardson-Extrapolation." Historia Mathematica, Vol. 6, 1979, pp. 280-293

[13] Jacques Dutka, "Richardson extrapolation and Romberg integration." *Historia Mathematica, Vol. 11, 1984, pp. 3-21

[14] Guido Walz, "The History of Extrapolation Methods in Numerical Analysis." Report No. 130, University of Mannheim, Oct. 1991, 11 pp.

[15] Claude Brezinski, "Convergence acceleration during the 20th century." Journal of Computational and Applied Mathematics, Vol. 122, No. 1-2, 2000, pp. 1-21

[16] Avram Sidi, Practical Extrapolation Method. Cambridge University Press 2003

Answer 1

TL;DR The convergence acceleration scheme used in the Electrologica X1's Algol compiler from 1961 is in fact identical to the one used in Carlson's algorithm. Carlson ^[1] expresses it in terms of the $a_{i}$ sequence, whereas it is expressed in terms of the $b_{i}$ sequence in the Algol compiler. It is remarkable that this work preceded Carlson's publication by a decade.

I took the "scenic route" to arrive at the above conclusion by working with numerically optimized coefficients. I figured that by comparing numerically optimized coefficients based on Carlson's algorithm and comparing the results with the data from the question, I might gain additional insights.

For consistency with the question, I will refer to the sequence of arithmetic means in Borchardt's algorithm by $a_{i}$ and the sequence of geometric means in Borchardt's algorithm by $b_{i}$, whereas Carlson uses $a_{i}$ and $g_{i}$. For the arctangent algorithm, Carlson specifies a recursively defined interpolation scheme as follows: $$a_{0}=1$$ $$d(0,n) = a_{n}, \ \ \ \ \ n= 0,1, 2, \ldots$$ $$d(k,n) =\frac{d(k-1,n) - 2^{-2k}d(k-1,n-1)}{1-2^{-2k}} \ \ \ \ \ k=1,2, \ldots, n$$ If we expand the recursion completely for increasing orders $n$ and substitute $a_{0}$ with $1$, we have $$d(1,1)=\frac{4a_{1}-1}{3}$$ $$d(2,2)=\frac{64a_{2}-20a_{1}+1}{45}$$ $$d(3,3)=\frac{4096a_{3}-1344a_{2}+84a_{1}-1}{2835}$$ $$d(4,4)=\frac{1048576a_{4}-348160a_{3}+22848a_{2}-340a_{1}+1} {722925}$$ $$d(5,5)=\frac{1073741824a_{5}-357564416a_{4}+23744512a_{3}-371008a_{2}+1364a_{1}-1}{739552275}$$ $$d(6,6)=\frac{4398046511104a_{6} - 1465657589760a_{5} + 97615085568a_{4} - 1543393280a_{3} + 5957952a_{2} - 5460a_{1} + 1}{3028466566125}$$

With increasing $n$, the most significant coefficient of Carlson's interpolation scheme converges to $\prod_{k=1}^{\infty} \frac{2^{2 k}}{2^{2 k}-1} \approx 1.452353642449597$.

I note that the above formulas were already given by Kommerell ^[2], but it appears that Carlson was unaware of this work. In Kommerell's construction, the numerators of the coefficients are equal to the coefficients of the polynomial $\varphi_{n}(z) = \prod_{i=1}^{n}(2^{2i}z-1)$, while their denominator is equal to $\varphi_{n}(1)$. Example for $n=4$: $\varphi_{4}(z) = (4z-1)(16z-1)(64z-1)(256z-1) = 1048576z^{4}-348160z^{3}+22848z^{2}-340z+1$ and $\varphi_{4}(1) = 722925$.

Using the above weights as a starting point, I employed a heuristic search to arrive at numerically optimized coefficients, and found the following:

$$ \begin{array}{|c|c|l|} \hline n & \text{maxrelerr} \ d'(n,n) \\ \hline 1 & 6.33\times 10^{-3} & 1.281346511270356600, -0.278171909605450910 \\ \hline 2 & 3.42\times 10^{-5} & 1.392343436221643600, -0.407467653127648710, 0.015109366353266122 \\ \hline 3 & 1.98\times 10^{-7} & 1.443749089539337900, -0.472600228629364570, 0.029164525114581032, -0.000313287367425835 \\ \hline 4 & 7.03\times10^{-10}& 1.450462790469922800, -0.481598983866556090, 0.031605112096719234, -0.000470284222909983, 0.000001365171792833 \\ \hline 5 & 1.92\times 10^{-13} & 1.451880901859910400, -0.483487683137185490, 0.032106603956075321, -0.000501665684048061, 0.000001844356501307, -0.000000001351158067 \\ \hline 6 & 8.60\times 10^{-18} & 1.4522354515313055956, -0.4839603006201718640, 0.0322325122092584540, -0.0005096286342431526, 0.0000019673164100937, -0.0000000018028887308, 0.0000000000003295997 \\ \hline \end{array} $$ Clearly, the maximum relative error found for $d'(n,n)$ matches up exceedingly well with the maximum relative error found for $s_{n-1}$ from the question, suggesting a direct connection. Based on the definition of Borchardt's algorithm for the arctangent, we have $a_{i} = 1 + \frac{1}{2^{i}}\sum_{k=1}^{i-1}(2^{k}b_{k})$ with $i=1,2,3,\ldots$.

For example $a_{1} = (1+b_{0})/2$, $a_{2} = (1+b_{0}+2b_{1})/4$, $a_{3} = (1+b_{0}+2b_{1}+4b_{2})/8$, and so on. Applying this to $d'(2,2)= 1.392343436221643600 (1+b_{0}+2b_{1})/4 - 0.407467653127648710 (1+b_{0})/2 + 0.015109366353266122$ results in $0.6961717181108218000 b_{1} + 0.144352032491586545 b_{0} + 0.159461398844852667$, which is pretty exactly $s_{1}$ from the question, and this similarly applies to the other pairs of $d'(n,n)$ and $s_{n-1}$.

Applying the substitutions to the $d(n,n$) above, we arrive at the following interpolation formulas based on the $b_{i}$ sequence:

$$s'_{0} = \frac{2b_{0}+1}{3}$$ $$s'_{1} = \frac{32b_{1}+6b_{0}+7}{45}$$ $$s'_{2} = \frac{2048b_{2}+352b_{1}+218b_{0}+217}{2835}$$ $$s'_{3} = \frac{524288b_{3}+88064b_{2}+55456b_{1}+27558b_{0}+27559}{722925}$$ $$s'_{4} = \frac{536870912b_{4}+89653248b_{3}+56698880b_{2}+28163936b_{1}+14082650b_{0}+14082649}{739552275}$$ $$s'_{5} = \frac{2199023255552b_{5}+366682832896b_{4}+232148959232b_{3}+115302782976b_{2}+57654370464b_{1}+28827182502b_{0}+28827182503}{3028466566125}$$

$s'_{4}$ here corresponds exactly to the coefficients found in the $\arctan$ algorithm of the 1961 Algol compiler for the Electrologica X1.

---

^[1] Bille C. Carlson, "An algorithm for computing logarithms and arctangents." Mathematics of Computation, Vol. 26, No. 118, Apr. 1972, pp. 543-549

^[2] Karl Kommerell, Das Grenzgebiet der elementaren und höheren Mathematik. Leipzig, K. F. Koehler Verlag 1936, pp. 34-38