SuperCon 2022 予選参加記

はじめに

3 回目のSuperConです. 今回の問題はここにあります.

https://www.gsic.titech.ac.jp/supercon/main/attwiki/index.php?plugin=attach&refer=SupercomputingContest2022&openfile=SuperCon2022Yosen20220606B2033.pdf

問題概要

 11 \times 11 のマス目それぞれに  0 から  9 までの数字が書き込まれており,座標  (i, j) の数字は  C _ {i, j} である.

時刻  0 には座標  (5, 5) にロボットがあって,あなたは時刻が  1 進むことに上下左右隣接するマスにロボットを進めることができる.(時刻  t でのロボットの位置を  (x _ t, y _ t) とする.)

各時刻においての目標の色のパターン  Cr _ t が与えられるので,ロボットの進路を好きに決め,コスト関数  \sum _ t \min (|C r _ t - C _ {x_t, y_t} |, 10 - |C r _ t - C _ {x_t, y_t} |) を最小化せよ.

書いたアルゴリズム

DP をして復元する方針とダイクストラ (または 0-1BFS っぽいもの) をする方針があると思いますが,DPの方針を取りました. 単純に時刻とその時のロボットの場所をキーにしてそこまでのコスト関数を持てばよいです.

ただこのままだとおもしろくなかったので,SIMD 化して,各行の処理を並列して実行するようにしました. 問題文の実行環境の欄にある「2.8 GHz Intel Core i5」という情報だけだとどの世代かがわからなかったので,sse4.2 相当対応の CPU で動くようにしています.

また, (i, j) \rightarrow (i+j, j) のように座標をいい感じに変換することで,やや計算量が下がります.(各行  \min の回数が  1 回減る.)

座標変換のイメージ

実行時間

手元で  1 ケースあたり平均  2.2 \mu 秒でした.(標準偏差  0.45 \mu 秒) 愚直な DP 解が  1 ケースあたり  110 \mu 秒だったのでかなり高速だと思います.

提出コード

#include <immintrin.h>

#include <algorithm>

#include "sc1.h"

inline int to_index(int x, int y) { return (x * SC_L + y); }

inline __m128i get_rotate_perm(int n, int r, int l) {
    alignas(32) unsigned char a[l];
    for (int i = 0; i < n; i++) {
        a[i] = (i + r + n) % n;
    }
    for (int i = n; i < l; i++) {
        a[i] = i;
    }
    return _mm_load_si128(reinterpret_cast<__m128i*>(a));
}

inline int mod_l(int x) {
    if (x >= SC_L) {
        return x - SC_L;
    } else {
        return x;
    }
}

