条件設定とか

itchyny さんと同じ点としては

主計算には Chudnovsky の公式を使います。
円周率の小数点以下1億桁の計算を行い、確認を行います。
GMP 使います。

逆に新しくしようと思う点は

C++ 使います。
確認法を変更します。
高速化したいです。できればマルチスレッド化します。

というあたりです。

イントロ

Chudnovsky の公式、Binary Splitting 法、分数の簡略化、計算項数の設定、などについては3年前の記事をご覧ください。(ちなみに少し話題の焦点がずれますが、Binary Splitting の説明に出てくる $P$, $Q$, $S$, $T$ と数式の対応、計算量削減の流れとか分かりにくくないですか？私は DRM のアプローチが分かりやすいのでそちらを毎度使っています。)

GMP の導入では configure するときに --enable-cxx をつけて gmpxx.h を使えるようにしましょう。

$ ./configure --enable-cxx
$ make
$ make check
$ sudo make install

先行研究

ということでまずは itchyny さんのコードを少し変更してコマンドラインから桁数指定できるようにしたり出力を少しわかりやすくしたりしました。git diff すると

diff --git a/itchyny.c b/itchyny.c
--- a/itchyny.c
+++ b/itchyny.c
@@ -8,6 +8,8 @@
 #include <math.h>
 #include <time.h>

+#define DEFAULT_DIGITS 10000
+
 mpz_t A, B, C, C3over24;

 void computePQT(mpz_t P, mpz_t Q, mpz_t T, int n1, int n2) {
@@ -59,11 +61,15 @@ void computePQT(mpz_t P, mpz_t Q, mpz_t T, int n1, int n2) {
 }

 int main (int argc, char* argv[]) {
-
   clock_t start, end;
   FILE* fp;

-  int digits = 10000;
+  int digits = 0;
+  if (argc > 1)
+    digits = atoi(argv[1]);
+  if (digits <= 0)
+    digits = DEFAULT_DIGITS;
+
   int prec = digits * log2(10);
   int n = digits / 14;

@@ -111,13 +117,13 @@ int main (int argc, char* argv[]) {
   mpf_div_ui(pi, pi, 12);

   /* mpf_out_str(stdout, 10, digits, pi); */
-  fp = fopen("output.txt", "w");
+  fp = fopen("pi_itchyny.out", "w");
   if (fp == NULL) {
-    printf("couldn't open output.txt");
+    printf("Could not open the file.\n");
     exit(EXIT_FAILURE);
   }
   end = clock();
-  printf("%.3f s\n",(double)(end - start) / CLOCKS_PER_SEC);
+  printf("Compute: %.3f s\n",(double)(end - start) / CLOCKS_PER_SEC);
   mpf_out_str(fp, 10, digits, pi);

   mpf_clear(pi);
@@ -131,7 +137,7 @@ int main (int argc, char* argv[]) {
   mpz_clear(C3over24);

   end = clock();
-  printf("%.3f s\n",(double)(end - start) / CLOCKS_PER_SEC);
+  printf("All    : %.3f s\n",(double)(end - start) / CLOCKS_PER_SEC);

   return 0;
 }

となりました。あとはテスト走行で計算時間を見ましょう。

桁数	計算時間(s)	合計時間(s)
100万	0.354	0.436
1000万	6.069	7.596
1億	108.527	134.432

ちょっと不思議な感じにも見えますが、多分キャッシュサイズとかメモリのswapとか関係してるんだと思います。元の記事に比べると6倍くらい速くなってるのは CPU 性能向上と GMP やコンパイラの最適化のおかげでしょう。

実装

ということで実装です。こいつが今後の高速化などなどのスタート地点になります。Gitはこちら。

#include <chrono>
#include <cmath>
#include <cstdint>
#include <cstdio>
#include <cstdlib>
#include <fstream>
#include <iomanip>
#include <iostream>
#include <memory>
#include <string>

#include <gmpxx.h>

using Integer = mpz_class;
using Float = mpf_class;

constexpr int kDefaultDigits = 10000;

class Timer {
  using Clock = std::chrono::system_clock;
  using Ms = std::chrono::milliseconds;

 public:
  Timer(const std::string& name) : name_(name), start_(Clock::now()) {}
  ~Timer() {
    auto duration = Clock::now() - start_;
    std::printf("%7s: %.3f\n", name_.c_str(),
                std::chrono::duration_cast<Ms>(duration).count() * 1e-3);
  }

 private:
  std::string name_;
  Clock::time_point start_;
};

void drm(const int64_t n0,
         const int64_t n1,
         Integer& x0,
         Integer& y0,
         Integer& z0) {
  if (n0 + 1 == n1) {
    static constexpr int64_t A = 13591409;
    static constexpr int64_t B = 545140134;
    static constexpr int64_t C = 640320;

    if (n0 == 0) {
      x0 = 1;
    } else {
      x0 = n0 * C;
      x0 *= n0 * C;
      x0 *= n0 * C / 24;
    }
    y0 = A + B * n0;
    z0 = -(6 * n0 + 5);
    z0 *= 6 * n0 + 1;
    z0 *= 2 * n0 + 1;
    return;
  }

  Integer x1, y1, z1;
  int64_t m = (n0 + n1) / 2;
  drm(n0, m, x0, y0, z0);
  drm(m, n1, x1, y1, z1);

  // y0 = x1 * y0 + y1 * z0;
  y0 *= x1;
  y1 *= z0;
  y0 += y1;

  x0 *= x1;
  z0 *= z1;
}

void compute(const int64_t digits, Float& pi) {
  const int64_t precision = digits * std::log2(10);
  mpf_set_default_prec(precision);
  pi.set_prec(precision);

  const int64_t n = digits / 14;
  Integer x, y, z;
  drm(0, n, x, y, z);

  x *= 426880;
  mpf_sqrt_ui(pi.get_mpf_t(), 10005);
  pi *= x;
  pi /= y;
}

void output(Float& pi) {
  std::ofstream ofs("pi_peria.out");
  if (!ofs.is_open()) {
    std::cerr << "Could not open the file\n";
    return;
  }
  int64_t precision = pi.get_prec();
  int64_t digits = precision / std::log2(10) + 1;
  ofs << std::setprecision(digits) << pi;
}

int main(int argc, char* argv[]) {
  int64_t digits = (argc > 1) ? std::atoi(argv[1]) : kDefaultDigits;
  if (digits <= 0)
    digits = kDefaultDigits;

  Timer all("All");
  Float pi;
  {
    Timer timer("Compute");
    compute(digits, pi);
  }
  output(pi);

  return 0;
}

時間計測としては

桁数	計算時間(s)	合計時間(s)
100万	0.400	0.488
1000万	6.071	7.562
1億	116.004	144.178

となり、itchyny さんのとほぼ同じ速さですね。これを少しずつ高速化…したいですがその前に計算が正しいことを確認できるようにしましょう。

円周率計算の高速化(1) 比較元の再構成と問題設定

条件設定とか

イントロ

先行研究

実装