循環的な構造の中で大小関係を推定する

3時より8時の方が後にあり、8時より11時の方が後にありますが、11時より3時の方が後にありそうです。 時計のように巡回している場合、局所的には大小関係のようなものがありそうでも、推移律が成り立つような大小関係はうまく定義できません。

その巡回周期の半分以下の差しか現れない場合のみ、局所的な大小関係が正しく求められます。 時計の例でいえば、差が六時間（十二時間の半分）以内の問題は正しく解けるということです。

特定の期間に含まれるかの判定

ある時刻qが、beginから始まりendを含まずに終わる期間に含まれるかの判定は、以下のように行えます。 問題設定よりbeginよりendが後にあることが保証されるので、期間が巡回周期に近いような場合でも正しく問題が解けます。

constexpr bool includes( size_t begin, size_t end, size_t q ) {
    if (begin > end) {
        return begin <=q || q < end;
    } else {
        return begin <= q && q < end;
    }
}

循環バッファを使ったキューの走査

全然関係ないですが、循環バッファを使ったキューを走査するとき、どこを走査するべきかは、先の問題と同一の解き方ができます。

template<class T, std::size_t N>
class Queue {
    std::array<T, N> buf;
    std::size_t head;
    std::size_t tail;
public:
    Queue() : buf{} head(0), tail(0) {}
    bool empty() const { return tail == head; }
    bool full() const { return (tail+1) % N == head; }
    void push(T elem) { assert( !full() ); buf[tail] = std::move(elem); tail = (tail+1) % N; }
    void pop() { assert( !empty() ); head = (head+1) % N; }
    const T& top() const { assert( !empty() );  return buf[head]; }
    void print() const;
};

こんな感じのキューを想定します。tailは次に要素が来た時に挿入するべき位置、headは最も古く入れた要素の位置です。 なお、この実装だと、空っぽか満杯かを区別できないため、N-1要素しか入れることができません。これを防ぐためには、空っぽであるか（あるいは、満杯であるか）のフラグを用意します。

このとき、このキューに入っている要素を全て書き出すには、以下のようにすればよいです（入れた順）。

template<class T, std::size_t N>
void Queue::print() const {
    if (head < tail) {
        for( std::size_t i = head; i < tail; ++i ) { std::cout << buf[i] << std::endl; }
    } else {
        for( std::size_t i =head; i < N; ++i ) { std::cout << buf[i] << std::endl; }
        for( std::size_t i = 0; i < tail; ++i ) { std::cout << buf[i] << std::endl; }
    }
}

最近使った固定要素数を保持する配列の作り方

よく使うデータへのアクセスを高速化するために、そういったデータを高速な記憶媒体に保持しておくことを、キャッシュすると言ったりします。

高速な記憶媒体はそれ相応に高い*1ので、大容量と両立することはできません。

そこで、よく使うデータを選別して、それだけをキャッシュにおいておくことが全体の高速化につながります。

具体的にどのデータが頻繁にアクセスされるかが事前にわかっている場合は、選別を最適化することができます。 しかし、オンライン的にアクセスがくる場合、未来予知をする必要があるため、最適解を求めることは不可能です。

最適解を近似するアルゴリズムに、LRUアルゴリズムがあります。 これは、キャッシュに入っていないデータへのアクセスが来た時、最も古いデータは捨てて代わりにそれを入れるというものです。

最近使っていないデータへのアクセスは、今後も来ないだろうと考えるのが常識的であり、それを採用したアルゴリズムになっています。 LRUというのは、least recent usedの略で、「最近使っていない」（ものを捨てる）という意味合いになっています。

実際には、一定数のデータを巡回的にアクセスする、などのパターンで最悪の性能を示すという弱点があります。 改善案が日々研究されているようですが、ここでは単純なLRUアルゴリズムを使うことにします。

実際にデータの順序を変更する方法

std::array<T, N> cacheについて、データが来るたびにそのデータを先頭に移動させるという方法です。

if( const auto it = std::find_if( cache.begin(), cache.end(), cond ); it != cache.end() ) {
    // cache hit
    std::rotate( cache.begin(), it, it + 1 );
    return cache[0];
} else {
    // cache miss
    // ....
}

のように、std::rotate関数を使うと非常に簡潔に書けます。

また、頻繁に使われるデータが先頭付近に存在する可能性が高くなるため、先頭からの線形探索は思ったより高速に動作する可能性が高いです（連結リスト - Wikipedia）。

