gauss_legendre_algorithm.cpp

#include<cmath>
#include<iomanip>
#include<iostream>
#include<gmpxx.h>

constexpr unsigned DIGIT = 10'000'001;

int main() {
    mpf_class a = 1_mpf;
    mpf_class b = 1_mpf / sqrt(2_mpf);
    mpf_class t = 1_mpf / 4_mpf;
    mpf_class p = 1_mpf;
    mpf_class d;
    mpf_class n;

    for (mpf_class n = 0_mpf, lim = log2(DIGIT); n < lim; n++) {
        mpf_class pa = a, pb = b, pt = t, pp = p;
        a = (pa + pb) / 2_mpf;
        b = sqrt(pa * pb);
        d = pa - a;
        t = pt - pp * d * d;
        p = 2_mpf * pp;
    }

    n = a + b;
    std::cout << std::setprecision(DIGIT) << n * n / (4_mpf * t) << std::endl;
}