はじめる前に

初投稿です
n番煎じです
新規性はありません
ほぼperiaさんの記事を~~パクり~~参考にしました。ありがとうございます

はじめに

3月14日は円周率の日でした。世間では『円周率は割り切れない』ことにあやかって、入籍や挙式をされる方もいらっしゃるそうです。

そこで、円周率を1億桁求めてお祝いすることにしました。

手法

ライブラリなど

普通のプログラミング言語の変数では、小数点以下1億桁も保持してくれません。そこで任意桁の数を扱うために、GNU MPを使います。バージョンは6.3.0です。

またプログラミング言語にはC言語¹ を、コンパイラにはGCCを採用します。

実行環境

Azureの仮想マシンを使います。サイズは「Standard D4s v3」で、vCPU数は4，RAMは16GiBです。²
またCPUのモデルは「Intel(R) Xeon(R) CPU E5-2673 v4」だそうです。

アルゴリズム

Chudnovskyの公式

円周率の計算には、現状では最速と言われるChudnovskyの公式を採用します。

\frac{1}{\pi} = \frac{12}{\sqrt{C^3}}\sum_{n=0}^{\infty}{\frac{(-1)^n (6n)! (A+Bn)}{(3n)! (n!)^3 C^{3n}}}

\begin{align}
  A &= 13591409  \\
  B &= 545140134 \\
  C &= 640320
\end{align}

1項あたり約14桁も精度が上がるという驚異の収束速度を誇り、Googleが円周率を100兆桁計算したときにも用いられました。

Binary Splitting Algorithm

Chudnovskyの公式の計算効率を上げるために、Binary Splitting Algorithm（通称BS法）を使います。
BS法は、periaさんの記事によると

3つの整数列$\{x_i\}, \{y_i\}, \{z_i\}$が定義されているときに
\frac{1}{x_0} \left(y_0 + \frac{z_0}{x_1} \left(y_1 + \frac{z_1}{x_2} \left(\cdots + \frac{z_{n-2}}{x_{n-1}}\left(y_{n-1}\right)\right)\right)\right)
という数式の計算コストを$O(n^2)$から$O(n(\log{n^3}))$まで下げるアルゴリズム

だそうです。

掛ける数と掛けられる数の桁数を揃えることで、効率的な掛け算アルゴリズムが使えるようになるらしいです。

Binary Splittingの式変形

円周率.jpさんの記事とperiaさんの記事では式の形がかなり異なりますが、単に式変形しただけです。

\begin{split}
  &  \sum_{n=0}^{N-1}{\frac{a(n)}{b(n)}\frac{p(0) \cdots p(n)}{q(0) \cdots q(n)}} \\
  &= \frac{a(0)}{b(0)}\frac{p(0)}{q(0)} + \frac{a(1)}{b(1)}\frac{p(0) p(1)}{q(0) q(1)} + \frac{a(2)}{b(2)}\frac{p(0) p(1) p(2)}{q(0) q(1) q(2)} + \cdots + \frac{a(N-1)}{b(N-1)}\frac{p(0) \cdots p(N-1)}{q(0) \cdots q(N-1)} \\
  &= \frac{1}{q(0)}\frac{a(0) p(0)}{b(0)} + \frac{p(0)}{q(0) q(1)}\frac{a(1)p(1)}{b(1)} + \frac{p(0) p(1)}{q(0) q(1) q(2)}\frac{a(2) p(2)}{b(2)} + \cdots + \frac{p(0) \cdots p(N-2)}{q(0) \cdots q(N-1)}\frac{a(N-1) p(N-1)}{b(N-1)} \\
  &= \frac{1}{q(0)} \left(\frac{a(0) p(0)}{b(0)} + \frac{p(0)}{q(1)} \left(\frac{a(1) p(1)}{b(1)} + \frac{p(1)}{q(2)} \left(\frac{a(2) p(2)}{b(2)} + \cdots + \frac{p(N-2)}{q(N-1)} \left( \frac{a(N-1) p(N-1)}{b(N-1)} \right)\right)\right) 
\cdots \right)
\end{split}

