嗚呼、帰りたい

プログラミングのことから日常のことまで。いわゆるごった煮というものだな。

トップ > C++ > サクッと出来ちゃう！お手軽πの求め方！

サクッと出来ちゃう！お手軽πの求め方！

C++

f:id:south0829:20190812175439p:plain

この記事はあくあたん工房お盆休みアドベントカレンダー3日目の記事です．

君がもし理系学生ならば，一つ聞きたいことがある．
理系学生の諸君に聞きたい．一度は円周率の夢を見たことがあるのではないだろうか．
誰しもが一度は円周率を算出したいと思ったことがあると思う．

良い出だしが思いつきませんでした．というわけで今回は円周率をラップトップPCで算出するお話です．記事の題名は気を引くために大げさにした誇張表現が含まれています．あと，この記事は試行錯誤の過程を残したものですので，どうか厳密さを期待しないでください．

どうして円周率をラップトップ PC で求めようと思ったん？

大学の教養科目で「円周率を求めよう」という課題があって面白そうだったからうんぬんかんぬん．

円周率を導出する

本題に入る前に

どんな環境使ったん？

計算結果は正確なん？

計算結果の検証はここに上げている chudnovsky のアルゴリズム以外の(アルゴリズム|公式)は全て chudnovsky のアルゴリズムの計算結果と一応照合しました．chudnovsky のアルゴリズムの計算結果の検証についてはあまり良い方法ではありませんが，ネット上に有志が上げている計算結果と照合しました．ただ，後述する理由で概ねあっていればよしという方針で進めています．

桁数は合ってるん？

浮動小数点数の丸めの関係で数桁ほど誤差がでている箇所があります．が，ただ単に「面白そう」という理由で着手したので厳密さを求めない方針で誤差はそのままになっています．

用いる言語の選定

円周率を導出するにあたってどうせならできるだけ多くの桁数を求めたかったので，自分が扱える言語の中から速度を重視して C++ を選びました．

多倍長整数ライブラリの選定

多くの桁数を求めるとなると C++ 標準の double 型では精度が足りません．そこで多倍長(整数|浮動小数点数)ライブラリの力を借りることにします．

Google 先生に「多倍長浮動小数点数」で聞いてみると Boost のページがでてきたので最初に Boost を使っていたのですが，Boost が内部的に用いている GNU MP に C++ インターフェースがあったのでこの記事で上げるプログラムは GNU MP を用いて書きました．公式をプログラムに起こす分にはとても使いやすかったです．

Machin の公式

円周率を導出するのに最初に用いた公式は「Machin の公式」と呼ばれる以下の公式です．

\displaystyle \frac{\pi}{4} = 4 \arctan{\frac{1}{5}} - \arctan{\frac{1}{239}}

Machin の公式は以下の「Leibniz の公式」が基になっています．

\displaystyle \frac{\pi}{4} = \arctan{1}

この式は非常に整っていますが，収束が非常に遅いです．実際に動かしてみます．プログラムで用いるのは次の級数の和の形です．

\displaystyle \frac{\pi}{4} = \sum_{n=0}^{\infty} \frac{(-1)^{n}}{2n+1}

実際には級数の和を無限項まで計算することはできないので今回は n = 10000 で計算を打ち切ります．

以下がコードです．

#include<cmath>
#include<iomanip>
#include<iostream>
#include<gmpxx.h>

int main() {
    mpf_set_default_prec(100 * log2(10));
    mpf_class pi = 0;

    for (unsigned n = 0; n <= 10000; n++) {
        pi += (n % 2 == 0 ? 4_mpf : -4_mpf) / (2 * n + 1);
    }

    std::cout << std::setprecision(100) << pi << std::endl;
}

以下の出力が得られます．

3.141692643590543213460768320877940222544825752138710733999805489190209879979564374094471751246362536

n = 10000 まで計算しても小数点以下3桁までの精度しか得られません．同様にして Machin の公式で円周率を導出してみることにしましょう．Machin の公式をプログラムで動かすためにまずは以下の級数の和の形のものを利用します．

\displaystyle \frac{\pi}{4} = \sum_{n=0}^{\infty} \left\{ 4 \frac{(-1)^{n}}{2n + 1} \left(\frac{1}{5}\right)^{2n+1} - \frac{(-1)^{n}}{2n+1} \left(\frac{1}{239}\right)^{2n+1} \right\}

