2018-08-23

分散プログラミングノート #1

勉強日記（バイトの記録）です。分散プログラミングやEthereumの読み物はネット上にたくさんありますが、具体的なアプリケーションのコーディングや冗長な説明が多すぎて読んでいられないので、同じようなことを感じる学習者のために、学習ノートを取っておきます。勉強する際にざっとノートを追ってくださると、余分な学習をしないで済むような構成にしたいと思っています。

8/23

Ethereum入門を全部読みました。

LinuxやOSX向けのgethのインストール手法から、ウォレットの実装まで載っている感じです。とりあえずなにか作ってみたい方は参考にできるかと思いますが、これを読んでも分散プログラミングやコントラクトの本質（理論）はまったくわからないのでhogeでした。プログラムも、JavaScriptの練習みたいな感じでしたね、、、

まあでもgethをまったく触ったことのない方でしたら一度遊んでみることをおすすめします。

CryptoZombiesのLesson 1-3をやりました。こちらも実装よりです。普段からプログラミングに慣れている人にとっては序盤はかなりテンポが悪くて辛いです。こちらに最初の２章の要点がまとまっているので読むだけのほうがいいです。第３章では、もうちょっと広い機能が扱われています。ゲーム形式になっていますが、クリアしなくてもURLを帰るだけで次のステージに進めることに気づきました。

全部合わせて7時間ぐらい勉強しました。

2018-08-07

等式論理の完全性

Infomation Maths

院試の勉強がてら教科書を読み直しているのですが、英語をちゃんと追うモチベが消失しかけそうなのでブログにまとめながら読み進めることによりモチベを回復させようと思います。

完全性の証明がメインなので他の説明は雑めです。また、意図的に数式を減らし言葉で説明しています（数式を読んで理解できる人は教科書だけで足りると思います）。一番いいのは教科書の数式と照らし合わせながら読むことですが、この記事だけでも理解ができるように書いたつもりです。

等式論理とは

たとえば $(x + y)^2=x^2 + y^2 + 2xy$ であることをコンピュータにどう証明させるか、ということを考えます。

話は次のような段階に分けられます。

文法的な点のみに着目して導出木を作る（そのために文法をちゃんと定義する）
導出木で文法的に導出される結果は実は意味的にも正しい（つまり、式の意味を考えず文法的な変換をしているだけなのに間違った式が導出されることがない）
正しい式には必ずそれを導出するような導出木が存在する

まずは1ですね。

Algebraic Specification

Algebraic Specificationとは、等式論理で用いられる演算子と公理をまとめたものです。この公理に加え、反射律や対称律などの推論規則がこの等式論理での導出規則になります（具体的に何があるかなどは省略します）

2の証明は、導出木の高さに関する帰納法（もしくは構成に関する帰納法）を使うことで比較的楽に示すことができます

ここまでの議論は、よく注意してみるとわかるように、式の意味を一切考えずになされています。式の意味というとわかりにくいかもしれませんが、要は $(x + y)^2=x^2 + y^2 + 2xy$ や $x + y = x + y^2$ のように正しい式と誤った式をまだ区別していない、という意味です。

そろそろ意味付けをしてあげましょう。

モデル

モデル（もしくは $\Sigma$ -algebra）とは、変数がとれる値の集合（これをキャリア集合と呼ぶことにする。いまの場合、たとえば実数など）と、演算子の解釈のことです。これまではまだ"+"という記号が足し算という意味をもつとは一言もいわれていませんでした（実は掛け算でした、みたいなことも否定できません）。

あるモデルのもとでの変数の評価とは、変数からキャリア集合への値の対応のことをいいます。具体的には、 $x = 1$ とする、という場合や $x = 2$ とする、という場合は（同じモデルの）別の評価ということになります。

変数の評価ができると、演算子の解釈が定義されていることから、項全体の値をも計算することができます。これを項の評価と呼ぶことにしましょう。

こうして、式が正しいとはどういうことか形式的に定義することができるようになりました。モデル $\mathbb{X}$ のもとである等式 $\mathbb{s} = \mathbb{t}$ が正しいとは、任意の変数の評価の組み合わせに対して、左辺（全体を項とみなしたとき）の評価と右辺の評価が一致する”が成り立つことをいいます。

Algebraic Specificationのもとで正しい式

