/* Calculate the arithmethic-geometric mean of two positive
 * numbers x and y.
 * Result will have d digits after the decimal point.
 */
define m(x, y, d) {
    auto a, g, o

    o = scale
    scale = d
    d = 1 / 10 ^ d

    a = (x + y) / 2
    g = sqrt(x * y)
    while ((a - g) > d) {
        x = (a + g) / 2
        g = sqrt(a * g)
        a = x
    }

    scale = o
    return(a)
}

scale = 20
m(1, 1 / sqrt(2), 20)