この方法は、sizeof(T)が大きい場合に大量のコピーが発生するために非効率です。

LRUオーダーを持つ方法

std::array<std::pair<T, size_t>, N>あるいは、std::pair<std::array<T, N>, std::array<size_t, N>>のような構造にする方法です。

前者はarray of struct (AoS, 構造体の配列)、後者はstruct of array (SoA, 配列の構造体）になっており、どちらが適切かは微妙なところです。

データの位置は変更せず、追加されたsize_t（実際にはもっと小さな型で十分であることが大半です）に順番（LRUオーダー）を記録します。

  • キャッシュに入っているデータにアクセスが来た場合、現在のLRUオーダーより小さい値のLRUオーダーを一ずつインクリメントし、自分のLRUオーダーを0にする
  • キャッシュに入っていないデータにアクセスが来た場合、LRUオーダーがN-1である要素を新たなデータで上書きし、LRUオーダーを0にする。他のLRUオーダーは一ずつインクリメントする

データの移動がないため、参照が無効になることがないのはうれしい点です。 どれを捨てるのかの判断にも線形探索が必要な点は「実際にデータの順序を変更する方法」に劣ります。

CPUのキャッシュは普通この方法で実装されています。

タイムスタンプを持つ方法

「LRUオーダーを持つ方法」では、アクセスのたびにLRUオーダーの書き換えが大量発生していましたが、代わりにタイムスタンプを持つ方法です。

アクセスがあるたびにタイムスタンプを更新し、どれを捨てるかの判断は単なる最小値の線形探索をすることによります。

タイムスタンプがオーバーフローしないようにするため、タイムスタンプ領域のサイズは「LRUオーダーを持つ方法」より大きくなります。 そのため、sizeof(T)が小さい場合にはそのオーバーヘッドが無視できなくなります。

以上の三手法は、sizeof(T)の大きさによって使い分けるとよいです。

sizeof(T)が小さい場合は「実際にデータの順序を変更する方法」が簡単です。一方、タイムスタンプ領域が無視できるほどsizeof(T)が大きいのであれば「タイムスタンプを持つ方法」が便利です。

なお、「タイムスタンプを持つ方法」は、タイムスタンプが最も古いものがどこにいるのかを判別するために最小ヒープを用いることで計算量を改善することが可能です。

ところで、「実際にデータの順序を変更する方法」は追加の領域が必要になっておらず一見不思議です。

これは、集合としておなじ配列に内在する冗長性（順番成分）がlog_2(N)ビット存在することが原因です。 LRUオーダーを覚えるためにこの情報量をかつようしているために追加の領域が必要になっていないのです。

連想配列に適用する方法

先の3つの手法は、キャッシュにデータがあるかの探索は線形検索でした。 しかし、通常の用途では連想配列を用いることが多いです。

さきの3つの手法のうち、「LRUオーダーを持つ方法」と「タイムスタンプを持つ方法」はそのまま連想配列に適用できます。 「実際にデータの順序を変更する方法」は連想配列の場合には使えません。

ところで、ECMAScript 2015で追加されたMapは、挿入順でイテレートするという動作をします。

これを利用すると、追加の記憶なしに「最近使ったことにする」touch（LRUオーダーを0にする）や「古いのを捨てて挿入する」insert関数が以下のように書けます。

性能がベストなのかはよくわかりませんが、それほど悪くないように思われます。

class FullAssociativeLRUTable extends Map {
    constructor(size) {
        super();
        this.maxSize_ = size;
    }

    touch(tag) {
        if (this.has(tag)) {
            // いったん削除して再挿入で、MRU position に移動
            let data = this.get(tag);
            this.delete(tag);
            this.set(tag, data);
        } else {
            console.log(`${tag}は含まれていないのにtouchした`);
        }
    }

    insert(tag, data) {
        if (this.has(tag)) {
            console.log(`${tag}は含まれているのにinsertした`);
        } else if (this.size < this.maxSize_) {
            this.set(tag, data);
            return null;
        } else {
            // keys().next()は最初の要素 (LRU position)
            let minTag = this.keys().next().value;
            let data = this.get(minTag);
            this.delete(minTag);
            this.set(tag, data);
            return data;
        }
    }
}

*1:パレート最適な記憶媒体を選んだら当然そうなります

累積和を使っても二分探索を使っても解ける問題

先週話題にした、ABC(AtCoder Beginner Contest)122のC問題は、累積和を用いて解くのが想定解のようですが、二分探索を用いても解けるという話を、少し一般化して考えてみます。

なお、以下では何も言わずに空間計算量だけ書いた場合、時間計算量も同じオーダーであるとします（それだけの量の空間を触るのに最低限必要な時間です）。

Θ(1)時間で計算可能な特性関数 \mathbb{N}→{0,1}が入力長Ο(N)で与えられる場合

元問題は、このクラスです。 元問題では、特性関数の種として長さNの文字列が与えられます。

累積和

  1. 【前処理】各点での特性関数の値のprefix sumになっている配列を構成します。
  2. 配列アクセスにより、端点での値を得ます。
  3. その差が区間内にある特性関数が1になる点の集合のサイズとなります。

この累積和を使った解法では、

  • 前処理にΘ(Nlog N)空間*1
  • クエリ一つあたりにΘ(log N)時間*2

で解けます。

二分探索を用いる解法

  1. 【前処理】特性関数が1になる点の集合を昇順ソート済み配列として構成します。
  2. 端点を二分探索します。下端（始端）はlower_boundのindex、上端（終端）はupper_boundのindexを得ます。
  3. その差が区間内にある特性関数が1になる点の集合のサイズとなります。

この二分探索を使った解法では、

  • 前処理にΘ(Nlog N)空間
  • クエリ一つあたりにΘ(log N)時間

かかります。後者の解法は定数項が大きく、特に優れた点は見受けられませんが、累積和手法を思いつけなかった場合には有用そうです。

各点での重みが0, 1とは限らない自然数の集合Sであり、それをΘ(1)時間で計算可能な \mathbb{N}→Sが入力長Ο(N)で与えられる場合

Sの要素に何が含まれるかが重要となります。各点での重みのとり得る集合Sの要素で最大のものをKとします。

この場合、累積和手法では何も複雑なことはありません。なお、空間計算量はlog K倍悪化します。

二分探索手法では、要素の重複を許す集合として構成すれば、同じ手法が使えます。しかし、前処理の空間計算量が最悪K倍悪化し、クエリ一つあたりにかかる時間もΟ(log K)倍悪化します。

つまり、KがΘ(N)くらいのオーダーであると、前処理の空間計算量がΘ(N2 log N)となり、N=105の時に2秒以内には完了しないアルゴリズムとなります。

特性関数が、特性関数が1になる点（N以下の自然数）のソートされていない列（長さK）として与えられる場合（K<N）

この場合に累積和手法を適用するには、前処理部分に以下の二つの考え方があります。

  1. ソートしてから累積和。
  2. いもす法。

1.の手法の場合、累積和を取る以上、Θ(Nlog N)の空間の確保は避けられません。 Kとしてあり得るのはN未満の自然数なのでソートはバケットソートとするのが計算量のオーダーとして最善になります。

結局両者のやる仕事は同一になります。

一方、二分探索の手法では、前処理として単にΘ(Klog K)（厳密には、さらにΟ(log N)×Θ(log K)倍の時間がかかるはずです）*3時間かかるソートを行うだけで済みます。

つまり、Kの大きさによっては、二分探索手法が生きる局面があり得ます。それは、KがNに比べて十分に小さい場合です。

特性関数が1となる点が自然数とは限らない（doubleで表されたり、有理数となる）場合

この場合には、累積和手法を適用するのは困難となります。一方、二分探索手法はそのまま適用可能です。

これは、Nが非常に大きく、Kがそれに比べて小さいパターンと同一視できます。

各点での重みが0, 1とは限らない自然数の集合Sであり、その点も自然数とは限らない場合

どちらの手法でも困難に思われます。セグメント木を使うと解けそうです。

Sが任意の有理数である場合も同様だと思います。

まとめ

  • 重みが非零となりうる点の数が少なければ、累積和手法は有効。
  • 各点の重みのとり得る値が非負で小さければ、二分探索手法は有効。
  • 元問題はその両方の性質「重みが非零となりうる点の数はΟ(N)」「各店の重みのとり得る値は0か1」を満たしていたためどちらでも解けた。

*1:通常考慮外とされますが、N未満の自然数を表すにはlogNビット必要な点に注意してください。

*2:N未満の自然数同士の引き算にはΘ(log N)時間かかる点に注意してください。

*3:Ο(log N)倍は、N未満の自然数（Nビットで表せる）の比較にかかる時間（入力がランダムなら平均Θ(1)ですが、最悪の場合はΘ(log N)です）、Θ(log K)項は、要素へのポインタの移動時に書き換えないといけないビット数です（要素の移動だとΘ(log N)ビットの書き換えが必要ですが、要素数が少ない場合要素へのポインタの移動にしてしまえばΘ(log K)ビットで済みます）。

AtCoder Beginner Contest 122に参加した記録

忘れないうちに書いておきましょう、と思って既に二週間が経過しました……。

この日はちょうど私の所属しているサークルの合宿の日で、サークルのメンバーのうち何人かがAtCoder Beginner Contest (ABC) に参加するという感じで盛り上がっていました。

私はコードを書くのをせかされる感じが嫌なのでコンテストは嫌いだったのですが、まぁ食わず嫌いなのもよくないので参加することにした感じです。

一応、数年前に、競技プログラミングの過去問を解く授業を受けたことはあったし、その授業の教科書（参考書だったかな）として買った蟻本はすべて読んだので、このくらいのレベルの問題は解けないとまずいなぁという思いもありつつ。

（以下、時系列で書いていきますが、少なくともコンテストに関しては初心者中の初心者なので、参考になるところはないと思われます……）

18:30～19:00

サークルメンバーのみんなで夜ご飯を食べました。ABCがあるので参加しようという話が出ました。

20:00～20:45

二分探索難しいとかいう話題で盛り上がりました。C#は書いたことがあるがC++は書いたことのないというメンバーに、最低限の知識（std::cin、std::cout、std::vector、std::string）を教えました。

20:55

やっぱり参加しようと決断し、ノートパソコンを取り出したものの、ファンが壊れていて異音がします……。うるさいので別室で受験することにしましょう。一時間しか寝ていないのでねむい……。

21:00

ABCスタート。とりあえずD問題をちらっと見てDPっぽいなぁと思いました。

21:01

A問題を解きましょう。特に難しいところはないようですが、if文を書くのはめんどくさいなぁとか、条件演算子を書くとみづらくなるなぁとか、そういうことを考えてやや時間を無駄にした感じがあります。 switch文を使ったコードを書いて正答になりました（3:04）。

なお、sed y/ATGC/TACG/でよいということに気付いたのは、コンテスト終了後でした（後述）。

21:03

B問題を解きましょう。制約が非常に緩いので、全ての部分文字列を列挙する手法が可能だなぁと思いました。この手法は、配列インデックス境界バグ*1を生みにくいので制約が緩い以上こっちの方針でやろうと思いました。実際、その方が読みやすく、つまりバグを埋めにくいようなコードになるはずです。

……と思ったら、C++のstd::string::substr関数は（何文字目から、何文字分）みたいな指定方法なので結局境界値バグを生みやすい形ですね……。

冷静に考えると、Θ(N)解法があることに気付きました。尺取り法（の退化したバージョン？）を使えばよさそうです。

私の実力では、このあたりになってくると頭の中だけでプログラミングするのは無理です。特に準備もしてこなかったので環境が整っているわけでもなく、wandboxでコンパイルして確かめることにしました。

とりあえずテストケースに成功したので提出してみます（8:52）。

21:08

C問題を解きましょう。クエリですか、なるほどです。クエリといえば、前処理です。これは累積和でできそうですね。 累積和でやるとすると、落とし穴になりそうなところ（コーナーケース）は、

  • 配列境界（-1要素目とか、N要素目とかにアクセスしてはいけない）
    • 'AC'という二文字を考えるので、str[i-1]とかstr[i+1]で失敗しないように気を付ける必要があります
  • 文字列の切れ目がちょうど'AC'にまたがる場合
    • 単純な累積和よりミスを誘うポイントが増えています
  • その複合
    • 文字列の切れ目が、ちょうど1文字目とか、末尾とか、そういうところになった場合、上記二つ同時に気を付けないといけません
    • （番兵を置いておけばこの問題はなくなりますが、このことに気付いたのはコンテスト終了後でした）

あたりでしょうか。

21:10

あれー、B問題がWAになっています。ちらっとコードを見た感じではバグっている箇所を見抜けませんでした。テストケースの後半が全滅しているのもよくわかりませんが……。とりあえず放置してC問題を解く作業に戻りましょう。

クエリの形式が1-originであることとか、いろいろ混乱し悪戦苦闘した結果、30分弱時間を溶かしましたが、なんとか正答です（35:10）。

サンプルケースがコーナーケースを網羅していたのはやさしさでしょうか。

ちなみに、紙とペンを用意していなかったので、以下のようなアスキーアートをメモ帳（notepad.exe）に描いて考えていました。

 A C A C T A C G 
| | | | | | | | |
0 0 1 1 2 2 2 3 3

21:35

B問題の提出コードをもう一度見てみましょう。文字列終端に達した時に、尺取虫のしっぽを引っ込める作業を忘れていることに気付きました。時間を置くと見つかるものですね。正答です（36:44）。

21:36

さてあと一時間ほどしかありませんが、D問題を解きましょう。【違反してない文字列】の末尾に文字をくっつけた時、【違反】が発生するかを考えればよさそうです。なぜなら、【違反している文字列】の末尾に何かを付け加えても、違反していない状態に戻るわけではないからです。よく考察していませんが、たぶん末尾3文字が「状態」であるDPでうまくいくのではないでしょうか。コピペで以下のコードを作り出しました。

std::uint64_t arr[4][4][4][101]
for( int i = 0; i < 4; ++i )
for( int j = 0; j < 4; ++j )
for( int k = 0; j < 4; ++k )
     if( i == 0 && j == 2 && k == 1 ) arr[i][j][k][3] = 0;
else if( i == 0 && j == 1 && k == 2 ) arr[i][j][k][3] = 0;
else if( i == 2 && j == 0 && k == 1 ) arr[i][j][k][3] = 0;
else arr[i][j][k][3] = 1;

std::size_t N;
std::cin >> N;
for( std::size_t t = 4; t <= N; ++t )
for( int i = 0; i < 4; ++i )
for( int j = 0; j < 4; ++j )
for( int k = 0; j < 4; ++k )
     if( i == 0 && j == 2 && k == 1 ) arr[i][j][k][t] = 0;
else if( i == 0 && j == 1 && k == 2 ) arr[i][j][k][t] = 0;
else if( i == 2 && j == 0 && k == 1 ) arr[i][j][k][t] = 0;
else arr[i][j][k][t] = (arr[i][j][0][t-1] + arr[i][j][1][t-1] + arr[i][j][2][t-1] + arr[i][j][k][3][t-1]) % 1'000'000'007;

std::uint64_t sum = 0;
for( int i = 0; i < 4; ++i )
for( int j = 0; j < 4; ++j )
for( int k = 0; j < 4; ++k )
sum = (sum+arr[i][j][k][N]) % 1'000'000'007;

エレガントさのかけらもありませんが、こちらも必死なので少なくとも自分にはぎりぎり読める感じの、わけのわからないコードになってきました。

いちおう、ケチらずにstd::uint64_tを選択したことにより、貰うDPとして書いてもオーバーフローの心配がないのは良かった点ですが、やっぱり配るDPで書いた方がよかったかなぁと思いつつ、しかし貰うDPとして書いた方が直感的だし(?)こっちでよかったかなとか思っていました。

このコードで3を入力すると61が出るのは当然ですが、4を入力すると235が出てきてテストケースと合いません。つまり、何か漏らしている条件があるのですが、まぁおおむね方針は立ちました。

22:01

arr[i][j][k] = (arr[0][i][j][t-1] + arr[1][i][j][t-1] + arr[2][i][j][t-1] + arr[3][i][j][k][t-1]) % 1'000'000'007;

の方が正しいことに気付きました。 233って出てきました。依然、何か見落としています。

22:04

"xxxxAGx"のような文字列の末尾にCを付け加えると条件に違反することに気付けました。

else if( i == 2 &&  k == 1 ) arr[i][j][k][t] = (arr[1][i][j][t-1] + arr[2][i][j][t-1] + arr[3][i][j][t-1]) %  1'000'000'007;

231って出てきました。もう一条件あるのを見落としていそうです。時間間に合わないかもなぁ、とか思いつつ、N=4のケースで違反する26通りを手で列挙するという作業を地道にこなすことにします（というかコード上で考えるのではなく、最初からそうするべきでした）。

22:18

末尾三文字が"AGC"、"ACG"、"GAC"になっている場合（"AGC"の原因が固まって存在している場合）のほかで違反してしまうパターンがなんとなくわかりました。 つまり、末尾四文字の中に"AGC"の原因が、関係ない一文字を挟んで存在している場合（"AGxC"と"AxGC"）というパターンでした。

else if( j == 2 &&  k == 1 ) arr[i][j][k][t] = (arr[1][i][j][t-1] + arr[2][i][j][t-1] + arr[3][i][j][t-1]) %  1'000'000'007;

ちゃんと230って出ました。N=100のケースを試してみたところ、ぴったり一致したため、ほぼあっている確信を得て提出、正答でした（79:28）。

22:19

この時点では、287位だったっぽいです。知り合いの順位を調べたり、Twitterをしたりして、暇をつぶします。

22:30ごろ

ふと順位表を見たら300位になっていました。なるほど、WAが一回あったから、WAなし80分完答みたいな参加者に抜かされたわけですね。

22:40

コンテストおわり。サークルメンバーと雑談します。

感想戦

D問題（DPのやつ）

  • やはり貰うDPではなく配るDPの方が筋がよかったらしい
  • あと、初期化のために「文字なし」っていう状態を付け加えた、53状態のDPをやるときれいに書けたらしい

C問題（クエリのやつ）

  • あるメンバーは累積和は思いつかなかったものの、（事前に二分探索の話をしていた影響もあって）二分探索で解く方法にたどり着いたらしい
    • 二分探索でも解けるということ自体が興味深い話に思えるので、後日記事を分けて書こうと思います
  • 累積和でどこを基準にするかみたいな話は、自分の中で何かルールを決めておくと、コンテストの限られた時間で思考をまとめやすいという話を伺いました

A問題（簡単なやつ）

  • こういうのはesolangで解く参加者もいるよね、あー普通にsedで解けるじゃん、みたいな話をしました
    • とはいっても、初めてのコンテストで緊張している中、勝手もわからないesolangを使う気にはなれないような……

コンテスト参加者のC問題の正答例を見て

どうやら、ACの区切れ目を気にするのではなく、Aが直前に存在するCっていう風にとらえてカウントしていけばよかったようです（解説ではCが直後に存在するAって書いてありますが、同じことです）。

さすがに配列境界の問題に警戒しすぎたっていうことなのかな？

反省

見通しを立てるのは数秒でできているのに実装に何十分もかかっているあたり、やはり細かいところを詰めるのが苦手だなぁと改めて実感しました。

B問題でのWAのようなものは、単なるケアレスミスと甘く見ず、本当に意識して練習しないとなおらないということを大学受験で経験している*2ので、ちゃんと反省したいです。

準備不足は、いつでも完璧なコンディションでできるわけではないので、それほど悔いていませんが、最低限紙と鉛筆はあった方がよさそうということは学びました。

結果

  • AtCoder Beginner Contest 122
  • 300th / 3024
  • パフォーマンス 1600
  • レーティング 400 (初参加)
  • 段級位 9級

*1:N要素の配列の（0-originで）N番目にアクセスしてしまうみたいなバグ。私はこのバグを発生がちなので、警戒していました

*2:化学の計算問題で計算ミスを多発させていましたが、「本質はわかっているのだから気を付ければ問題ない」と甘くスルーしてしまっていました。幸い、受験までの残り時間が十分にある時期に、この考えが誤りであることに気付いてなおせたので事なきを得ましたが。

近況とか今年度の反省とか

今週は無理しすぎて風邪をひいてしまいました。つらい。

来年度からはもっと規則正しい生活をするようにしたいです。

それから、「いつまでにやるか明確に決まっていないけれど、いずれやらないといけないこと」は3月に片づけようと思っていたのに、いろいろ自分に言い訳して結局片付きませんでした。 こういうのを先延ばしにする癖を直したいです。

そういえば、今年度はなんだか、Twitterでは知っているけれどあったことのない人間にたくさん会った一年間でした。 そういう機会を増やしていきたいですね。

最大公約数をもっと高速に求める（その4）

先週の記事の続きです。

最大公約数をもっと高速に求める（その3） - よーる

前回示したコードは、以下のようなものでした。

uint64_t gcd_impl( uint64_t n, uint64_t m ) {
        constexpr uint64_t K = 5;
        for( int i = 0; i < 80; ++i ) {
                uint64_t t = n - m;
                uint64_t s = n - m * K;
                bool q = t < m;
                bool p = t < m * K;
                n = q ? m : t;
                m = q ? t : m;
                if( m == 0 ) { return n; }
                n = p ? n : s;
        }
        return gcd_impl( m, n % m );
}

uint64_t gcd_pre( uint64_t n, uint64_t m ) {
    for( int i = 0; i < 4; ++i ) {
        uint64_t t = n - m;
        bool q = t < m;
        n = q ? m : t;
        m = q ? t : m;
        if( m == 0 ) { return n; }
    }
    return gcd_impl( n, m );
}

uint64_t gcd( uint64_t n, uint64_t m ) {
        return n > m ? gcd_pre( n, m ) : gcd_pre( m, n );
}

分岐を外していい場合

（その1）で示したコードでは、互減法と互除法のどちらか適切な方を選ぶ分岐をするより、10回の固定回ループの方が速いことを使って高速化しました。 これは、分岐予測ミスのペナルティ時間で10回くらいのループは実行できるということによっています。 ちなみに、10回というのをちゃんと最適化すると、8回か9回がベストなようです。

今回は、80回の固定回ループとなっているため、それと同じ論理では固定回数ループの方が速いとは言えません。 つまり、互減法と互除法のどちらか適切な方を選ぶ分岐をすることで高速化する可能性が残されています。

nがmの何倍大きいときに互除法が適切だと判断するべきかという問題は、実験してみると32倍あたりに最適解があるようです。 これはループ6回分に相当します。

（その1）で示したコードでは、ループ一周当たり3cycleで8周がベストという結果でした。つまり、分岐ミスペナルティと引き換えに無駄な仕事をする価値があるのは24cycleくらいだということになります。 今回のコードでは、ループ一周当たり4cycleで6周がベストという結果だったので、やはり24cycleくらいだということが支持されます*1

そこで、以下のようなコードとしてみます。

uint64_t gcd_impl( uint64_t n, uint64_t m ) {
    constexpr uint64_t K = 5
    if( m == 0 ) { return n; }
    for( ; n / 32 < m; ) {
        uint64_t t = n - m;
        uint64_t s = n - m * K;
        bool q = t < m;
        bool p = t < m * K;
        n = q ? m : t;
        m = q ? t : m;
        if( m == 0 ) { return n; }
        n = p ? n : s;
    }
    return gcd_impl( m, n % m );
}

このようにすることで、さらに一割ほど高速化しました。

今週はちょっと忙しくて記事を書く時間が取れなかったのでちょっと少ないですが、ここまで！

ここまでで、std::gcd関数の2.4倍ほど速いgcd関数を作れています。

*1:もちろん、分岐ミスペナルティの確率が違う＆引き換えにできる仕事の内容が違うため、単純な比較はできないのですが、20-30cycleというのは一つの目安になります。

最大公約数をもっと高速に求める（その3）

先週の記事の続きです。

最大公約数をもっと高速に求める（その2）【cmova命令は遅い】 - よーる

前回示したコードは、以下のようなものでした。

uint64_t gcd_impl( uint64_t n, uint64_t m ) {
        for( int i = 0; i < 10; ++i ) {
                uint64_t t = n - m;
                bool q = t < m;
                n = q ? m : t;
                m = q ? t : m;
                if( m == 0 ) { return n; }
        }
        return gcd_impl( m, n % m );
}

uint64_t gcd( uint64_t n, uint64_t m ) {
        return n > m ? gcd_impl( n, m ) : gcd_impl( m, n );
}

このコードの機械語命令列

gcd_impl関数のループ一回分の機械語命令列は、以下のようになっています（実際にはループアンローリングされており、レジスタ割り付けだけが異なる命令列が3回周期で現れます）。

sub    rdi, rsi    # t = n - m;
cmp    rdi, rsi    # q = t < m;
mov    rax, rdi
cmovb  rax, rsi    # n = q ? m : t;
cmovnb rdi, rsi    # m = q ? t : m;
test   rdi, rdi    # m == 0
jz     Label       # if( m == 0 )

ループ一回分は7命令になっているわけですが、だからといって7 cycleかかるわけではありません。近年のCPUでは、これくらいの命令列は並列に実行されます。

今回の命令列で、並列実行のボトルネックになるのは依存の連鎖です。 具体的には、t = n - mを計算し、その結果を使ってq = t < mを計算し、その結果によって次ループのnとmが算出される部分です。

qはtが算出されなければ計算できませんし、nとmはqが算出されなければ計算できません。また、次のループのtはnとmが算出されなければ計算できません。

つまり、どんなに並列に実行できるCPUを持ってきても、この部分は各々1cycleかかってしまいます。

よってこのコードでは、ループ一回分に3cycleかかることになります。

なお、mov命令が気になるかもしれませんが、最近のCPUでは実質的に0cycleで実行されると考えてよいようです。

互減法をさらに高速化する

互減法の一ステップでは、以下のように計算が進みます。

GCD(n, m) = if n-m >= m then GCD(n-m, m) else GCD(m, n-m)

これを拡張し、以下のようにしてみます。

GCD(n, m) = if n-Km >= m then GCD(n-Km, m) else if n-m >= m then GCD(n-m, m) else GCD(m, n-m)

もしn-Km >= mであれば、互減法のKステップを一気に進めることができるため、互減法の効率が上がることが予想されます。 また、それにより、アルゴリズム全体に占める除算の割合が減少し、高速化すると予想されます。 ただし、以下の点に注意する必要があります。

  1. 途中結果n-Kmが負になる可能性がある。
  2. Kmがオーバーフローする可能性がある。

1.については、移項することで回避可能です。 2.については、前処理でmを小さくしておけば回避可能です。

改良されたメインループ

メインループの内部を以下のようにしてみます。メインループの前に前処理を置くことでうまくmは小さくしているとして、m*Kはオーバーフローしないと仮定します。

uint64_t t = n - m;
uint64_t s = n - m * K;
bool q = t < m;
bool p = t < m * K;
n = q ? m : t;
m = q ? t : m;
if( m == 0 ) { return n; }
n = p ? n : s;

n = p ? n : sが追加されたため、さらに依存の連鎖が伸び、このループ一回には最低でも4cycleかかるようになりました。 もし、m*Kの計算に4cycle以上かかる場合、そのせいで依存の連鎖が伸び、ループ一回にかかるサイクル数はもっと長くなります。 あるいは、m*Kの計算が他の命令の実行の裏で並列に実行できない場合も、ループ一回にかかるサイクル数はもっと長くなります。

最近のIntel社製CPUであれば、MUL命令は3cycleで完了するため、m*Kの計算のせいで依存の連鎖が伸びるということはありません（Kの値によってはコンパイラ最適化によりLEA命令を駆使した命令列となり、もっと速く計算されます）。 また、m*Kの計算は、他の命令の実行の裏で行うことができるため、そのせいでループ一回にかかるサイクル数が長くなるということはありません。

元のコードのループ一回には3cycleかかっていたところ、4cycleかけることで、Kステップ一気に実行できる可能性が追加されたことになります。 nとmの大きさに大きな差がなかった場合でも、少なくとも1ステップは進むため、それほど損したことにはなりません。

Kとして何をとるか

Kとして2や3などの小さな数をとれば、n-Km >= mである可能性がそれなりに高いため、改良により互減法のステップを効率よく進めることができるようになります。 Kとして大きな数をとる場合、n-Km >= mである可能性が低く、意味のある仕事をできる確率（pがfalseになり、n = sが実行される確率）が下がります。しかし、意味のある仕事をできた場合、互減法のステップを一気に進めることができます。

このトレードオフを考慮に入れ、いくつかの実験を行った結果、以下にようにすることでstd::gcdより二倍以上高速なgcd関数を作れることがわかりました。

uint64_t gcd_impl( uint64_t n, uint64_t m ) {
        constexpr uint64_t K = 5;
        for( int i = 0; i < 80; ++i ) {
                uint64_t t = n - m;
                uint64_t s = n - m * K;
                bool q = t < m;
                bool p = t < m * K;
                n = q ? m : t;
                m = q ? t : m;
                if( m == 0 ) { return n; }
                n = p ? n : s;
        }
        return gcd_impl( m, n % m );
}

uint64_t gcd_pre( uint64_t n, uint64_t m ) {
    for( int i = 0; i < 4; ++i ) {
        uint64_t t = n - m;
        bool q = t < m;
        n = q ? m : t;
        m = q ? t : m;
        if( m == 0 ) { return n; }
    }
    return gcd_impl( n, m );
}

uint64_t gcd( uint64_t n, uint64_t m ) {
        return n > m ? gcd_pre( n, m ) : gcd_pre( m, n );
}

gcd_pre関数は、m*Kがオーバーフローしないようにするための前処理です。互減法を4回繰り返すことで、mは大きくとも0xffff'ffff'ffff'ffff / 5であると保証できます。

Kとして5をとったことにより、32のような大きな商が立つ場合でも、6回程度ループを回れば互除法の一ステップを行ったことになります。これ以上大きな商が立つ場合は、除算を行った方がよさそうです。 そのようなことの発生頻度が低いため、単純な互減法であれば10回程度で除算に切り替えたほうがよかったところ、80回まで増やすことが可能になりました。 なお、gcd_implのループ回数80回は目安であり、必ずしも最適な値とは限りません。