ある式が、Algebraic Specificationの正しいとは、あらゆる妥当なモデルのもとでその式が正しいことをいいます。妥当なモデル（ $(\Sigma, E)$ -Algebraのこと）とは、Algebraic Specificationの公理にある任意の等式がすべてそのモデルのもとで正しいようなモデルです。

完全性の証明

いよいよ本題です。次のことを証明します。

（目標）導出できないような等式について、ある妥当なモデルが存在して、ある評価が存在し、そのもとでの左辺の評価と右辺の評価が異なる

証明は、実際にその評価が異なるようなモデルを構築することにより行われます。読みやすさのため概略を先に言っておくと、構築するモデルのキャリア集合は、これまでの例で見た実数全体のような集合ではなく、”あらゆる項全体からなる集合”を”項1 = 項2が導出できる”の同値関係で割った同値類、とします。評価は、項からその同値類への写像とします。そうすると導出できない等式について、左辺と右辺は別の同値類に分類されることになり、それがとる値（キャリア集合の要素）は自明に異なるので、左辺と右辺の評価が異なるようなモデルが実際にでき、証明完了ということになります。

騙されたような感じがする証明ですが、その感覚は正しく、実はもっとも重要なポイントは、こうして定義されたモデルが”妥当”であることの証明パートです。それを何も書かなかったため概略だけピンとこないのは当然です。

さて、細かく見ていきましょう。

まず、導出できる等式 $\mathbb{s} = \mathbb{t}$ について、同値関係 $\mathbb{s} \sim_E \mathbb{t}$ を考えます（これが同値関係となる証明は省略）。また、さらに演算子の解釈の定義がこの同値関係についてwell-definedであることを示すことができます。すなわち、演算子の引数をその同値類で置き換えても結果は変わらない、ということです。たとえば、 $(x + y)^2 + z$ の $(x + y)^2$ を $x^2 + y^2$ に置き換えても大丈夫ですよね。

あらゆる項をこの同値関係で割った同値類を考えましょう。各同値類をキャリア集合としたモデルを考えます（演算子の解釈は任意ですが、妥当なモデルであるためにはもちろんいくらかの制約がかかります）。

さて、次のような変数の評価 $J_c$ を考えましょう。すなわち、変数の評価をその変数が属する同値類、ということにします。このとき、項の構成に関する帰納法により、項の評価は、項が属する同値類であることが証明できます。ここで目標をもう一度見てみましょう。

（目標）導出できないような等式について、ある妥当なモデルが存在して、ある評価が存在し、そのもとでの左辺の評価と右辺の評価が異なる

モデルが妥当であると仮定すれば、評価として、 $J_c$ を考えると、導出できない等式の両辺では明らかに対応する値が異なるので目標を達成できます。

さて、こうして定義されたモデルが妥当であることを言えば、上述したとおり、まさにこれが”導出できない等式”の左辺と右辺の評価を（自明に）異なるものにするモデルとなります。復習ですが、モデルが妥当であることを証明するには、定義よりAlgebraic Specificationにある任意の公理がそのモデルのもとで正しいことを言えばよいです。 $J$ を任意の”変数の評価”とします。すなわち、 $J(x)$ は同値類のいずれかをとります。ここで、変数 $x_i$ が属する同値類の代表元を $u_i$ とおくことにします。このとき、左辺の $J$ での評価は、左辺に現れる $x_i$ をすべて $u_i$ で置き換えたもの、を $J_c$ で評価したものとなることが証明できます。代入に関する導出規則により、新たな等式、「公理の左辺に現れる $x_i$ をすべて $u_i$ で置き換えたもの = 右辺に現れる $x_i$ をすべて $u_i$ で置き換えたもの」が成り立ち、右辺に現れる $x_i$ をすべて $u_i$ で置き換えたものの評価は、右辺の $J$ での評価になるので、結局公理にある任意の等式について、左辺の評価と右辺の評価は一致します。つまりいまのモデルは妥当です。

こうして、長くなりましたが、等式論理が完全であることが示されました。

2018-08-01

ARC068 E - Snuke Line

Sport Programming

解説がよくわからなかったので自分用にまとめ直します。

問題

数列 $1, 2, ..., m$ があり、これ上に $n$ 個の区間があります。区間は $[l_i, r_i]$ で表されます。

$m$ 行出力しなさい。 $i = 1, \ldots, m$ 行目には、 $i$ の倍数が何種類の区間に含まれているかを答えなさい。

制約： $1 \le n \le 3 \times 10^5, 1 \le m \le 10^5, 1 \le l_i \le r_i \le m$

