Skip to content

Instantly share code, notes, and snippets.

@shouth
Last active July 16, 2020 06:15
Show Gist options
  • Save shouth/e03a95a4c7aef9bc2bdd2937897f3203 to your computer and use it in GitHub Desktop.
Save shouth/e03a95a4c7aef9bc2bdd2937897f3203 to your computer and use it in GitHub Desktop.
#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;
}
@HansRobo
Copy link

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

@shouth
Copy link
Author

shouth commented Jul 31, 2019

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

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment