function agm(atom a, atom g, atom tolerance=1.0e-15)
    while abs(a-g)>tolerance do
        {a,g} = {(a + g)/2,sqrt(a*g)}
        printf(1,"%0.15g\n",a)
    end while
    return a
end function
?agm(1,1/sqrt(2))   -- (rounds to 10 d.p.)