Here is an adaptation of MATLAB code from Trefethen, Ten Digit Algorithms (2005), based on Borwein & Borwein, The Arithmetic-Geometric Mean and Fast Computation of Elementary Functions (1984) that calculates $\pi$ via the AGM method.
ClearAll[npi];
npi[digits_] := Block[{two, iter},
iter[{x_, y_, p_}] :=
With[{s = Sqrt@x},
{(s + 1/s)/2,
(y*s + 1/s)/(1 + y),
p*(1 + x)/(1 + y)}];
two = SetPrecision[2, 1 + digits];
With[{y = Sqrt@(Sqrt@two), eps = 10.^(-digits/2)},
NestWhile[
iter,
{(y + 1/y)/2, y, two + Sqrt@2},
Abs[Last@#1 - Last@#2] > eps &,
2]
]
];
Examples:
ClearSystemCache["Numeric"]ClearSystemCache[] (* clears cached values of Pi *)
digits = 10^6;
N[Pi, digits]; // AbsoluteTiming
pi = Last@npi[digits]; // AbsoluteTiming
pi - Pi // Abs
(*
{0.393234, Null}
{6.24918, Null}
0.*10^-1000000
*)
ClearSystemCache["Numeric"] (* clears cached values of Pi *)
digits = 7000;
N[Pi, digits]; // AbsoluteTiming
pi = Last@npi[digits]; // AbsoluteTiming
pi - Pi // Abs
(*
{0.001243, Null}
{0.008355, Null}
0.*10^-7000
*)
The AGM algorithm is asymptotically quadratically convergent. Below are the number of digits of accuracy and the ratio with the previous step. It takes 12 iterations to reach 7000 digits and 19 to reach a million digits (in fact, about 1.43 million); in general it will take around Log2[digits] - 1 iterations.