そしてこれをプログラムに起こしたものがこちらです．今回は n = 10 まで計算します．

#include<cmath>
#include<iomanip>
#include<iostream>
#include<gmpxx.h>

int main() {
    mpf_set_default_prec(100 * log2(10));
    mpf_class pi = 0;

    mpf_class p5 = 5, p239 = 239;
    for (unsigned n = 0; n <= 10; n++) {
        pi += (n % 2 ? -16_mpf : 16_mpf) / ((2 * n + 1) * p5)
            + (n % 2 ? 4_mpf : -4_mpf) / ((2 * n + 1) * p239);
        p5 *= 25, p239 *= 57121;
    }

    std::cout << std::fixed << std::setprecision(100 + 1) << pi << std::endl;
}

以下の出力が得られます．

3.14159265358979329474737485771534354337874722102783990394158875360366006245603313666170009334530001440

n = 10 までの計算で小数点以下16桁までの精度が得られました！単純に比較することはできませんが，Machin の公式の方が Leibniz の公式よりも収束が圧倒的に速いのは確かなようです．

さて，上記で示した Machin の公式の級数の和の形ですが，中括弧の中身について考えてみます．

\displaystyle 4 \frac{(-1)^{n}}{2n + 1} \left(\frac{1}{5}\right)^{2n+1} - \frac{(-1)^{n}}{2n+1} \left(\frac{1}{239}\right)^{2n+1}

\left(\frac{1}{5}\right)^{2n+1} と \left(\frac{1}{239}\right)^{2n+1} を見比べてみましょう．この二者のうち， n が大きくなるときより速く減少するのは 5 \lt 239 より後者の \left(\frac{1}{239}\right)^{2n+1} であることがわかります．これより，上記に示した級数の和の収束する速さは \left(\frac{1}{5}\right)^{2n+1} に影響を受けることがわかります．そこで二者を別々に扱って，級数の和が収束する速さを大体揃えてやります．具体的には， n が 1 大きくなると \left(\frac{1}{239}\right)^{2n+1} は概ね 10^{-5} 倍されるのでそれに対応する分 \sum_{n=0}^{\infty} \left(\frac{1}{5}\right)^{2n+1} を計算してやります．これらのことを考慮して得られる式が以下になります．

\displaystyle \frac{\pi}{4} \simeq 4 \sum_{n=0}^{3m+2} \frac{(-1)^{n}}{2n+1} \left(\frac{1}{5}\right)^{2n+1} - \sum_{n=0}^{m} \frac{(-1)^{n}}{2n+1} \left(\frac{1}{239}\right)^{2n+1}

この式は m が 1 大きくなるたびに少なくとも4桁ずつ精度が上がっていきます．

そしてこちらがプログラム．

#include<cmath>
#include<chrono>
#include<fstream>
#include<iomanip>
#include<iostream>
#include<gmpxx.h>

int main() {
    auto s = std::chrono::system_clock::now();

    unsigned digit = 1000000;
    mpf_set_default_prec((digit + 1) * log2(10));
    mpf_class pi = 0;
    unsigned m = digit / 4;
    
    mpz_class p5 = 5;
    for (unsigned n = 0; n < 3 * m + 2; n++) {
        pi += (n % 2 == 0 ? 16_mpf : -16_mpf) / ((2 * n + 1) * p5);
        p5 *= 25;
    }

    mpz_class p239 = 239;
    for (unsigned n = 0; n < m; n++) {
        pi += (n % 2 == 1 ? 4_mpf : -4_mpf) / ((2 * n + 1) * p239);
        p239 *= 57121;
    }

    auto e = std::chrono::system_clock::now();
    std::ofstream ofs("machin_pi.txt");
    ofs << std::fixed << std::setprecision(digit + 1) << pi << std::endl;
    std::cout << std::chrono::duration_cast<std::chrono::milliseconds>(e - s).count() << std::endl;
}

コンパイルコマンド

g++ machin_like_formula.cpp -o machin_like_formula -lgmp -lgmpxx -O3 -mtune=native -march=native -mfpmath=both

digit に大体の桁数を指定できるようにしました．「大体」と言ってるのは，前述したように浮動小数点数の丸めの都合で数桁ほど誤差が出るからです．出力については演算結果はファイルに出力し，標準出力には実行時間をミリ秒で出力するようにしました．ということで実際に実行時間を測ってみました．

