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 }

.

using Printf
 
agm1step(x, y) = (x + y) / 2, sqrt(x * y)
 
function approxπstep(x, y, z, n::Integer)
    a, g = agm1step(x, y)
    k = n + 1
    s = z + 2 ^ (k + 1) * (a ^ 2 - g ^ 2)
    return a, g, s, k
end
 
approxπ(a, g, s) = 4a ^ 2 / (1 - s)
 
function testmakepi()
	setprecision(512)
	a, g, s, k = BigFloat(1.0), 1 / √BigFloat(2.0), BigFloat(0.0), 0
	oldπ = BigFloat(0.0)
	println("Approximating π using ", precision(BigFloat), "-bit floats.")
	println("   k     Error  Result")
	for i in 1:100
		a, g, s, k = approxπstep(a, g, s, k)
		estπ = approxπ(a, g, s)
		if abs(estπ - oldπ) < 2eps(estπ) break end
		oldπ = estπ
		err = abs(π - estπ)
		@printf("%4d%10.1e%68.60e\n", i, err, estπ)
	end
end
 
testmakepi()
 