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 2(q) and 3(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 10n 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 22(q4) = 4*10-500 and 32(q4) = 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 32(q4) = 1 wasn't good to quite 2000 digits. Using a smaller value of q fixed that.
Related posts[1] The Borwein brothers, Pi and the AGM. arXiv
The post Computing logs with the AGM first appeared on John D. Cook.