digit 実行時間(ミリ秒)
10 (<1)
100 (<1)
1,000 2
10,000 687
100,000 262932

1,000,000桁の計算を始めて4時間半経っても計算結果がでなかったのでそこで打ち止めにしました．100,000桁でもそこそこな桁数ですが，私は満足できませんでした．「比較的短時間で1,000,000桁…いや，1,000,000,000桁計算したいな？」こんな思いに駆られた私はまた別のアルゴリズムを探しました．そこで新たに出会ったのが次のアルゴリズムです．

Gauss-Legendreのアルゴリズム

次に用いた「Gauss-Legendre のアルゴリズム」は Machin の公式と同様に単純明快です．

\displaystyle a_0 = 1

\displaystyle b_0 = \frac{1}{\sqrt{2}}

\displaystyle t_0 = \frac{1}{4}

\displaystyle p_0 = 1

として4つの反復式

\displaystyle a_{n+1} = \frac{a_{n} + b_{n}}{2}

\displaystyle b_{n+1} = \sqrt{a_{n} b_{n}}

\displaystyle t_{n+1} = t_{n} - p_{n} (a_{n+1} - a_{n})^{2}

\displaystyle p_{n+1} = 2 p_{n}

を解き進めたときに

\displaystyle \pi \simeq \frac{(a + b)^{2}}{4t}

で \pi の近似値を得ることができます． このアルゴリズムの背景には完全楕円積分と呼ばれる以下に示す積分があります．

\displaystyle K(k) = F(\frac{\pi}{2}, k) = \int_{0}^{\frac{\pi}{2}} \frac{d \theta}{\sqrt{1 - k^{2} \sin^{2}{\theta}}}

\displaystyle E(k) = E(\frac{\pi}{2}, k) = \int_{0}^{\frac{\pi}{2}} \sqrt{1 - k^{2} \sin^{2}{\theta}} d \theta

正直なところよく理解できていない詳細を書くと長くなるので気になる方はこちらをご覧ください．

このアルゴリズムは n が 1 大きくなると精度が2倍になるという驚きの収束の速さを誇ります．プログラムに起こすと以下のようになります．

#include<cmath>
#include<chrono>
#include<fstream>
#include<iomanip>
#include<iostream>
#include<gmpxx.h>

int main() {
    auto s = std::chrono::system_clock::now();

    const unsigned digit = 1'000'000'000;
    mpf_set_default_prec((digit + 1) * log2(10));
    mpf_class a = 1_mpf;
    mpf_class b = 1_mpf / sqrt(2_mpf);
    mpf_class t = 1_mpf / 4_mpf;
    mpf_class p = 1_mpf;
    mpf_class c;

    for (unsigned n = 0, lim = log2(digit + 1); n < lim; n++) {
        mpf_class pa = a, pb = b, pt = t, pp = p;
        a = (pa + pb) / 2_mpf;
        b = sqrt(pa * pb);
        c = pa - a;
        t = pt - pp * c * c;
        p = 2_mpf * pp;
    }

    mpf_class n = a + b;
    mpf_class pi = n * n / (4_mpf * t);

    auto e = std::chrono::system_clock::now();
    std::ofstream ofs("gauss_legendre_pi.txt");
    ofs << std::fixed << std::setprecision(digit + 1) << pi << std::endl;
    std::cout << std::chrono::duration_cast<std::chrono::milliseconds>(e - s).count() << std::endl;
}

コンパイルコマンド

g++ gauss_legendre_algorithm.cpp -o gauss_legendre_algorithm -lgmp -lgmpxx -O3 -mtune=native -march=native -mfpmath=both

こちらも同じように実行時間を計測してみました．

digit 実行時間(ミリ秒)
10 (<1)
100 (<1)
1,000 (<1)
10,000 2
100,000 74
1,000,000 1482
10,000,000 27329
100,000,000 526097

Machin の公式と比べると圧倒的に速いです．収束の速さは伊達ではありません．しかしながら，1,000,000,000桁を求めるには1時間以上時間が掛かりそうです．「もっと速いアルゴリズムはないものか？」そこで見つけたのが次のアルゴリズムとなります．

Chudnovsky のアルゴリズム

最後にたどり着いたのは「Chudnovsky のアルゴリズム」です．このアルゴリズムは以下の等式で表されます．

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