ここで

\begin{align}
  x_n &= q(n) \\
  y_n &= \frac{a(n)p(n)}{b(n)} \\
  z_n &= p(n) \\
\end{align}

とおくと

\frac{1}{x_0} \left(y_0 + \frac{z_0}{x_1} \left(y_1 + \frac{z_1}{x_2} \left(y_2 + \cdots + \frac{z_{n-2}}{x_{n-1}}(y_{n-1})\right)\right) \cdots \right)

となり、一致します。

どうしてBinary Splitting？

紹介されている式だけではBinary Splitting2つに分割する感がないですが、こちらのサイトに、階乗に対してBinary Splittingを適用する例がありました。

普通に階乗を計算する場合、$n! = (n - 1)! \cdot n$は$大きい数 \times 小さい数$になりがちですが、隣接する項どうしを掛けていくことで、桁数を揃えることができるようです。

では紹介されている式と連続する積にどんな関係があるのかというと、行列$M_k$を

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

と定義したとき、

M_0 M_1 \cdots M_{n-1}

の1行1列目の成分が分子、1行2列目の成分が分母に対応するようです。
BS法は、このような行列積を効率よく計算する手法と言えますね。

公式の変形

Chudnovskyの公式を式変形して、BS法が適用できるようにします。periaさんの記事に載っている方法で変形すると

\frac{1}{\pi} \sim \frac{1}{426880 \sqrt{10005}} \left(A + \frac{5Y(1, n)}{X(1, n)}\right)

\begin{cases}
  X(l, u) = X(l, m) X(m, u) \\
  Y(l, u) = Y(l, m) X(m, u) + Z(l, m) Y(m, u) \\
  Z(l, u) = Z(l, m) Z(m, u)
\end{cases}

\begin{align}
  X(k, k+1) &= -k^3 C^3 / 24 \\
  Y(k, k+1) &= A + kB \\
  Z(k, k+1) &=
  \begin{cases}
    (2k + 1)(6k + 1)(6k + 5) & \text{ if } k \neq n - 1 \\
    0 & \text{ if } k = n - 1
  \end{cases}
\end{align}

\begin{cases}
  A = 13591409 \\
  B = 545140134 \\
  C = 640320
\end{cases}

となりました。
$X(k, k+1)$と$Z(k, k+1)$の正負が本家とは異なりますが、これは$-C^3 / 24$を定数とし、$Z(k, k+1)$の計算において$-1$を掛けるステップを短縮する狙いがあります。

なおitchynyさんの記事では異なる手順で式変形をされていますが、両方とも本質的には同じです。
あえて言うなら、periaさんは公式に則って変形し、itchynyさんは「2つに分割する」という面から変形しているようです。

結果の確認

結果の正しさを確認する方法はいくつかありますが、今回はGNU MPFRのmpfr_const_pi()を使って、小数点以下1億10桁まで計算させます。

$ ./a.out > result_mpfr

比較用に使うので、結果をファイルに出力させておきました。

ソースコード

MPFRのバージョンは4.2.1です。

#include <stdio.h>
#include <gmp.h>
#include <mpfr.h>

int main(void) {
    mpfr_t pi;

    mpfr_init2(pi, 332192844);
    mpfr_const_pi(pi, MPFR_RNDN);

    mpfr_printf("%.100000010Rf\n", pi);

    mpfr_clear(pi);
    mpfr_free_cache();

    return 0;
}

コーディング

導出した公式を素直に実装していきましょう。
工夫としては、構造体を用いて記述を簡潔にしたことや、一時変数を可能な限り使いまわしたことです。

ソースコード

#include <stdio.h>
#include <stdlib.h>
#include <gmp.h>