ポイント

今 $i$ の倍数を見ているとして、区間の長さが $i$ 以上の区間は確実に含まれます。では区間の長さが $i$ 未満の区間については何種類含まれるでしょうか？これは実は次のように考えても大丈夫です。

長さ $m$ の配列 $cnt$ を用意し、 $j = 1, \ldots m$ 番目を、その区間を長さ $i$ 未満の区間が何種類カバーしているか、を表すことにします。
長さ $i$ 未満の区間を幅が短い順に見ていき（そのような区間のインデックスを $j$ とする）、 $l_j \le k \le r_j$ となるような $k$ について $cnt[k]$ をインクリメントします。この操作はセグメント木やBITを使うことで $O(\log m)$ で可能です。

このとき、 $cnt[i], cnt[2i], \ldots$ の合計が、区間の長さが $i$ 未満の区間のうち、 $i$ の倍数が含まれる区間の種類になります。なぜなら、長さが $i$ 未満の区間は、複数の $cnt[ij], cnt[ik]$ にまたがってカウントされることはありえないからです。もちろん、上記の作業をすべての（長さが $i$ 以上の区間をも含む）区間についてに行ってしまうと重複が発生するのでうまく数えることができません。

さて、このアルゴリズムを愚直に行った場合の計算量は、 $i$ についてのループと、調べる区間の数と、区間についての加算があり、さらに最後に $i$ の倍数を順に見ていく作業があるので $O(nm\log m + n\log m) = O(nm\log m)$ です。まだダメです。

高速化

さて、先程述べたアルゴリズムには明らかに高速化できる点があります。

長さ $i$ 未満の区間を順に見ていき（そのような区間のインデックスを $j$ とする）、 $l_j \le k \le r_j$ となるような $k$ について $cnt[k]$ をインクリメントします。この操作はセグメント木やBITを使うことで $O(\log m)$ で可能です。

の部分で、 $i$ を1から $m$ まで昇順に調べるとすると、 $i$ 番目を調べているとき、 $i - 1$ 番目を調べたときにインクリメントした区間を再びインクリメントしています。これは明らかに無駄であり、 $i - 1$ 番目を調べたときの続きから該当する区間をインクリメントしていくことで、それぞれの区間についてインクリメントされる回数を1回に抑えることができます。こうすると、計算が間に合います。

実装

みやすさのため、セグメント木部分を隠蔽します。

struct range {
    int w, l, r;
    bool operator<(const range& rhs) const {
        return w < rhs.w;
    }
};

int main() {
    int n, m;
    cin >> n >> m;
    
    vector<range> vs;
    for (int i = 0; i < n; i++) {
        int l, r;
        cin >> l >> r;
        vs.push_back({r - l, l, r});
    }
    
    sort(vs.begin(), vs.end());
    SegmentTree cnt(vector<int>(m + 1));
    
    int small_range_begin = 0, small_range_end;
    for (int d = 1; d <= m; d++) {
        int cnt = 0;
        int long_range_cnt = n - int(lower_bound(vs.begin(), vs.end(), range({d, 0, 0})) - vs.begin());
        cnt += long_range_cnt;
        small_range_end = n - long_range_cnt;
        for (int i = small_range_begin; i < small_range_end; i++) {
            cnt.add(vs[i].l, vs[i].r + 1, 1);
        }
        small_range_begin = small_range_end;
        for (int i = d; i <= m; i += d) {
            cnt += cnt.getsum(i, i + 1);
        }
        cout << cnt << endl;
    }
}

2018-07-25

SIMD命令使ってみた

Low Layer Sport Programming

SIMD命令とは

なんか、普通のアセンブリ命令に加えて、128bitの演算とかができる特別な命令（使用できるかはプロセッサに依存しますが、だいたいできそう）

この問題を $O(n^{2})$ で通せると解説に書いてあったので試してみようと思います。

できること

32bitの命令を4回適用する代わりに、128bitの命令を1回適応するようなコードを書きます。 128bitの変数を定義して、4つ分のintをロードし、演算して、その後、4つ分のintの領域にストアしてあげるという流れです。

準備

#include <immintrin.h>

変数や関数

関数名はIntel Intrinsics GuideやSIMD: Integer: Logical operationsを見ればいろいろ載っていますが今回は最小限のものを

