Help us understand the problem. What is going on with this article?

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

円周率計算の高速化のその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}

で計算できるはずです。必要計算項数は桁数を $\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 *= k * C;
  }
  y = A + B * k;
  z = 4 * k + 1;
  z *= 2 * k + 1;
  z *= 4 * k + 3;
  z *= k + 1;
}

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

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

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

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
Comments
No comments
Sign up for free and join this conversation.
If you already have a Qiita account
Why do not you register as a user and use Qiita more conveniently?
You need to log in to use this function. Qiita can be used more conveniently after logging in.
You seem to be reading articles frequently this month. Qiita can be used more conveniently after logging in.
  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
ユーザーは見つかりませんでした