Compute 2/pi

chudovsky.cpp

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

static mpq_class _rational(mp_bitcnt_t iterations)
{
    mpz_class K = 6;
    mpz_class M = 1;
    mpz_class L = 13591409;
    mpz_class X = 1;
    mpz_class numerator = L;

    for (mp_bitcnt_t k = 1; k <= iterations; ++k) {
        M = M * K * (K * K - 16) / (k * k * k);
        K += 12;
        L += 545140134;
        X *= -262537412640768000;

        numerator = M * L - 262537412640768000 * numerator;
    }

    mpq_class result(numerator, X);

    result.canonicalize();

    return result;
}

static mpf_class _2opi(mp_bitcnt_t precision)
{
    mpf_class result = sqrt(mpf_class(mpq_class(1, 10005), precision + 64));

    result *= _rational(precision / 45 + 1) / 213440;

    result.set_prec(precision);

    return result;
}

int main()
{
    enum Order { Ascending = -1, Descending = 1 };
    enum Endian { Little = -1, Native = 0, Big = 1 };

    typedef std::size_t Word;

    const mp_bitcnt_t precision = 1024;
    const mp_bitcnt_t length = precision / std::numeric_limits<Word>::digits + !!(precision % std::numeric_limits<Word>::digits);

    unsigned long buffer[length];

    mpf_class object = _2opi(precision);

    mpf_mul_2exp(object.get_mpf_t(), object.get_mpf_t(), precision);

    mpz_export(buffer, nullptr, Descending, sizeof(Word), Native, 0, mpz_class(object).get_mpz_t());

    std::hex(std::cout);
    std::cout.fill('0');

    for (mp_bitcnt_t k = 0; k < length; ++k)
        std::cout << std::setw(2 * sizeof(Word)) << buffer[k] << '\n';
}