procedure main(A)
    a := real(A[1]) | 1.0
    g := real(A[2]) | (1 / 2^0.5)
    epsilon := real(A[3])
    write("agm(",a,",",g,") = ",agm(a,g,epsilon))
end

procedure agm(an, gn, e)
    /e := 1e-15
    while abs(an-gn) > e do {
       ap := (an+gn)/2.0
       gn := (an*gn)^0.5
       an := ap
       }
    return an
end