// Struct definition
typedef struct _XYZ_t {
    mpz_t X;
    mpz_t Y;
    mpz_t Z;
} XYZ_t;

// Function declaration
int calc_pi_f(mpf_t, unsigned long);
int bsa(mpz_t, mpz_t, unsigned long);
void bsa_core(XYZ_t *, unsigned long, unsigned long, mpz_t, mpz_t, mpz_t, unsigned long);

// Inline function definition
inline void XYZ_init(XYZ_t* xyz) {
    if (xyz) {
        mpz_inits(xyz->X, xyz->Y, xyz->Z, NULL);
    }
}

inline void XYZ_clear(XYZ_t* xyz) {
    if (xyz) {
        mpz_clears(xyz->X, xyz->Y, xyz->Z, NULL);
    }
}

// Main
int main(void) {
    mpf_set_default_prec(332192844);

    unsigned long n_bsa = 7051368;
    mpf_t pi;
    mpf_init(pi);

    calc_pi_f(pi, n_bsa);

    gmp_printf("%.100000010Ff\n", pi);

    mpf_clear(pi);

    return 0;
}

int calc_pi_f(mpf_t pi, unsigned long n_bsa) {
    mpq_t pi_q, tmp;
    mpf_t sqrt10005;

    // Init
    mpq_inits(pi_q, tmp, NULL);
    mpf_init(sqrt10005);

    // Do BSA
    mpz_t X_z, Y_z;
    mpz_inits(X_z, Y_z, NULL);

    int ret = bsa(X_z, Y_z, n_bsa);    // X(1,n) , Y(1,n)

    if (ret != 0) {
        mpz_clears(X_z, Y_z, NULL);

        return ret;
    }

    mpz_mul_ui(Y_z, Y_z, 5UL);     // 5Y(1,n)

    mpq_set_z(pi_q, Y_z);
    mpq_set_z(tmp, X_z);

    mpz_clears(X_z, Y_z, NULL);

    // Calc Pi
    mpq_div(pi_q, pi_q, tmp);      // 5Y(1,n) / X(1,n)

    mpq_set_ui(tmp, 13591409UL, 1UL);
    mpq_add(pi_q, tmp, pi_q);      // A + (5Y(1,n) / X(1,n))

    mpq_inv(pi_q, pi_q);           // 1 / (A + (5Y(1,n) / X(1,n)))

    mpq_set_ui(tmp, 426880UL, 1UL);
    mpq_mul(pi_q, tmp, pi_q);      // 426880 / (A + (5Y(1,n) / X(1,n)))

    mpf_set_q(pi, pi_q);
    mpf_sqrt_ui(sqrt10005, 10005UL);
    mpf_mul(pi, sqrt10005, pi);    // 426880 * sqrt(10005) / (A + (5Y(1,n) / X(1,n)))

    // Free
    mpf_clear(sqrt10005);
    mpq_clears(pi_q, tmp, NULL);

    return 0;
}

int bsa(mpz_t X, mpz_t Y, unsigned long n) {
    if (n < 3) {
        return 2;
    }

    XYZ_t xyz;
    XYZ_init(&xyz);

    mpz_t A, B, mC3over24;
    mpz_init_set_si(A, 13591409L);
    mpz_init_set_si(B, 545140134L);
    mpz_init_set_si(mC3over24, -10939058860032000L);

    bsa_core(&xyz, 1, n, A, B, mC3over24, n);

    mpz_swap(X, xyz.X);
    mpz_swap(Y, xyz.Y);

    mpz_clears(A, B, mC3over24, NULL);
    XYZ_clear(&xyz);

    return 0;
}