__attribute__((target("sse4.2"), optimize("unroll-loops"))) int main() {
    SC_input();

    constexpr int SIMD_L = 16;

    __m128i rotate_right = get_rotate_perm(SC_L, -1, SIMD_L), rotate_left = get_rotate_perm(SC_L, 1, SIMD_L);

    alignas(32) unsigned char dp[SC_Nt + 1][SC_L][SIMD_L];
    std::fill(dp[0][0], dp[0][SC_L], 127);

    constexpr int center = (SC_L - 1) / 2;
    dp[0][0][center] = 0;

    __m128i C[SC_L];
    for (int i = 0; i < SC_L; ++i) {
        alignas(32) unsigned char tmp[SIMD_L];
        for (int j = 0; j < SC_L; j++) {
            tmp[j] = SC_C[to_index(mod_l(i + j), j)];
        }

        C[i] = _mm_load_si128(reinterpret_cast<__m128i*>(tmp));
    }

    for (int i = 0; i < SC_Nt; i++) {
        __m128i min[SC_L];
        for (int j = 0; j < SC_L; j++) {
            min[j] = _mm_load_si128(reinterpret_cast<__m128i*>(dp[i][j]));
        }

        for (int j = 0; j < SC_L; j++) {
            min[j] = _mm_min_epu8(min[j], _mm_shuffle_epi8(min[j], rotate_left));
        }

        for (int j = 0; j < SC_L; j++) {
            __m128i abs = _mm_abs_epi8(_mm_sub_epi8(C[j], _mm_set1_epi8(SC_Cr[i + 1])));
            __m128i dist = _mm_min_epu8(abs, _mm_sub_epi8(_mm_set1_epi8(SC_Nq), abs));
            __m128i res = _mm_min_epu8(min[mod_l(j - 1 + SC_L)], _mm_shuffle_epi8(min[mod_l(j + 1)], rotate_right));
            res = _mm_add_epi8(res, dist);
            _mm_store_si128(reinterpret_cast<__m128i*>(dp[i + 1][j]), res);
        }
    }

    int min_x = 0, min_y = 0, min = 255;
    for (int i = 0; i < SC_L; i++) {
        for (int j = 0; j < SC_L; j++) {
            if (dp[SC_Nt][i][j] < min) {
                min_x = i;
                min_y = j;
                min = dp[SC_Nt][i][j];
            }
        }
    }

    unsigned char next_min = 255;
    for (int i = SC_Nt - 1; i >= 0; i--) {
        SC_ans[i + 1] = to_index(mod_l(min_x + min_y), min_y);

        int next_x, next_y;

        for (int s = 0; s < 2; s++) {
            int x = mod_l(min_x + (s == 0 ? SC_L - 1 : 1));
            for (int j = 0; j < 2; j++) {
                int y = mod_l(min_y + (s == 0 ? j : SC_L - j));
                if (next_min >= dp[i][x][y]) {
                    next_x = x;
                    next_y = y;
                    next_min = dp[i][x][y];
                }
            }
        }

        min_x = next_x;
        min_y = next_y;
    }

    SC_output();
}

Barrett Reduction について考えたこと

Barrett Reduction とは

えびちゃんさんの記事 がわかりやすかったです. アイデアは固定された除数 m に対して 適当な 2 冪である定数 B を使って \left \lceil \frac{B}{m} \right \rceil を前計算しておくことで,\left \lfloor \frac{N}{m} \right \rfloor\left \lfloor \frac{\left \lceil \frac{B}{m} \right \rceil \times N}{B} \right \rfloor で近似すれば除算を乗算とビット演算で置き換えられるという感じです. また, N \lt m ^ 2 の時近似の誤差が高々 1 であることが知られています.(ACLの実装も (実質)  \lt m ^ 2 を要求しています)

考えたことその1 : N \lt m ^ 2 とは限らない場合

B = qm - r \ (0 \leq r \lt m) とします.この時 \left \lceil \frac{B}{m} \right \rceil = \frac{B + r}{m} です. \left \lfloor \frac{\left \lceil \frac{B}{m} \right \rceil \times N}{B} \right \rfloor = \left \lfloor\left \lceil \frac{B}{m} \right \rceil \times  \frac{N}{B} \right \rfloor = \left \lfloor \frac{B + r}{m} \times  \frac{N}{B} \right \rfloor = \left \lfloor \frac{N}{m} +  \frac{rN}{mB} \right \rfloor で,求めたい値が \left \lfloor \frac{N}{m} \right \rfloor であることに留意すると,誤差が高々 1 になる十分条件\frac{rN}{mB} \lt \frac{mN}{mB} \leq 1 より N \leq B であればよさそうということがわかります.N \lt m ^ 2 よりも扱いやすい制約になった気がします.

考えたことその2 : 誤差を 0 にする方法

前の章と同じようなことをすると,\left \lfloor \frac{N}{m} +  \frac{rN}{mB} \right \rfloor \leq \left \lfloor \frac{N + 1}{m} \right \rfloor であればよさそうなので,こうなる条件を考えることにします.これの十分条件として \frac{rN}{mB} \lt \frac{1}{m} \Leftrightarrow rN \lt B があるので, N の最大値を  N _ { \max } として mN _ {\max} \leq B \lt 2mN _ {\max} となるように B をとることとします. この時オーバーフローが若干心配になりますが,\left \lceil \frac{B}{m} \right \rceil \times N \lt  (\frac{2mN _ {\max}}{m} +1) \times N = 2N _ {\max}^2 + N _ {\max} なので途中で出てくる数はそこまで大きくないことも確認できました.

