Last active
July 16, 2020 06:15
-
-
Save shouth/e03a95a4c7aef9bc2bdd2937897f3203 to your computer and use it in GitHub Desktop.
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
#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; | |
} |
僕のGCC(7.4.0)だと出ないですね…
GCCのバージョン教えてもらってもいいですか?
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
-std=c++14
か-std=c++1z
のオプション付けないとコンパイルエラーにならない?