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 }

.

 
program Calculate_Pi;
 
{$APPTYPE CONSOLE}
 
uses
  System.SysUtils,
  Velthuis.BigIntegers,
  System.Diagnostics;
 
function IntSqRoot(value, guess: BigInteger): BigInteger;
var
  term: BigInteger;
begin
  while True do
  begin
    term := value div guess;
    if (BigInteger.Abs(term - guess) <= 1) then
      break;
    guess := (guess + term) shr 1;
  end;
  Result := guess;
end;
 
function ISR(term, guess: BigInteger): BigInteger;
var
  value: BigInteger;
begin
  value := term * guess;
  while (True) do
  begin
    if (BigInteger.Abs(term - guess) <= 1) then
      break;
    guess := (guess + term) shr 1;
    term := value div guess;
  end;
  Result := guess;
end;
 
function CalcAGM(lam, gm: BigInteger; var z: BigInteger; ep: BigInteger): BigInteger;
var
  am, zi, v: BigInteger;
  n: UInt32;
begin
  n := 1;
  while True do
  begin
    am := (lam + gm) shr 1;
    gm := ISR(lam, gm);
    v := am - lam;
    zi := v * v * n;
    if (zi < ep) then
      break;
    z := z - zi;
    n := n shl 1;
    lam := am;
  end;
  Result := am;
end;
 
function BIP(exp: Integer; man: UInt32 = 1): BigInteger;
begin
  Result := man * BigInteger.Pow(10, exp);
end;
 
function Compress(val: string; size: Integer): string;
begin
  result := val.Remove(size, val.Length - size * 2).Insert(size, '...');
end;
 
const
  DEFAULT_DIGITS = 25000;
 
var
  d: Integer;
  am, gm, z, agm, pi: BigInteger;
  StopWatch: TStopwatch;
  s: string;
 
begin
  StopWatch := TStopwatch.Create;
 
  d := DEFAULT_DIGITS;
  if (ParamCount > 0) then
  begin
    d := StrToIntDef(ParamStr(1), d);
 
    if ((d < 1) or (d > 999999)) then
      d := DEFAULT_DIGITS;
  end;
 
  StopWatch.Start;
 
  am := BIP(d);
 
  gm := IntSqRoot(BIP(d + d - 1, 5), BIP(d - 15, Trunc(Sqrt(0.5) * 1e+15)));
 
  z := BIP(d + d - 2, 25);
  agm := CalcAGM(am, gm, z, BIP(d + 1));
 
  pi := (agm * agm * BIP(d - 2)) div z;
  s := pi.ToString.Insert(1, '.');
 
  StopWatch.Stop;
  Writeln(Format('Computation time: %.3f seconds ', [StopWatch.ElapsedMilliseconds
    / 1000]));
 
  Writeln(Compress(s, 20));
  readln;
end.
 