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 }

.

/*REXX program calculates the value of  pi  using the  AGM  algorithm.                  */
parse arg d .;   if d=='' | d==","  then d= 500  /*D  not specified?  Then use default. */
numeric digits d+5                               /*set the numeric decimal digits to D+5*/
z= 1/4;                a= 1;     g= sqrt(1/2)    /*calculate some initial values.       */
n= 1
       do j=1   until a==old;       old= a       /*keep calculating until no more noise.*/
       x= (a+g) * .5;               g= sqrt(a*g) /*calculate the next set of terms.     */
       z= z - n*(x-a)**2;  n= n+n;  a= x         /*Z  is used in the final calculation. */
       end   /*j*/                               /* [↑]  stop if  A  equals  OLD.       */
 
pi= a**2 / z                                     /*compute the finished  value of  pi.  */
numeric digits d                                 /*set the numeric decimal digits to  D.*/
say pi / 1                                       /*display the computed value of  pi.   */
exit 0                                           /*stick a fork in it,  we're all done. */
/*──────────────────────────────────────────────────────────────────────────────────────*/
sqrt: procedure; parse arg x;  if x=0  then return 0;  d=digits();  numeric digits;  h=d+6
      numeric form; m.=9; parse value format(x,2,1,,0) 'E0' with g "E" _ .; g=g *.5'e'_ %2
            do j=0  while h>9;        m.j=h;                 h=h%2+1;          end  /*j*/
            do k=j+5  to 0  by -1;    numeric digits m.k;    g=(g+x/g)*.5;     end  /*k*/
      numeric digits d;     return g/1