1
0

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?

More than 3 years have passed since last update.

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

Last updated at Posted at 2020-02-11

円周率計算の高速化シリーズの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段です。

1
0
0

Register as a new user and use Qiita more conveniently

  1. You get articles that match your needs
  2. You can efficiently read back useful information
  3. You can use dark theme
What you can do with signing up
1
0

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?