実際この方法のほうが ACL の実装より幾分か高速な様です.(N _ {\max} = 2 ^ {64} としたとき \left \lceil \frac{B}{m} \right \rceil2 ^ {64} で切り下げ除算した結果が常に 1 なことを利用した高速化もしています.)

https://wandbox.org/permlink/2lTdpwYkMnUw3MkG

追記  N _ {\max} = 2 ^ {63} とした方がやや楽です.

https://wandbox.org/permlink/ozHZS9fb8BleVpt2

参考文献

ridiculousfish.com

rsk0315.hatenablog.com

AtCoder ARC128 C - Max Dot の双対解

問題概要

整数 N, M, S と数列 A_1, A_2, \dots, A_N が与えられる.

非負実数列  x_1, x_2 ,\dots x_N0 \leq x_1 \leq x_2 \leq \dots \leq x_N \leq M \sum _ {i=1}^{N} x_i = S を満たすとき,\sum _ {i=1}^{N} A_i x_i の最大値を求めよ .

問題へのリンク

解法

 y_i = x_i - x_{i-1}とします.(ただし  x_0 = 0とします)

 y_i を使って問題を言い換えると以下のようになります.

制約
  • y_i \geq 0
  •  \sum _ {i=1}^{N} y_i \leq M
  •  \sum _ {i=1}^{N} (N + 1 - i) \times y_i = S
最大化
  •  \underset{y_i}{\max} \sum _ {i=1}^{N} \left( \sum _ {j=i}^{N} A_j \right ) y_i

これは線形計画問題なので双対を取っても最適値は同じです.

制約
  •  \lambda _ 1 \geq 0
  •  {} ^ \forall i,  \sum _ {j=i}^{N} A_j \leq \lambda _ 1 + (N + 1 - i) \times \lambda _ 2
最小化
  •  \underset{\lambda _ 1, \lambda _ 2} {\min}M \lambda _ 1 + S \lambda _ 2

うれしいことに最適化すべき変数が 2 つになりました.

 \lambda _ 1 を固定すると  \lambda _ 2 の最小値は  O(N) で出せるので, \lambda _ 1 で三分探索して M \lambda _ 1 + S \lambda _ 2 を最小化すれば計算量は  O(N \log (M A_{\max})) です.

提出コード

おまけ 三分探索で最適解が出ることの証明

 \lambda _ 1 = p での制約を満たす  \lambda _ 2 の最小値を  f(p) と書くことにします.

ここで  (\lambda _ 1, \lambda _ 2) = (p_1, f(p_1)), (p_2, f(p_2)) は制約を満たすので,任意の 0 \leq s \leq 1 を満たす  s について  (s p_1 + (1 - s) p_2, s f(p_1) + (1 - s) f(p_2)) も制約を満たします.

 \lambda _ 1 = s p_1 + (1 - s) p_2 での  \lambda _ 2 の最小値は  f(s p_1 + (1 - s) p_2) であるので  f(s p_1 + (1 - s) p_2) \leq s f(p_1) + (1 - s) f(p_2)) であり,f(p) は凸関数です.

よって M \lambda _ 1 + S f(\lambda _ 1) も凸関数となるので,三分探索などで最小値が求められます.

supercon2021 予選

はじめに

問題

www.gsic.titech.ac.jp

githubリポジトリ

github.com

途中経過

6月2日

問題を読む。

判定問題を考える。
N^2個のマスを頂点としてそれぞれのマスからどの操作を使えるかを前計算して辺と見なせば、辺の個数は自明にO(MN^2)なのでDFSなどを使えばO(MN^2)で解ける。
ところが辺を情報を生成を適当にやると O(M^2N^2) になってTLE。
ではどうするかというと辺が使える条件を「操作の途中で盤面からはみ出さない」かつ「操作の途中で障害物にぶつからない」と言い換える。前者は各操作についてy,xの変位の\max, \minを持っておけばok。後者は操作の途中まで決めた時ちょうど障害物に当たるような場所を逆算して求めることで調べられる。 \sum |S_i| = MNなので、全体でO(MN^2)になる。
これ割と難しくないですか?diff1700ぐらいあります。

