11

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?

enter image description here

kipf
  • 2,531
Leo B.
  • 253
  • 2
  • 7
  • @njuffa Or it could have been his co-author Jaap Zonneveld who was more into numerical math. That curious result alone should have been worth an article, if one of them had indeed discovered it. It is peculiar why they didn't simply do +/- pi/2 -/+ arctan(1/x) for |x| > 1, and didn't approximate arctan(x) on [0,1] with a polynomial which would have been quite a few times faster. – Leo B. Aug 10 '24 at 07:44
  • 2
    I think the advantage of the approach proposed here is that it leads to a very compact solution (as measured by total words required for code +data). BTW, it seems the same algorithm was used in the English Electric KDF9, at least I found the same constants used: http://www.findlayw.plus.com/KDF9/Arctan.html – njuffa Aug 10 '24 at 08:13
  • 2
    If not Dijkstra, Peter Wynn seems the more likely source. He worked at CWI between 1960 and 1964, and authored "An arsenal of Algol procedures for complex arithmetic", BIT, Vol. 2, No. 4, 1962 (scan). He could also have been the link to the KDF9 implementation of $\arctan$? – njuffa Aug 10 '24 at 08:40
  • The handwriting in the note does not match Dijkstra's or Wynn's. ALD refers to the double-precision arithmetic complex. In the 1950s and 1960s Wynn published on convergence acceleration, which appears to be the purpose of the denominator here. The $\arctan$ implementation used in the Electrologica X8 is of the conventional type: F.J.M. Barning, Standaardfuncties in X8 Assembler-code ELAN. Report, Mathematical Center, Computing Department, Nov. 1965 (scan) – njuffa Aug 11 '24 at 19:55
  • 1
    From the Kruseman Aretz's paper about the compiler: While Dr. E.W. Dijkstra and J.A. Zonneveld were developing the compiler Miss M.J.H. Ro¨mgens and Miss S.J. Christen started work on the organisational and arithmetic subroutines The handwriting, giving its neatness, is likely Miss Ro¨mgens's or Miss Christen's, then. – Leo B. Aug 11 '24 at 20:10
  • B. C. Carlson, "An Algorithm for Computing Logarithms and Arctangents." Mathematics of Computation, Vol. 26, No. 118, April 1972, pp. 543-549, (scan), computes $\arctan$ with the same AGM-like iteration used here and refers to it as Borchardt's algorithm. He notes its slow convergence and proposes a scheme to accelerate the convergence, however his approach is very different from the convergence acceleration scheme used here. – njuffa Aug 12 '24 at 06:06
  • 1
    The relevant publication by Borchardt is: 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 (scan from Borchardt's collected works). – njuffa Aug 12 '24 at 06:29

2 Answers2

5

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.

I have not been 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, but proposes an entirely different way to accelerate the convergence.

I note that the weight coefficients used in the 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 $< 1.8 \cdot 10^{-13}$ when using the following set of coefficients:

$0.725940450930031340$, $0.121226383896505420$, $0.076666493926177651$,
$0.038082414121027114$, $0.019042129238739429$, $0.019042127887694091$

njuffa
  • 2,168
  • 1
  • 20
  • 21
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
albert chan
  • 2,204