chudnovsky_algorithm.cpp

#include<chrono>
#include<cmath>
#include<iomanip>
#include<iostream>
#include<fstream>
#include<string>
#include<utility>
#include<gmpxx.h>
#include<omp.h>
#include<unistd.h>

const mpz_class A = 13591409;
const mpz_class B = 545140134;
const mpz_class C = 640320;
const mpz_class CP3 = C * C * C;
const mpz_class CP3D24 = CP3 / 24;

void compute_pqt(mpz_class n1, mpz_class n2, mpz_class &p, mpz_class &q, mpz_class &t) {
    if (n1 + 1 == n2) {
        mpz_class n2p2 = n2 * n2;
        mpz_class n2p3 = n2p2 * n2;
        p = -72 * n2p3 + 108 * n2p2 - 46 * n2 + 5;
        q = n2p3 * CP3D24;
        t = (A + B * n2) * p;
    } else {
        mpz_class m = (n1 + n2) >> 1;
        mpz_class p1, q1, t1;
        mpz_class p2, q2, t2;

        #pragma omp task shared(p1, q1, t1)
        compute_pqt(n1, m, p1, q1, t1);
        #pragma omp task shared(p2, q2, t2)
        compute_pqt(m, n2, p2, q2, t2);
        #pragma omp taskwait
        
        p = p1 * p2;
        q = q1 * q2;
        t = t1 * q2 + p1 * t2;
    }
}

int main(int argc, char* argv[]) {
    bool show = false, fout = false;
    std::string file = "pi.txt";
    unsigned digit = 100'000'000;
    unsigned prec;
    opterr = 0;

    for (char opt; (opt = getopt(argc, argv, "d:sf:F")) != -1;) {
        switch (opt) {
        case 'd':
            if (optarg[0] == '-') goto show_usage;
            digit = atoi(optarg);
            break;
        case 's':
            show = true;
            break;
        case 'f':
            if (optarg[0] == '-') goto show_usage;
            fout = true;
            file = optarg;
            break;
        case 'F':
            fout = true;
            break;
        default:
        show_usage:
            std::cout << "Usage: " << argv[0] << " [-d argument] [-f argment] [-F] [-s]" << std::endl;
            return -1;
        }
    }
    
    prec = (digit + 1) * log2(10);
    
    #ifdef _OPENMP
    std::cout << "[INFO] OpenMP : Enabled (Max # of threads = " << omp_get_max_threads() << ")" << std::endl;
    #else
    std::cout << "[INFO] OpenMP : Disabled" << std::endl;
    #endif
    std::cout << "[INFO] Calculation begin" << std::endl;
    
    mpf_set_default_prec(prec);
    mpz_class p, q, t;
    mpf_class a, b;
    auto s = std::chrono::system_clock::now();
    #pragma omp parallel shared(p, q, t)
    {
        #pragma omp single
        compute_pqt(0, ceil(digit / 14.0), p, q, t);

        #pragma omp sections
        {
            #pragma omp section
            a = mpf_class(q) / (12 * (t + A * q));
            #pragma omp section
            b = sqrt(mpf_class(CP3));
        }
    }
    mpf_class ans = a * b;
    auto e = std::chrono::system_clock::now();
    
    if (show) {
        std::cout << std::fixed << std::setprecision(digit) << ans << std::endl;
    }

    std::cout << "[INFO] Calculation time : " << std::chrono::duration_cast<std::chrono::milliseconds>(e - s).count() << " ms" << std::endl;    

    if (fout) {
        std::cout << "[INFO] Printing the calculation result on " << file << "" << std::endl;
        std::ofstream ofs(file);
        ofs << std::fixed << std::setprecision(digit) << ans << std::endl;
        std::cout << "[INFO] Printed" << std::endl;
    }

    return 0;
}

To run this program, make sure that gcc, gmp and gmpxx are installed on your computer and compile this program with the command below

g++ chudnovsky_algorithm.cpp -o chudnovsky_algorithm -lgmp -lgmpxx -fopenmp -O3 -mtune=native -march=native -mfpmath=both

and run ./chudnovsky_algorithm. By default, this program calculates π up to 100 million decimal places and take about 110 seconds to have it done (as far as I run this program on my laptop whose CPU is Intel Core i5 8250U). Without any options, this program doesn't print π.

Options:

• -d Calculate π up to the decimal places given as argument.
• -f Print the calculation result on the file given as argument.
• -F Print the calculation result on pi.txt.
• -s Print the calculation result on stdout.

Reference:
https://bellard.org/pi/pi2700e9/pipcrecord.pdf

-std=c++14か-std=c++1zのオプション付けないとコンパイルエラーにならない？

僕のGCC(7.4.0)だと出ないですね…
GCCのバージョン教えてもらってもいいですか？