嗚呼、帰りたい

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

2020年になってしまったので現状確認をする

あけました.おめでとうございました.2020年になってしまった.自分が成人式に出席するような年なのがびっくりです.光陰矢のごとしとはよく言ったものだなと.

昨年は色々ありました.1年間の浪人生活を終え,第一志望だった阪大に前期で落とされてしまったものの,なんとか後期で国立の京都工芸繊維大学に拾ってもらえた(=合格した)のは非常に大きい出来事でした.私は臆病な人間ですから,後期の合格発表直前までもし滑り止めの私立に入学したときに確実に発生するであろう奨学金の負債に心を悩ませていました.

続きを読む

ScalaでSHA256を実装する

この記事はあくあたん工房アドベントカレンダー12月22日の記事です

はじめに

趣味のプログラミングの一環としてSHA256を実装したくなることは誰にでもあると思います.僕もその一人で,SHA256を実装すべく仕様書に目を通していたところ,Scalaで実装したら楽しいんじゃないかと思い,実装してみました.本当はアドベントカレンダーで投稿する本来の日までに完成させようと思っていたのですが,アドベントカレンダーの最終日である12月25日が迫るなか,どうしてもバグが取れないので未完成ですが投稿してしまうことにします.バグが取れ次第更新する予定です.バグは取れました.

実装

早速ですがコードを貼ります.

仕様書を読むとどうやら大雑把な流れは初期値6A09E667 BB67AE85 3C6EF372 A54FF53A 510E527F 9B05688C 1F83D9AB 5BE0CD19に対して入力を用いてfoldを行うものだったのでfoldLeftを用いた実装にしてみました.またその内部の二つのイテレーションfoldLeftを使って実装しました.全体としてだいぶスッキリした印象です.

おわりに

ものすごく雑な記事となってしまいましたが,大学の課題で思っていたよりも時間が取れなかったので以上にしたいと思います.最後までお読みいただきありがとうございました.

競プロで使えるかもしれないC++標準機能

f:id:south0829:20191215163918p:plain

この記事はあくあたん工房アドベントカレンダー12月15日の記事です.

adventar.org

12月14日の記事は けんけんヾ(๑╹◡╹)ノ” さんの Android JetpackのNavigationコンポーネントにChrome Custom Tabsを追加する でした.

はじめに

今年4月にAtCoderで競プロを始めてどハマりして以来,定期的にcpprefjpを読むようになった帰宅志願者です.cpprefjpを読んでいると競プロに使えそうなC++の標準機能に出会うのですが,あまり知られていないように思われるものがいくつか存在します.そこでこの記事ではそういった「知られざる便利(かもしれない)機能」を紹介したいと思います.なお僕はC++の規格自体に精通しているわけではないので,紹介内容は必ずしも正確ではないかもしれません.ご了承ください.なおこの記事では using namespace std を用いることを前提とします.

pairtuple から値を取り出す

C++14の場合

tie()を使うとこんな風に書けます.

auto p = make_pair(111, "aaa");
int x;
string y;
tie(x, y) = p;
// x = 111, s = "aaa"
auto t = make_tuple("ABC", 123, 'D');
string a;
int n;
char c;
tie(a, n, c) = t;
// a = "ABC", b = 123, c = 'D'

auto s = make_tuple(1.23, 4.56, 7.89);
double = d, e;
// 要らない値の位置には標準で用意されている
// ignoreというプレースホルダを指定すると無視できる
tie(d, ignore, e) = s;
// d = 1.23, e = 7.89

C++17の場合

構造化束縛を使いましょう.

auto p = make_pair(111, "aaa");
auto [ x, y ] = p;
// x = 111, y = "aaa"
auto t = make_tuple("ABC", 123, 'D');
auto [ a, n, c ] = t;
// a = "ABC", b = 123, c = 'D'

auto s = make_tuple(1.23, 4.56, 7.89);
// 要らない値は無視出来ない
auto [ d, unused, e ] = s;
// d = 1.23, unused = 4.56, e = 7.89

map<string, int> m;
// ここにmへの挿入処理
for (const auto& [ key, value ] : m) {
// このように範囲for文でも使える
}

変数の型を取得する

正確には式の型ですが decltype で取得出来ます.

vector<vector<pair<int, int>>> g;
// のような変数があったとして
decltype(g) h;
// とすると(実際にこんな使い方をするかは別として)
// hをvector<vector<pair<int, int>>>型の変数として宣言できる.

// 実際の使い方としては
vector<pair<int, int>> v;
auto f = [&](auto a, auto, b) -> decltype(v[0]) {
// みたいにラムダ式で明示的に返り値の型を指定する
// 必要があるときとかは使える…かもしれない
};

二乗和の平方根

std::hypot() という関数がありまして,

double x1, y1;
double x2, y2;
double dist = hypot(x1 - x2, y1 - y2);
// sqrt(pow(x1 - x2, 2) + pow(y1 - y2, 2)) と同じ結果が得られる

みたいに二点間の距離の計算が簡潔に書けるのではないかと.

同じ文字で一定の長さの文字列を作る

stringのコンストラクタには,第一引数に文字数,第二引数に文字を取って文字を文字数回繰り返した文字列を生成するものがあります.

string a(1, 'A'), b(2, 'B'), c(3, 'C');
// a = "A", b = "BB", c = "CCC"

配列の各要素に添字番号を代入する

iota()を使うと配列の各要素に添字番号を代入することが出来ます.

vector<int> v(10);
iota(begin(v), end(v), 0);
// a[0] = 0, a[1] = 1, ... , a[9] = 9

