Computing logs with the AGM

In the previous post I said that Jonathan and Peter Borwein figured out how to use the rapid convergence of the AGM to compute various functions, including logarithms. This post will show how to compute logarithms using the AGM.

First I need to define the Jacobi theta functions

Note that if q is very small, the series above converge very quickly. We will pick q so small that θ₂(q) and θ₃(q) are trivial to compute with sufficient accuracy.

Suppose we want to compute the natural log of 10 to a thousand digits. We can use the following equation from [1].

We want q to be small, so we want 1/q to be large. We set q = 10^–n to compute log 10ⁿ and divide the result by n.

We can show from the definition of the theta functions that

Since we want 1000 digit accuracy, we can set q = 10^-250. Then θ₂²(q⁴) = 4×10^-500 and θ₃²(q⁴) = 1 to about 2000 digits of accuracy.

We can calculate the result we want with just 20 iterations of the AGM as the following bc code shows.

    define agm(a, b, n) {
        auto i, c, d
        for (i = 1; i <= n; i++) {
            c = (a + b)/2
            d = sqrt(a*b)
            a = c
            b = d
        }
        return a
    }

    scale=2000
    # a(1) = arctan(1) = pi/4
    x = a(1)/(agm(4*10^-500, 1, 20)*250) - l(10) # error
    l(x)/l(10) # log_10(error)

This returns -999.03…, so we were basically able to get 1000 digits, though the last digit or two might be dodgy. Let’s try again with a smaller value of q, setting 10^-300.

    x = a(1)/(agm(4*10^-600,1, 20)*300) - l(10); l(x)/l(10)

This returns -1199.03… and so we have about 1200 digits. We were using 2000 digit precision in our calculations, but our approximation θ₃²(q⁴) = 1 wasn’t good to quite 2000 digits. Using a smaller value of q fixed that.

Related posts

• Computing π with the AGM and mpmath
• Three surprises with bc

[1] The Borwein brothers, Pi and the AGM. arXiv