判定問題を解いたら、リポジトリを立てて、テストケースを100個生成した。
テストケースを眺めていると操作に不自然に同じ文字が続いていたのでgenerator.cのコードを読んだ。原因を理解したが、問題にはあまり関係なさそうだったのでそのあとは一切気にしなかった。

f:id:Mitarushi:20210618212254p:plain
テストケースの偏り

6月3日

判定問題のコードを書き上げた。

到達したかの判定やある場所である操作が使えるかの情報を持つための配列にvector<vector<char>>を使うかvector<vector<bool>>を使うかで迷ったが結局charにした。ちなみにあとでbool版も書いたが遅かった。

ここで頂点や辺の数を集計したファイルを作った。あとで便利だった。

github.com

このファイルを作った目的の一つにグラフの疎性の確認があって、疎なグラフでDFSをたくさん回してもいいのではという期待があったが見事に打ち砕かれた。
密なグラフということで何回もDFSを再計算するのは避け、順方向の遷移だけで済むビームサーチ方針にしようかななどと考えた。

6月4日

DFSの再計算をせずに辺を削除する方法を考えた。

最短経路木に含まれる u \rightarrow vの辺を削除するとき、到達できることの必要十分条件はは頂点(0,0)から vへのパスがあることであることは明らかなので、それをうまく管理したいと思った。
動的木のアルゴリズムを調べてみたが応用できそうに見えないうえ、実装もしたくないのであきらめた。

6月5日

蜜蜂にgit/githubの使い方を教えた。
環境構築大変だよね、わかる。

6月6日

簡単な貪欲を書いてみる。

操作の個数を0の状態から、その操作を追加したときに到達できるようになる頂点の個数が最大のものを選ぶ、というのを繰り替えすように実装した。
そのままのDFSでは複数回走らすには遅かったので、操作を追加する前の到達済みの配列を流用して、途中からDFSを回すようにすると多少早くなった。

6月7日

DFSを毎回計算して操作の良しあしを決めるのは効率が悪いので、以下の指標を導入した。
(操作を追加する前に到達できる頂点の数) + (その操作によって増える辺u \rightarrow vのうち、uにすでに到達出来ていて、v に到達できていないものの数)

前者は終了判定の時に入手可能で、後者は辺の数(\leq N^2)のオーダーで求められるので、速度が上がる。

6月11日

visualizerを書いた。

少なくともC++で書くのは苦行で、JSはよく解らないので、Pythonで書くことにした。GUIを持っているものを作るときにはOpenCVPyGameを使うのが好みなのだけど、ライブラリを使うのは避けたかったのでtkinterを使った。
visualizerを書いたことでアルゴリズムが改善したとかはなかったが、出力形式に関するバグの発見に使えたので書いてよかったと思う。 このバグを見つけた時はmainかvisualizerのバグかわからなかったが、蜜蜂先生にcheckerを書いてもらったことでバグの場所を特定できた。感謝。

f:id:Mitarushi:20210618221706p:plain
visualizer

6月12日

local checkerを書いた。

pythonだとshellの実行とか並列化がやりやすくて良き。 この時の結果は貪欲だけで 1788/100ケースだった。

6月13日

ビームサーチの実装をした。

素直に貪欲をビームサーチにした形で、貪欲の時の指標に沿って上何個かを保持する。 ビーム幅は30ぐらいしかない割に 1587/100になった、すごい。

6月14日

終点の入次数/出次数に応じて辺のスコアを変えること思いついたので、いろいろ実験を始めた。
AHCで話題に上がっていて気になっていたoptunaで入次数/出次数の重みの最適化をすることにした。実行には10並列でも半日ぐらいかかる(遅い)。

