(41/42) AGM - Arithmetic-Geometric Mean

[Albert Chan]

Thanks, Werner

Great idea on testing convergence of arithmetic mean.
I think this always work (no infinite loops).

My argument is like this (AM, GM for arithmetic and geometric means)

We know AM ≥ GM, but with rounding error, it is possible the rule fail.
At the final stage of AGM, with rounding error, we may have AM ≥ GM - 1ULP
(gap cannot be more than 1 ULP, otherwise next sequence will reduce it)

Example, with IEEE double, final stage of AGM:

AGM(a,b) = AGM((a+b)/2, √(ab)) = ...

AGM(1,4) ⇒ AGM(2.2430285802876027, 2.2430285802876031) ; AM = GM - 1 ULP
AGM(1,5) ⇒ AGM(2.6040081905309407, 2.6040081905309402) ; AM = GM + 1 ULP
AGM(1,6) ⇒ AGM(2.9513287423905727, 2.9513287423905727) ; AM = GM

This shows that testing convergence AM = GM will not work.

Note that, even with rounding error, AM sequence is decreasing.
(We required GM ≥ AM + 2 ULP to reverse this trend)

If we assume non-zero AGM arguments, AM sequence always converge.
This applied to intermediates too. Example, AGM(-1,1) should not converge (*)

Above (max of ± 1 ULP) may apply to complex numbers too.
Again, using IEEE double for real and imaginary parts.

AGM(1+1j, 1+5j):
AM ⇒ (1.1429851941908882 + 2.680170008949454j)
GM ⇒ (1.142985194190888 + 2.6801700089494545j) ; AM = GM + (1,-1) ULP

This suggest your setup also work for complex numbers

(*) Actually, AGM(-1,1) converged for IEEE double, when AM underflow to 0.0
I wonder testing for GM convergence may be better, quiting much earlier.