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 }

.

import std.bigint;
import std.conv;
import std.math;
import std.stdio;
 
BigInt IntSqRoot(BigInt value, BigInt guess) {
    BigInt term;
    do {
        term = value / guess;
        auto temp = term - guess;
        if (temp < 0) {
            temp = -temp;
        }
        if (temp <= 1) {
            break;
        }
        guess += term;
        guess >>= 1;
        term = value / guess;
    } while (true);
    return guess;
}
 
BigInt ISR(BigInt term, BigInt guess) {
    BigInt value = term * guess;
    do {
        auto temp = term - guess;
        if (temp < 0) {
            temp = -temp;
        }
        if (temp <= 1) {
            break;
        }
        guess += term;
        guess >>= 1;
        term = value / guess;
    } while (true);
    return guess;
}
 
BigInt CalcAGM(BigInt lam, BigInt gm, ref BigInt z, BigInt ep) {
    BigInt am, zi;
    ulong n = 1;
    do {
        am = (lam + gm) >> 1;
        gm = ISR(lam, gm);
        BigInt v = am - lam;
        if ((zi = v * v * n) < ep) {
            break;
        }
        z -= zi;
        n <<= 1;
        lam = am;
    } while(true);
    return am;
}
 
BigInt BIP(int exp, ulong man = 1) {
    BigInt rv = BigInt(10) ^^ exp;
    return man == 1 ? rv : man * rv;
}
 
void main() {
    int d = 25000;
    // ignore setting d from commandline for now
    BigInt am = BIP(d);
    BigInt gm = IntSqRoot(BIP(d + d - 1, 5), BIP(d - 15, cast(ulong)(sqrt(0.5) * 1e15)));
    BigInt z = BIP(d + d - 2, 25);
    BigInt agm = CalcAGM(am, gm, z, BIP(d + 1));
    BigInt pi = agm * agm * BIP(d - 2) / z;
 
    string piStr = to!string(pi);
    writeln(piStr[0], '.', piStr[1..$]);
}