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

円周率計算の高速化シリーズの2本目…という位置付けですが、まだ高速化はしません。今回は別公式の利用を可能にして計算結果の確認ができるようにしましょう。ちなみに自分で他の公式やらで計算する事に興味がなければ Google の Pi delivery サービスなどを使って「正解」が手元に用意できてるなら今回のは必要ありません。

公式

今回確認用に使う公式は Ramaujan の公式です。

\frac{1}{\pi} = \frac{2\sqrt{2}}{9801}\sum_{k=0}^{\infty}\frac{(4k)!(1103+26390k)}{(k!)^4 396^{4k}}

何となくメイン公式の Chudnovsky の公式と似てますが、同じ理論背景から導かれる公式なのでまぁそういうものかもしれません。モジュラーとか保型形式とかよくわかりません。これを Chudnovsky の公式と同様に変形すると $A=1103$、$B=26390$、$C=396$ と置換して

\begin{eqnarray}
\dfrac{1}{\pi} &=& \dfrac{2\sqrt{2}}{9801}\sum_{k=0}^{\infty}\frac{(4k)!(A+Bk)}{(k!)^4 C^{4k}}\\
&=& \dfrac{2\sqrt{2}}{9801}\left(A+\dfrac{4!(A+B)}{(1!)^4C^4}+\dfrac{8!(A+2B)}{(2!)^4C^8}+\cdots\right)\\
&=& \dfrac{2\sqrt{2}}{9801}\left(A+\dfrac{1\cdot2\cdot3\cdot4}{1^4C^4}\left(A+B+\dfrac{5\cdot6\cdot7\cdot8}{2^4C^4}\Big(A+2B+\cdots\right)\right)\\
\end{eqnarray}

となるので

\begin{eqnarray}
x_k &=& k^4C^4 \quad (ただし x_0=1)\\
y_k &=& A+Bk\\
z_k &=& (4k+1)(4k+2)(4k+3)(4k+4)=8(4k+1)(2k+1)(4k+3)(k+1)
\end{eqnarray}

で計算できます。ついでに自明な部分の約分ができるので約分しておくと

\begin{eqnarray}
x_k &=& k^3C^4/8 \quad (ただし x_0=1)\\
y_k &=& A+Bk\\
z_k &=& (4k+1)(2k+1)(4k+3)
\end{eqnarray}

となります。必要計算項数は桁数を $\log_{10} \frac{C^4}{8\cdot4\cdot2\cdot4}=7.9825407783902$ で割れば求まります。ということでコーディングしましょう。[前回]用意したフレームワークを利用すると計算変更が簡単です。

class Ramanujan : public Computer {
public:
  Ramanujan(int64_t digits) : Computer(digits) {}
  ~Ramanujan() override = default;

private:
  void setXYZ(int64_t k, mpz_class& x, mpz_class& y, mpz_class& z) override;
  void postProcess(mpz_class& x, mpz_class& y) override;
  int64_t terms(int64_t digits) const override { return digits / 7.98; }
};

void Ramanujan::setXYZ(int64_t k, mpz_class& x, mpz_class& y, mpz_class& z) {
  static constexpr int64_t A = 1103;
  static constexpr int64_t B = 26390;
  static constexpr int64_t C = 396;

  if (k == 0) {
    x = 1;
  } else {
    x = k * C / 2;
    x *= k * C / 2;
    x *= k * C / 2;
    x *= C;
  }
  y = A + B * k;
  z = 4 * k + 1;
  z *= 2 * k + 1;
  z *= 4 * k + 3;
}

void Ramanujan::postProcess(mpz_class& x, mpz_class& y) {
  x *= 9801;
  y *= 4;
  mpf_sqrt_ui(pi_.get_mpf_t(), 2);
  pi_ *= x;
  pi_ /= y;
}

あとは呼び出し元を少し変更して、計算式を簡単に変更できるようにする…つもりでしたが面倒なのでやめました。コメントアウト方式で変更です。

  std::unique_ptr<Computer> computer(new Chudnovsky(digits));
  // std::unique_ptr<Computer> computer(new Ramanujan(digits));

ついでに第2引数をもらった時は誤差計算してするようにしましょう。

void Computer::check(const char* answer_filename) {
  if (!answer_filename)
    return;

  const int64_t precision = digits_ * std::log2(10) + 10;
  std::ifstream ifs(answer_filename);
  if (!ifs.is_open())
    return;
  mpf_class answer(0, precision);
  ifs >> answer;
  answer -= pi_;

  static constexpr int kBase = 10;
  static constexpr int kDigits = 10;
  mpf_out_str(stderr, kBase, kDigits, answer.get_mpf_t());
}

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

  std::unique_ptr<Computer> computer(new Chudnovsky(digits));
  // std::unique_ptr<Computer> computer(new Ramanujan(digits));
  {
    Timer all("All");
    {
      Timer timer("Compute");
      computer->compute();
    }
    computer->output();
  }
  computer->check(answer_file);

  return 0;
}

これで例えば pi_1000.txt に 1020 桁くらいの円周率をテキスト形式で保存してたとすると

$ ./pi_drm 1000 pi_1000.txt
Compute: 0.000
    All: 0.000
0.1696594446e-1006

という形で計算誤差表示してくれます。

さて、そろそろコードが長くなってきたのでファイル分けしたセットをこちらに置いておきます。これでいつでも誤差を含めて計算できます。

これで前準備が整ったので、次はようやく高速化の第1段です。

円周率計算の高速化 (2) 計算の正しさの確認

公式