僕はUnionFind木を構造体で持つときにコンストラクタで配列を初期化するときに使ってます.

オブジェクトを構築しつつコレクションへの挿入する

emplace 系の関数を使うと簡潔に書けます.

vector<pair<int, int>> v;

// emplace 系の関数を使わない書き方
v.push_back(make_pair(12, 34));

// emplace 系の関数を使った書き方
v.emplace_back(12, 34);

vector 以外にも queuestack などのコレクションには emplace 系の関数が生えているので詳しく調べてみるといいかもしれません.

累積和を取る

partial_sum() という関数があります.

vector<int> a;
// ここに初期化処理
partial_sum(begin(a), end(a), begin(a));
// aの中身は元の配列の累積和になる.

この関数はデフォルトでは累積和の計算をしますが,bit_xor()gcd() と組み合わせることが出来ます.

vector<int> a;
vector<int> res1, res2, res3;
// ここに初期化処理
partial_sum(begin(a), end(a), begin(res2), bit_xor<>()); // 累積xor
partial_sum(begin(a), end(a), begin(res3), gcd<int, int>); // 累積gcd (C++17)

条件を満たす隣接要素を検索する

adjacent_find() を使います.

vector<int> v;
// ここに初期化処理
auto itr1 = adjacent_find(begin(v), end(v));
// 最初の等しい隣接要素の最初のイテレータが入る
auto itr2 = adjacent_find(begin(v), end(v), not_equal_to<>());
// 最初の異なる隣接要素の最初のイテレータが入る

この関数ぶっちゃけ使わない

n個の要素から最大・最小となる値を得る

n個の要素が「配列」の場合

max_element() という関数が便利です.

vector<int> v;
// ここに初期化処理
auto itr = max_element(begin(v), end(v));
// イテレータが返ってくる

n個の要素が「変数」だったり「定数」だったりする場合

max()initializer_list を突っ込みましょう.

int a = 25, b = 51, c = 19;
int m = max({ a, b, c, 38 });
// m = 51

条件によって操作する変数を変える

三項演算子は左辺値に対しても使えます.左辺値がなんのことかわからなくてもまあ例を見てみてください.

int a = 100, b = 1;
bool flag = false;
(flag ? a : b) = 25;
// a = 100, b = 25

僕は二分探索法で使ってます.

int64_t ok = 0, ng = 1e9 + 1;
while (abs(ng-ok) > 1) {
    bool flag = false;
    // なんらかの処理
    (flag ? ok : ng) = n;
}

代入処理以外にも使えます.

// 入力された数が奇数ならaへ,偶数ならbへ突っ込む
vector<int> a, b;
for (int i = 0; i < 10; i++) {
    int tmp;
    cin >> tmp;
    (tmp % 2 ? a : b).push_back(tmp);
}

立っているbitの数を数える

bitset を使うと楽チンです.

int b = 31415926;
int c = bitset<32>(b).count();

複数の同じ型の変数に対して同一の処理を施す

範囲for文は initializer_list に対しても使えることを利用します.

vector<int> a, b;

for (auto& e : { a, b }) {
// なんらかの処理
}

おわりに

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() は最初はpqtをタプルに入れて返すようにしてました.最近覚えたばかりの構造化束縛を使ってみたかったのです.しかしながら速さはでませんでした.そこで参照渡しの形に書き換えてやると先程とは打って変わって非常に速くなったのです.また,アルゴリズムを変えると速度が向上するのも非常に興味深いことでした.

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

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

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

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

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

授業が始まって一週間が経った

 なんとか受験をくぐり抜け、入学式を終え、授業が始まって一週間が経った。現時点での所感を書いていこうと思う。

履修登録

 なんで履修登録ってあんなに難しいん?なあ?この一週間は暇さえあれば履修登録について考える日々やったんやけど。…とにもかくにも、履修登録は本当に困難を極めた。必修の授業は何も考えずに入れられるので楽だが、選択必修科目を卒業に必要な単位数に照らし合わせて取るというのは本当に困難な作業だった。それに、年間50単位制限で前期だけでなく後期の授業も考慮しなければならない。そしてなんといっても教養科目という不確定要素。禿げるかと思ったわ。まあ幸いなことに抽選の方はすべて第一希望が通ったので運が良かった。第一希望で通るとは思っていなかったのでそのあたりを調整していたにもかかわらず第一希望で通ってしまったことで逆に狂わされてしまったが、そのことについては気にしないことにする。

 個人的に教養科目の単位は早いうちに取り切ってしまいたいと考えているので、単位取得上限に含まれない時間外科目を取ることにした。おかげで土曜日の午前と交通費9,000円は犠牲になったのだが、安定というものは何事にも代えがたいものだ。これで2回生以降が楽になるのならその程度の犠牲は受け入れる。

 上記のなんやかんやで履修登録を終えた今はひとまず一安心している。あとは出席と提出物とテストをしっかりこなすのみだ。これがうまくできるかどうかは別の話だが…。

課外活動

 ロボコンプロジェクトのForteFibre、大学非公認団体のあくあたん工房に入り、そして競技プログラミング部を設立した。こちらでは本格的な活動はまだ行っていないのでまだどうということもできないが、とりあえずは楽しみにしている。

これからについて

 とりあえずいったん落ち着いたのでバイトに応募した。うまいことバイトと勉強しつつ学生生活を謳歌できたらいいなと考えている。漠然とした不安感はまだ残るが、これから活動していくうちに消えていくものだと信じている。以上。