__m128i fuga = _mm_set1_epi32(hoge); // 128ビット領域を同じint32_t hogeを4つ並べて初期化
__m128i fuga = _mm_load_si128((__m128i*)&hoge); // hogeのアドレスから128ビット分で初期化（hogeはvector<int>の要素など）

_mm_store_si128((__m128i*)dst, hoge); // 128ビット領域hogeを、dst（配列とか）に戻す

__m128i piyo = _mm_add_epi32(hoge, fuga); // 32ビット整数4つの足し算を同時に
__m128i piyo = _mm_xor_si128(hoge, fuga); // 32ビット4つのxorを同時に（ふつうに128ビットのxorです）

実装

for文のブロック化も込で実装しましたが、なかなかうまくいきません。。。

#include <iostream>
#include <immintrin.h>
using namespace std;

int as[202020];
int bs[202020];
int block = 4000;
int main() {
    int n;
    cin >> n;
    for (int i = 0; i < n; i++) {
        scanf("%d", &as[i]);
    }
    for (int i = 0; i < n; i++) {
        scanf("%d", &bs[i]);
    }

    __m128i tmp = _mm_set1_epi32(0);
    __m128i a, b;
    for (int i = 0; i < n; i+=block) {
        for (int j = 0; j < n; j+=block) {
            for (int ii = i; ii < i + block; ii++) {
                a  = _mm_set1_epi32(as[ii]);
                for (int jj = j; jj < j + block; jj+=4) {
                    b = _mm_load_si128((__m128i*)(bs + jj));
                    b = _mm_add_epi32(a, b);
                    tmp = _mm_xor_si128(tmp, b);
                }
            }
        }
    }

    int ans[4] = {};
    _mm_store_si128((__m128i*)ans, tmp);

    int ret = 0;
    for (int i = 0; i < 4; i++) {
      ret ^= ans[i];
    }

    cout << ret << endl;
}

2018-07-25

ICFPCに参加した（参加してない）

Sport Programming

去年は"ラムダの川渡り"みたいなタイトルの問題で、学科の友達と参加したのですが、たしかデータの送信形式が想定と違っていて修正が間に合わなくて提出失敗、という残念な結果でした。

今年はなんとか提出だけでもACしようと思っていたのですが、提出失敗してしまいました。以下原因。来年は絶対提出して優勝するぞ

unzipコマンドとGUIの展開の違い

簡単に言うと、GUIで展開すると、zip名のフォルダが作られます。しかし、フォルダをzipしたファイルだと、そのフォルダが展開されるだけです。

このせいで提出に失敗しました。提出フォームには、配布ファイルと同じ形式で提出と書いていたのですが、配布ファイルをGUIで展開してしまったため、配布ファイルはフォルダを圧縮したものだと勘違いしてしまいました。

Submit Acknowledgeみたいなページで提出状況を確認できるのですが、不幸にも開始直後に開いてしまい、何もないページだと勘違いしてしまいました（このページは開始２時間後以降にavailableになる）。

まあ実は考察も微妙だったので、今回はhogeみたいなコードを世に知らしめられず良かったとポジティブに考えることにしました。

ICFPC中に同時開催されたFHC Round 1は通過したのでRound 2頑張ろうと思います。

2018-07-14

ICPC2018参加記

Sport Programming

結果としては4完で50位でした。決していいとは言えない成績ですが、個人的には良い大会だったなと思っています。

問題の解法とかは調べたらいくらでも出てくるので、個人的な感想だけを書き残します。

当初のイメージでは、ABを10分程度で片付けて、Cを20分ぐらいで片付けて、そのあと、開始1時間程度で4完して残りの2時間で1,2完して20位ぐらいには入りたいねという感じでした。

しかし、今回のICPCはBを読み始めてすぐ、過去最高難易度のBだ、という感想を持ちました。別に"やるだけ"なんですが、Bのレベルじゃないんですよね。もう終了してから1週間以上経ちますが、思い出せるだけでも

折ると必ずしも厚さが2倍にならない
余分に右にインデックスをとっておかないと右におってセグフォする可能性がある
余分に上にインデックスをとっておかないと同様にセグフォする可能性がある

という点に気をつけないといけなく、最初は初期サイズちょうどの配列をとったためセグフォ、次に、上下左右に余裕をもたせようとした結果、インデックスの管理が大変になって混乱。よく考えてみると右と上だけでいいし、問題文に記載されている座標のとり方に従えば自然な実装ができるので、かなり回り道をしました。

