円周率計算の高速化 (3) 無駄な計算の除去

• 円周率計算の高速化(1) 比較元の再構成と問題設定
• 円周率計算の高速化(2) 計算の正しさの確認
• 円周率計算の高速化(3) 無駄な計算の除去。約7％高速に
• 円周率計算の高速化(4) 並列化。60％の時間削減

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

前置き

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)

次は並列化をしてもっと高速にしましょう。