A = 13591409

B = 545140134

C = 640320

このまま計算するとこのアルゴリズムは速さがでませんが，最適化を施すことができます．まず，

\displaystyle S(n_1, n_2) = \sum_{n = n_1 + 1}^{n_2} a_n \prod_{k = n_1+1}^{n} \frac{p_k}{q_k}

\displaystyle a_n = (-1)^{n}(A + Bn)

\displaystyle p_n = (2n - 1)(6n - 5)(6n - 1)

\displaystyle q_n = \frac{n^{3}C^{3}}{24}

とおいて，

\displaystyle P(n_1, n_2) = \prod_{n = n_1 + 1}^{n_2} p_k

\displaystyle Q(n_1, n_2) = \prod_{n = n_1 + 1}^{n_2} q_k

\displaystyle T(n_1, n_2) = P(n_1, n_2) Q(n_1, n_2)

とすると，

\displaystyle P(n_1, n_2) = P(n_1, m)P(m, n_2)

\displaystyle Q(n_1, n_2) = Q(n_1, m)Q(m, n_2)

\displaystyle T(n_1, n_2) = T(n_1, m)Q(m, n_2) + P(n_1, m)T(m, n_2)

と再帰的に評価することができます．binary splitting algorithm というようです．これらから，

\displaystyle \pi \simeq \frac{Q(0, N)}{12T(0, N) + 12AQ(0, N)} C^{3/2}

を得ることができます．

実のところを言うと，このアルゴリズムの節については先駆者様の記事とほぼ同内容になります．ですから，詳細についてはそちらを見ていただくとして，この記事及び記事中で挙げておられる論文を基に少し僕の要素を加えた僕のプログラムを紹介しようと思います．

#include<chrono>
#include<cmath>
#include<iomanip>
#include<iostream>
#include<fstream>
#include<string>
#include<utility>
#include<gmpxx.h>
#include<omp.h>
#include<unistd.h>

const mpz_class A = 13591409;
const mpz_class B = 545140134;
const mpz_class C = 640320;
const mpz_class CP3 = C * C * C;
const mpz_class CP3D24 = CP3 / 24;

void compute_pqt(const mpz_class &n1, const mpz_class &n2, mpz_class &p, mpz_class &q, mpz_class &t) {
    if (n1 + 1 == n2) {
        mpz_class n2p2 = n2 * n2;
        mpz_class n2p3 = n2p2 * n2;
        p = -72 * n2p3 + 108 * n2p2 - 46 * n2 + 5;
        q = n2p3 * CP3D24;
        t = (A + B * n2) * p;
    } else {
        mpz_class m = (n1 + n2) >> 1;
        mpz_class p1, q1, t1;
        mpz_class p2, q2, t2;

        #pragma omp task shared(p1, q1, t1)
        compute_pqt(n1, m, p1, q1, t1);
        #pragma omp task shared(p2, q2, t2)
        compute_pqt(m, n2, p2, q2, t2);
        #pragma omp taskwait
        
        p = p1 * p2;
        q = q1 * q2;
        t = t1 * q2 + p1 * t2;
    }
}

int main(int argc, char* argv[]) {
    bool show = false, fout = false;
    std::string file = "pi.txt";
    unsigned digit = 100'000'000;
    opterr = 0;

    for (char opt; (opt = getopt(argc, argv, "d:sf:F")) != -1;) {
        switch (opt) {
        case 'd':
            if (optarg[0] == '-') goto show_usage;
            digit = atoi(optarg);
            break;
        case 's':
            show = true;
            break;
        case 'f':
            if (optarg[0] == '-') goto show_usage;
            fout = true;
            file = optarg;
            break;
        case 'F':
            fout = true;
            break;
        default:
        show_usage:
            std::cout << "Usage: " << argv[0] << " [-d argument] [-f argment] [-F] [-s]" << std::endl;
            return -1;
        }
    }
    
    const unsigned prec = (digit + 1) * log2(10);
    
    #ifdef _OPENMP
    std::cout << "[INFO] OpenMP : Enabled (Max # of threads = " << omp_get_max_threads() << ")" << std::endl;
    #else
    std::cout << "[INFO] OpenMP : Disabled" << std::endl;
    #endif
    std::cout << "[INFO] Calculation begin" << std::endl;
    
    auto s = std::chrono::system_clock::now();
    mpf_set_default_prec(prec);
    mpz_class p, q, t;
    mpf_class a, b;
    #pragma omp parallel shared(p, q, t)
    {
        #pragma omp single
        compute_pqt(0, ceil(digit / 14.0), p, q, t);

        #pragma omp sections
        {
            #pragma omp section
            a = mpf_class(q) / (12 * (t + A * q));
            #pragma omp section
            b = sqrt(mpf_class(CP3));
        }
    }
    mpf_class pi = a * b;
    auto e = std::chrono::system_clock::now();
    
    if (show) {
        std::cout << std::fixed << std::setprecision(digit) << pi << std::endl;
    }

    std::cout << "[INFO] Calculation time : " << std::chrono::duration_cast<std::chrono::milliseconds>(e - s).count() << " ms" << std::endl;    

    if (fout) {
        std::cout << "[INFO] Printing the calculation result on " << file << "" << std::endl;
        std::ofstream ofs(file);
        ofs << std::fixed << std::setprecision(digit) << pi << std::endl;
        std::cout << "[INFO] Printed" << std::endl;
    }

    return 0;
}

