Pi

pi.cc

#include <iostream>
#include <cmath>
#include <omp.h>
#include <gmp.h>

#define LOG2_10 3.321928094887362

std::string pi(size_t places) {
    size_t scaled = ceil(static_cast<double>(places) / 14);
    mpf_set_default_prec(LOG2_10 * places);

    mpf_t fpi; mpf_init(fpi);

    omp_lock_t write_lock;
    omp_init_lock(&write_lock);

#pragma omp parallel
    {
        mpz_t ztmp, den;
        mpz_inits(ztmp, den, nullptr);
        mpf_t ftmp[2];
        mpf_t partial; mpf_init_set_ui(partial, 0);
        mpf_inits(ftmp[0], ftmp[1], nullptr);
#pragma omp for
        for (size_t k = 0; k < scaled; k++) {
            mpz_fac_ui(ztmp, 6 * k);
            mpf_set_z(*ftmp, ztmp);
            mpf_mul_ui(*ftmp, *ftmp, 545140134 * k + 13591409);
            mpz_fac_ui(den, k);

            mpz_pow_ui(den, den, 3);
            mpz_fac_ui(ztmp, 3 * k);
            mpz_mul(den, den, ztmp);
            mpz_ui_pow_ui(ztmp, 640320, 3 * k);
            if (k & 1) mpz_neg(ztmp, ztmp);
            mpz_mul(den, den, ztmp);

            mpf_set_z(ftmp[1], den);
            mpf_div(ftmp[0], ftmp[0], ftmp[1]);
            mpf_add(partial, partial, ftmp[0]);
        }
        omp_set_lock(&write_lock);
        mpf_add(fpi, fpi, partial);
        omp_unset_lock(&write_lock);
        mpf_clear(partial);
        mpf_clears(ftmp[0], ftmp[1], nullptr);
        mpz_clears(ztmp, den, nullptr);
    }
    omp_destroy_lock(&write_lock);

    mpz_t ztmp, den;
    mpz_inits(ztmp, den, nullptr);
    mpf_t ftmp[2];
    mpf_inits(ftmp[0], ftmp[1], nullptr);

    mpf_set_ui(*ftmp, 426880);
    mpf_div(fpi, *ftmp, fpi);
    mpf_sqrt_ui(*ftmp, 10005);
    mpf_mul(fpi, fpi, *ftmp);

    char buf[places + 3];
    int off = gmp_snprintf(buf, places + 2, "%.*Ff\n", places, fpi);
    buf[off] = '\0';

    mpf_clears(ftmp[0], ftmp[1], fpi, nullptr);
    mpz_clears(ztmp, den, nullptr);

    return buf;
}

int main() {
    auto pi_str = pi(50000);
    // auto pi_str = pi(10000);
    std::cout << pi_str << std::endl;
}