Almkvist Berndt 1988 begins with an investigation of why the agm is such an efficient algorithm, and proves that it converges quadratically. This is an efficient method to calculate 



π


{\displaystyle \pi }

.

With the same notations used in Arithmetic-geometric mean, we can summarize the paper by writing:





π
=



4


a
g
m

(
1
,
1

/



2



)

2




1
−

∑

n
=
1


∞



2

n
+
1


(

a

n


2


−

g

n


2


)





{\displaystyle \pi ={\frac {4\;\mathrm {agm} (1,1/{\sqrt {2}})^{2}}{1-\sum \limits _{n=1}^{\infty }2^{n+1}(a_{n}^{2}-g_{n}^{2})}}}



This allows you to make the approximation, for any large   N:





π
≈



4


a

N


2




1
−

∑

k
=
1


N



2

k
+
1


(

a

k


2


−

g

k


2


)





{\displaystyle \pi \approx {\frac {4\;a_{N}^{2}}{1-\sum \limits _{k=1}^{N}2^{k+1}(a_{k}^{2}-g_{k}^{2})}}}



The purpose of this task is to demonstrate how to use this approximation in order to compute a large number of decimals of 



π


{\displaystyle \pi }

.

library(Rmpfr)
 
agm <- function(n, prec) {
  s <- mpfr(0, prec)
  a <- mpfr(1, prec)
  g <- sqrt(mpfr("0.5", prec))
  p <- as.bigz(4)
  for (i in seq(n)) {
    m <- (a + g) / 2L
    g <- sqrt(a * g)
    a <- m
    s <- s + p * (a * a - g * g)
    p <- p + p
  }
  4L * a * a / (1L - s)
}
 
1e6 * log(10) / log(2)
# [1] 3321928
 
# Compute pi to one million decimal digits:
p <- 3322000
x <- agm(20, p)
 
# Check
roundMpfr(x - Const("pi", p), 64)
# 1 'mpfr' number of precision  64   bits 
# [1] 4.90382361286485830568e-1000016
 
# Save to disk
f <- file("pi.txt", "w")
writeLines(formatMpfr(x, 1e6), f)
close(f)