コンパイルコマンド

g++ chudnovsky_algorithm.cpp -o chudnovsky_algorithm -lgmp -lgmpxx -fopenmp -O3 -mtune=native -march=native -mfpmath=both

今回はコマンドっぽく扱えるようにしました．アルゴリズム部分は先駆者様のものと概ね同様ですが，僕はこれに並列化を施しました．コード中にある #pragma omp ... というのが該当します．論文の方には並列化についても述べられていましたが，そこまで理解して実装するのが大変手間アルゴリズムの方に興味をもってやっていながらもできれば桁数のために速さが欲しい程度だったので軽い並列化としてOpenMPを用いる程度にとどまりました．OpenMPを無効化した状態と有効化した状態でそれぞれ実行時間を計測してみました．

digit OpenMP無効化下での実行時間(ミリ秒) OpenMP有効化下での実行時間(ミリ秒)
10 (<1) 14
100 (<1) 15
1,000 (<1) 15
10,000 1 17
100,000 26 60
1,000,000 456 590
10,000,000 8202 8261
100,000,000 174354 108325
1,000,000,000 2836657 1918724

1,000,000桁までは並列化のコストの方が並列化して得られる速さを上回っていますが，10,000,000桁になると並列化して得られる速さの方が支配的になってきます．これで目標だった円周率1,000,000,000桁を算出することができました！！

円周率算出を得て

やりごたえがありました．最初に適当にプログラムを書いて改良を加えていくたびに少しずつ速くなっていくのはやはり楽しい過程でした．例えばですが，最初の Machin の公式の部分はまさにその試行錯誤の過程がそのまま反映されています．また，累乗については最初はループごとに最初から計算し直す形式を取っていたのですが，ループごとに前のループまでの計算結果に対して演算を施すように変更したところ，こちらも劇的に速度が向上しました．最後の Chudnovsky のアルゴリズムで書いた compute_pqt() は最初はp，q，tをタプルに入れて返すようにしてました．最近覚えたばかりの構造化束縛を使ってみたかったのです．しかしながら速さはでませんでした．そこで参照渡しの形に書き換えてやると先程とは打って変わって非常に速くなったのです．また，アルゴリズムを変えると速度が向上するのも非常に興味深いことでした．

空間計算量について考える機会にもなりました．1,000,000,000桁を算出するべくプログラムを実行していた際，8GBのRAMでは足りずに

GNU MP: Cannot allocate memory (size=537133072)
fish: './gauss_legendre_algorithm' terminated by signal SIGABRT (Abort)

というメッセージと共にプログラムが落とされてしまったので，脳筋的な手法ではありますが手動で16GB分のスワップファイルを用意して解決するなんてことがありました．普段，時間計算量は意識していたものの空間計算量は意識したことがなかったのでややビビりました．

更に並列化の効果とトレードオフについて学ぶ機会にもなりました．「計算量が小さいと並列化のコストが大きすぎる」ということを知ってはいたのですが，いざ実際に計測してみるとそのとおりで軽く感動を覚えました．また計算量が大きいと並列化が効果的になってくるのも非常に面白いと思いました．

この記事はこれで以上となります．最後まで読んでいただきありがとうございました．