function m = AGM(a,b)

m = min(a,b);
M = max(a,b);

while max((M-m)./m) > eps
    m0 = sqrt(m.*M);
    M = (m+M)/2;
    m = m0;
end

return