void bsa_core(XYZ_t *xyz, unsigned long l, unsigned long u, mpz_t A, mpz_t B, mpz_t mC3over24, unsigned long n) {
    if (l + 1 == u) {
        // X(k, k + 1) = mC3over24 k^3
        mpz_ui_pow_ui(xyz->X, l, 3UL);
        mpz_mul(xyz->X, xyz->X, mC3over24);

        // Y(k, k + 1) = A + kB
        mpz_set(xyz->Y, A);
        mpz_addmul_ui(xyz->Y, B, l);

        if (l != n - 1) {
            // Z(k, k + 1) = (2k + 1)(6k + 1)(6k + 5) if (k != n - 1)
            mpz_set_ui(xyz->Z, 2 * l + 1);
            mpz_mul_ui(xyz->Z, xyz->Z, 6 * l + 1);
            mpz_mul_ui(xyz->Z, xyz->Z, 6 * l + 5);
        } else {
            // Z(k, k + 1) = 0 if (k == n - 1)
            mpz_set_ui(xyz->Z, 0UL);
        }
    } else {
        XYZ_t xyz_next;
        XYZ_init(&xyz_next);

        int m = (l + u) / 2;

        bsa_core(xyz, l, m, A, B, mC3over24, n);
        bsa_core(&xyz_next, m, u, A, B, mC3over24, n);

        // X(l, u) = X(l, m) X(m, u)
        mpz_mul(xyz->X, xyz->X, xyz_next.X);

        // Y(l, u) = Y(l, m) X(m, u) + Z(l, m) Y(m, u)
        mpz_mul(xyz->Y, xyz->Y, xyz_next.X);
        mpz_addmul(xyz->Y, xyz->Z, xyz_next.Y);

        // Z(l, u) = Z(l, m) Z(m, u)
        mpz_mul(xyz->Z, xyz->Z, xyz_next.Z);

        XYZ_clear(&xyz_next);
    }
}

実行と計測

実行結果を、さっき作ったファイルの内容と比較してみましょう。

$ ./a.out | cmp result_mpfr
$

分かりにくいですが、cmpは比較元と比較先に差がない場合は何も出力しません。
従って、円周率1億10桁が正確に求められているようです。

次に計測です。

計測にはmultitimeコマンドを使い、10回の実行時間の平均と標準偏差を測定しました。
またgmp_printf()は、10進数に変換して出力する部分で時間を消費してしまうため、コメントアウトしました。

$ multitime -n 10 ./a.out
===> multitime results
1: ./a.out
            Mean        Std.Dev.    Min         Median      Max
real        417.435     1.472       415.373     417.588     419.323
user        413.849     1.420       411.848     413.985     415.706
sys         3.537       0.078       3.380       3.560       3.652

実行時間の平均は413.849秒、標準偏差は1.420秒でした。
だいたい7分くらいかかってますね。

まとめ・今後の展望

当初の目的通り、円周率1億桁を求めることができました。ですが、あまりにも時間がかかりすぎています。
というのも、GNU MPの開発元が提供している「gmp-chudnovsky.c」を利用すると、なんと2分ちょっとで円周率1億10桁が求まってしまうからです。

これより速くするのは難しいかもしれませんが、並列化やアルゴリズムの見直しによって、現状の実行時間を短くすることができるはずです。

参考資料

円周率.jp 「Ramanujan系公式」 http://xn--w6q13e505b.jp/formula/ramanujan.html
円周率.jp 「binary Splitting」 http://xn--w6q13e505b.jp/method/binary.html
peria 「Chudnovskyの公式を用いた円周率の計算用メモ」 https://qiita.com/peria/items/c02ef9fc18fb0362fb89
itchyny 「円周率を1億桁計算しました！ - その試行錯誤の詳しい経緯と結果 -」 https://itchyny.hatenablog.com/entry/20120304/1330870932
peria 「円周率計算の高速化 (3) 無駄な計算の除去」 https://qiita.com/peria/items/f959b3a1a3f5b8d939ef

C++など他の言語から利用することもできますが、慣れているのでC言語にしました ↩
HPC用の仮想マシンもあるそうですが、よく分かりません ↩