ビームサーチの重複を避けるためのzobrist hashの存在をここで思い出したので実装した。スコアが少し上がった。 \mathrm{mod} 2^{32}でのhashではあったが意外と衝突しても平気そうだった。

6月15日

探索の序盤と終盤で評価値の計算速度を縮めた。
探索の序盤、終盤はすでに到達できるuと到達できないvの候補が少ないので、少ないほうを選べば楽に全探索できる。
個の高速化でビーム幅を70に変えた。

また、時間を食うビームサーチで時間を使いきった対策に、最初に高速な貪欲出力をすることにした。

6月16日

実装も詰めに差し掛かる。
まずTLE対策に9.5秒を過ぎた時点でビーム幅を10にするようにした。ここの境界の決定には気を使っていて、ジャッジのCPUがi3-8109Uであることを予想し手元環境の約1.4倍ぐらいとみて、それに合わせて時間を決めた。

またビームサイズを増やした時のスコアの動きを見て、ビーム幅を500ぐらいに出来ないかなと考えていた。

6月17日

実質最終日。
ここで天啓に導かれ、bitsetを使うことを思いつく。

ボトルネックである操作のスコアの計算を高速化する方法を考える。
ます操作での移動量shiftと、頂点に到達したか(visited)、頂点で操作を使えるか(edge)をbitsetとして持つ。
すると求める値は popcount(((visited & edge) \lt \lt shift) & (\neg visited))となり、これはO(N^2/w)の計算量となる。

ビーム幅500とまではいかないものの128にすることができた。 1530/100

-algo.txt

初めに、答えがYESかNOであるかを求めます。
入力を読み取った後、各マスについてどの操作が使えるかを判定します。このために各操作について、操作の中でy,x座標が最大でどれだけ増える/減るかをそれぞれ求めることで範囲外に出るかどうか、また操作を「途中まで進めた時ある障害物とぶつかる」⇔「ある障害物から操作の途中まで逆算した位置から操作をする」であることを利用し判定しています。この判定にはO(N^2 M)かかります。

その後、DFSの要領で各マスに到達できるかを判定しています。グラフの頂点をマス目の各マスに見立て、グラフの有向辺をマスXから1回の操作によって移動できるマスYに対し、XからYへ貼ることでできるグラフでDFSを行っています。このグラフの頂点の数はN^2、辺の数はO(N^2 M)なのでDFSで十分高速に判定できます。

探索ではビームサーチで操作を一つずつ追加して、もし全ての頂点に達成できたなら終了するという方針にしました。終了判定にはYES/NO判定と同様にDFSを用いましたが、すでに見た頂点を探索で使いまわすことで高速に処理しています。
ビームサーチでは、ある状態に操作を追加したときの評価関数は次のものを用いました。
(追加する前に到達できる頂点数) + (追加する操作の辺のうち、端の一方がすでに到達でき、もう一方がまだ到達できないものの個数)

第一項は終了判定の結果が流用できるため、第二項を高速に処理することを考えます。
まず到達できる頂点、または到達できない頂点の個数が少ないときは、少ないほうの頂点から出る/入る辺を全て調べることで効率よく処理できます。
上のどちらでもない場合はbitsetにより高速化しました。到達できる頂点をvisitedで、その頂点で辺が使えるかをedgeとし、これをbitsetとして管理し、操作の移動量をshiftとして書くと求める値は
popcount(((edge & visited) << shift) & (\neg visited)) となり、これはワードサイズをwとしてO(N^2 / w)で処理できます。

また、ビームサーチにおいて使用する候補の操作を重複して含むことを防ぐために、zobrist hashをもちいてビームの情報を管理しています。

ビーム幅は最初は128ですが、9.5秒を超えた場合は10に設定しなおし、すぐに終了させます。

感想

終了後他の人を見て煮た方針がなくびっくりした。 最初の貪欲から劇的な改善というほどではなかったので個人的には心残りがあったが、解法ガチャはSRぐらいは引けているようで良かった。

追記(2021/06/23)

予選通過した、やったー!