Qiita Teams that are logged in
You are not logged in to any team

Log in to Qiita Team
Community
OrganizationAdvent CalendarQiitadon (β)
Service
Qiita JobsQiita ZineQiita Blog
Help us understand the problem. What is going on with this article?

円周率計算の高速化 (3)

More than 1 year has passed since last update.

円周率計算の高速化からの3本目です。前回と違って今回はちゃんと高速化します。

TL, DR;

不要な計算をスキップしたら 7 %高速化しました。

前置き

Chudnovsky の公式を用いた円周率の計算用メモで、BS 法の計算では

M_k=\begin{pmatrix} z_k & y_k \\ 0 & x_k \end{pmatrix}

のもとで

\begin{pmatrix}Y\\X\end{pmatrix}
= M_0M_1M_2\cdots M_{n-1} \begin{pmatrix}0&1\end{pmatrix}

を計算するという説明をしました。ところで$Z$はどこへ行ったのでしょう?

ということで

そう、$Z$は最後の計算では計算する必要がないのです。ということはその前の段階で$Y$を計算するのに使う$Z$も必要なくて…と再帰的につながります。これをプログラムに導入しましょう。

void Computer::drm(const int64_t n0,
                   const int64_t n1,
                   mpz_class& x0,
                   mpz_class& y0,
                   mpz_class& z0,
                   bool need_z) {  // ここに追加
  if (n0 + 1 == n1) {
    setXYZ(n0, x0, y0, z0);
    return;
  }

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

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

  x0 *= x1;
  if (need_z)  // この if 分岐追加
    z0 *= z1;
}

void Computer::compute() { 
  const int64_t precision = digits_ * std::log2(10) + 10;
  pi_.set_prec(precision); 

  const int64_t n = terms(digits_);
  mpz_class x, y, z;
  drm(0, n, x, y, z, false);  // 最終的に z は要らない
  postProcess(x, y);
}

この3行くらい(と対応する.hファイルの関数宣言)を編集するだけです。

実験

Chudnovsky 公式で 1000 万桁を計算する時間をそれぞれ5回ずつ計測しました。基底変換と出力はGMPライブラリの中の問題なので、それらの時間は含んでいません。

1回目 2回目 3回目 4回目 5回目 平均
変更前 6.214 6.488 6.545 6.433 6.467 6.429
変更後 5.886 5.922 6.004 6.024 6.046 5.976

ということで7%くらい短時間で計算できているようです。(Git)

Why not register and get more from Qiita?
  1. We will deliver articles that match you
    By following users and tags, you can catch up information on technical fields that you are interested in as a whole
  2. you can read useful information later efficiently
    By "stocking" the articles you like, you can search right away