もともと、一つの配列だけを保持して実装していたのですが、この計算量なら折るたびにコピーしても間に合いますし、絶対そっちのほうがミスが減るので早めにそうするべきでした。

自分にはABCぐらいをさっさと通して他のチームメンバーにDEを考察させる余裕をもたせてあげないといけないと思っていたのでなかなかに焦りました。個人戦だと、「まあ今回たまたま実装で詰まったけどまあ気をとりなおしてこ」みたいに思えるのですが、さすがに団体戦だと自分のミスでチームの振る舞いを壊してしまう可能性があるので、かなり責任を感じます。それでも、自分がBをデバッグ中にCを通してくれたのは本当に感謝ですし、さらにはDの実装（大きさhogeの部分集合を列挙していくタイプの問題）までも他メンバーがやってくれたのは本当に本当に感謝です。

練習時点では、5完前後だったのですが、正直このときは、全部の問題を自分が解いて実装しても（時間はもちろん遅れるのですが）同じ問題数解けるんじゃないかと思っていたのであまり全員でやっている感がなかったのですが、今回の予選は本当にチームメンバー全員の力を集結してこの順位だったと思います（ぼくがむしろ足を引っ張ってしまいました）。

次回に向けての反省としては

テストケースがちゃんと書かれた問題文を印刷する（Chromeで印刷するとテストケースが印刷されない謎の仕様があるようです。練習時点で気づいていたのですが本番はプリンターが遠いこともあり妥協してしまいました）
不調時の精神コントロール(⌒,_ゝ⌒)

まあ、競プロがんばります

2018-06-13

Shorの素因数分解アルゴリズム

Infomation Maths

いくつか文献を見ましたが、多分初見で一番わかりやすいのはこれです。

流れを説明すると

Hadamardゲートでレジスタ１について座標変換し $$ | 0\rangle + | 1\rangle + \cdots + | n - 1\rangle $$ という正規直交基底を作り出します。
Unitaryゲートで $$ | 0 \rangle | a^{0} \,\mathrm{mod}\, N \rangle + | 1 \rangle | a^{1} \,\mathrm{mod}\, N \rangle + \cdots + | n - 1 \rangle | a^{n-1} \,\mathrm{mod}\, N \rangle $$ というレジスタ１とレジスタ２の合成状態を作り出します。
2の状態で、レジスタ２の状態は $| a^{k} \,\mathrm{mod}\, N \rangle$ で、これは数通りしかありません。この状態で２を観測することによって状態を固定したとすると、状態１は収縮し、 $r$ を位数として $| a + nr \rangle$ の形のものの重ね合わせに制約されます。しかし、状態１を観測するのは１回しかできないので、なんとかして１度の観測で $r$ を知る必要があります。
そこで量子フーリエ逆変換が登場します。これを行うことで、 $r$ を一発の測定で得ることが可能になります。量子フーリエ変換後のレジスタの状態は $$\frac{1}{N} \sum_{j = 0}^{n - 1}\sum_{k = 0}^{n - 1} \exp(2\pi ijk/n) | k \rangle | a^j \,\mathrm{mod}\, N\rangle$$ となっています。この状態でレジスタ１を観測するわけですが、本来観測結果は $n$ 通りあります。しかし、いま $k = k_0$ が観測されたとして、その確率を考えてみると、そのときのレジスタ２の値（観測していないが収縮しているはず）を $a^{j_0} \,\mathrm{mod}\, N$ とすると $$\Biggl|\frac{1}{N} \sum_{j: x^j \equiv x^{j_0}} \exp(2\pi ijk_0/n) \Biggr|^2$$ となります。和の条件を満たす $j$ が $j = j_0 + mr$ のようにかけるので、これをもとに式を整理すると（省略しますが、誤差が $O(1/N)$ の積分に帰着します。完璧な導出は元論文にあります）、 $|rk_0 \,\mathrm{mod}\, N| \le \frac{r}{2}$ を満たすとき、この値が大きくなり、観測されやすいことがわかります。さて、実はこれを満たす $r$ は１つに絞られます。具体的な操作としては、この不等式を変形すると整数 $d$ を用いて $$-\frac{r}{2} \le rk_0 - dN \le \frac{r}{2}$$ つまり $$\Bigl|\frac{k_0}{N} - \frac{d}{r}\Bigr| \le \frac{1}{2N}$$ なので、 $\frac{k_0}{N}$ を連分数分解することにより $r$ が求められます。